The Pencil Code User Meeting 2008 Official Blog

Large-scale dynamos in turbulent convection with shear - Petri

2008-08-22T00:43:00.000-07:00

Good introductory link to the Solar Dynamo: scholarpedia article

http://www.scholarpedia.org/article/Solar_dynamo

Related publication: http://arxiv.org/abs/0806.0375

Need shear to avoid catastrophic quenching. Dual role of shear: generation and helicity loss.

New Users

2008-08-21T08:01:00.001-07:00

Andrej

Relative Velocities and Collision Speeds of Particles in PPDs

Voelk et al. (1980) -> long standing result for formulas for relative velocities of two particles with different stopping times.

Particles of different stopping times drift differently and meet at higher velocities. Growth is hindered away from the diagonal.

First Science: shearing boxes with the Pencil code, forced hydro turbulence, MRI.

Results of collisional speeds.

Martin

Magnetic field dissipation in a jet

Originally for Poynting-flux dominated flows. If the field decays faster than 1/r it produces acceleration even beyond the Alfven radius (The radius at which magnetic pressure equals the thermal pressure, nearly equal to the magnetospheric radius).

Naturally produces radiation (gyromagnetic, ciclotron, sincrotron)

MHD in expanding coordinates. Interplay jet/box. IC: steady state wind 1D (1.5D solution)

BC: the biggest possible in the radial direction... (?)

Raphael

Chemical Symmetry breaking. One of the processes that was quite important for the emergence of life. Applications to astrobiology (explains why there is a "astrobiology_data" subroutine in chemistry.f90)

Plan to use the Pencil code in this problem: Inhomogeneities in the energy source. Put some localized sources of energy. Chemicals coming from crust, or light coming from above.

Excitation of Inertial-Acoustic Waves in PPDisks (Tobi)

2008-08-21T05:20:00.000-07:00

Properties of Keplerian flows that develop in the accretion disk. Shear, angular momentum transport, etc, check Balbus & Hawley 1998 for a thorough review.

Propagation of waves in the MRI turbulence. High-time resolution needed to see the waves. Some people have seen it before but mistook it as numerical artifacts.

The waves are trailing.

Radial wavenumber depends linearly in time. The fourier modes swing successively from leading (kx/ky<0)>0).

The power shows shear in k-space.

Inertial acoustic wave equation

d2p/dt2+[k^2c^2+Omega^2]p=-c^2 i kx(t) ksi + (q-2)Omega*Nx + dNy/dt

ksi=vortensity

N=div(T), T=stress tensor (nonlinear terms)

WBKJ analysis reveals that there is a preferred wavenumber for which the amplitude of the waves is maximum. The preferred wavenumber is very low.

The steady damping of the waves seen towards high wavenumbers may be due to turbulent diffusion.

A Monte-Carlo coagulation-fragmentation model for the Pencil Code (Andras)

2008-08-21T02:50:00.000-07:00

Link to the talk (pdf)

Try overcoming the meter sized barrier by coagulation instead of self-gravity

1g -> 1e12 monomers inside a particle. Not feasible with current computers.

Need something more statistical.

Critics to Ander's nature paper: starts with particles that are already too big (40-60-80)cm, and without dust physics (coagulation-fragmentation; CF). Can planet formation work with smaller particles and detailed CF models?

Method. Add virtual particles. Real particles grow by colliding with virtual particles. Real particles do not meet, the rate of collision real-real is much smaller than real-virtual. The model is just real-virtual. There are no virtual-virtual collision either.

Test: Smoluchowski equation that describes coagulation. The test reproduces the mass x number density distribution. For comparison, see Ormel et al. 2008, Wetherill 1990.

Important energies: Eroll-> Energy needed to rotate a monomer by 90 degrees.

Porosity model:

Ecoll < 5Eroll (Hit and stick)

If 5Eroll < Ecoll < Efrag (compaction)

Fragmentation always lead to a monomer (limitation)

With porosity, the radius gets bigger by a factor 2 compared to just using turbulence + brownian
motion.

For the model showed, the grains reach complete compactness at 0.1cm size. The MMSN reaches compactness just at m size.

The code seems fast. 1 million timesteps takes minutes in a single processor. Seems no overhead for the pencil code.

Simulations of the geodynamo (Graeme)

2008-08-20T04:42:00.000-07:00

Link to talk (pdf)

Thermal convection in a rotating spherical shell of weakly compressible fluid (motivated by geodynamo)

Cartesian domain (equivalent to lsphere_in_a_box)

Damping outside boundaries (for velocity and entropy)

Standard MHD with Laplacian diffusion (what's the Reynolds number?)

Polytropic ideal gas with a background state (density and temperature)

Non-equidistant grid outside the shell (step-linear)

Full sphere and half-sphere runs. Mag field boundaries: Btan=0, d(Bn)/dn=0. (Has to be adapted for spherical coordinates.)

Comparison with Kageyama & Sato, Phys Plasmsas, 1995 (KS95). Early results (McMillan and Sarson, PEPI, 2005) are not correct due to scaling glitch. (cp=1 and R=2/5 were hard-coded in the Pencil Code at that time, before the EquationOfState module was developed.)

Weak dependence of vorticity on the polar direction for non-magnetic convection.

When putting the magnetic field with small magnetic diffusivity (high Rm), dynamo action is obtained. At higher Rm, the simple convection roll pattern is disrupted.

There is a dynamo benchmark in the geodynamo community (Christensen et al., PEPI, 2001), but just for the Boussinesq case . Unfortunately this is not a stable state for a compressible gas.

Qualitatively speaking, the results are similar to another dynamo benchmark.

Future work:

tidy up and publish (good plan)
further investigate reversing solutions

Boussinesq and anelastic fluids
spherical coordinates
insulating boundary conditions
make Coriolis force implicit

Local box simulations of geodynamo

Testfield method (Axel)

2008-08-20T02:10:00.001-07:00

Links to the talk ( ppt, pdf)
Prolific output in 2008 after quiescent period after 2005

In many astrophysical systems a large scale magnetic field is generated, but the mechanism for their formation is still elusive.

Helical turbulence: Homogeneous and isotropic. However, all helices have the same handiness.

Large scale magnetic fields usually appear in spacetime plots of fields as a finite amplitude field after dimensional average.

Source of turbulence in forced dynamos: a body force. Non-natural, closest analog being supernova explosions.

Convection (with shear) and MRI are two sources of turbulence that can generate large magnetic fields. Structure visible in the z-direction.

Low Prandlt number dynamos (Pm=visc/resist):

Helical turbulence: inject energy at intermediate wavelengths. Energy cascades to the right (smaller wavelengths) and to the left (larger scales). This is due to helicity. At non-helical turbulence at other Pradtl numbers, one can excite the dynamos at progressively lower Reynolds numbers as the Prandlt number increases. In the Sun the Prandt number is 1e-5. So, very large Reynolds numbers are needed.

With helical turbulence the dynamo is still excited at low Prandlt numbers. The only difference is that the viscous range moves towards larger scales (still smaller than the scale where energy in being forcely supplied). The inertial range shrinks. So, astrophysics don't have any problems with small Prandlt numbers. The only thing we need is stratification to drive helicity loss.

Testfield equations:
perturb the original induction equation and separate them in mean field and fluctuations. The eq. for the fluctuations cannot induce a large scale field.

Validation: Roberts flow
Non-linear alpha and eta tensors
Turbulent diffusion as a function of Reynolds number. Quenching goes as a factor 5. Turbulent diffusion is not catastrophically quenched. Several consistency checks.

Turbulent Combustion (Natalia)

2008-08-20T01:02:00.000-07:00

Link to the talk (pdf)

Outline of theory, challenges and preliminary results.

Many species, every single one of them contributes to the fluid equations.

The pressure is

P=m/mu * R * T

1/mu=Sum_k (Y_k/m_k)

And this now is a pencil.

Initial condition: air mixture at atmospheric pressure and high temperature (1200K)

Reaction mechanisms (8 equations) are used instead of the 19 of the Li mech Nils talked about. They should give the same results, but the code is not yet at that stage.

Chemistry Modules (Nils)

2008-08-20T00:25:00.000-07:00

Link to the talk (pdf)

Rationale: 90% of the World energy comes from combustion

Challenges of combustion:
Pollution (NOx and C02)

NOx: produces acid rain
CO2: greenhouse gas

Efficiency:

Safety: Combustion of hydrogen. Hydrogen's energy density is not high, one has to compress hydrogen quite a lot for combustion. Highly compressed hydrogen is quite dangerous.

Challenges:
Rapid mixing of air and fuel
Mixing quality -> emissions
Air at high p and T=> auto-ignition
High flame speed -> flashback
Fuel and operational flexibility

Needs:
Better fundamental understanding of physics/kinetics at relevant conditions
Experimental data at the parameter range
Simulated models to predict outcomes.

Hydrogen is a hard fuel to handle (high flame speed, very inflamable, etc)

Mixing in turbulent reactions takes place at the smaller scales of turbulence
How to resolve them in numerical simulations?
And how to couple with the combustion reactions?

Need to resolve smaller than the Kolmogorov length scale.

Reations: example -> hydrogen and oxygen making steam (simplest reaction)

19 subreactions have to happen to produce steam

Methane + plus takes 500 rections with more than 50 different species. Burning longer hydrocarbons takes even more reactions and catalysts.

Two basic geometries:
jet in co-flow (mixed, pre-mixed0

Anchored channel flow
premixed

Solid walls produce turbulence, so high resolution is needed near walls if one wants to have them. Avoiadable.

Boundaries:
Inner: turbulent non-reflecting inlet (NSCBC)
Outer: non-reflective outlet conditions (NSCBC)

Navier-Stokes Characteristic Boundary COnditions (NSCBS), Poinsot & Lele (1992)

Set the time derivative at the boundary (required a hook in pde)

The subroutine rewrites the boundary condition before the runge-kutta adds to the f-array

Timestep controlled by the fastest reaction.

Sample: H2 combustion with Lithium mechanism (PhD thesis of Andera Gruber)
resolution: 630x200x280

Box size: 2cm! Not applicable to industrial scales yet.

250 k CPU hours.

Tuesday Discussion

2008-08-19T07:43:00.000-07:00

Usernames in svn:

We should have usernames as close as possible to first initial+last name

Backwards compatibility: the older usernames should stay for consistency with the thousands of line of code through the year that are stored in the cvs repository.

Everyone with check-in permission should have an entry at the developers file. Otherwise, as the email notification does not show everyone, we sometimes cannot contact someone that made a suspicious change in the code.

Notification:

There are 152 people with check in permission, but a lot of these people not active on the code anymore. We could remove their writing permission and shrink the list to 30-40 people only.

There are people who used the code for some time, then they moved to another project and come back to the code after an extended period of time. Should these people be bothered by tons of commit messages? Graeme says this is not a problem, because in a worst case scenario he just deletes them. Andre tends to focus on the modules he works on. Most of the active people read the messages (no matter how many) and delete them. Axel mentioned that it would be good if people could choose which modules they want to receive messages about. Wolfgang described it as a bookkeeping nightmare.

Axel's suggestion

SVN notification
= SVN permission
= phone list
= bug tracker (trac)

A person in any of these "lists" should get notification. They lose their notification when logging out of the three lists.

Discussion: Putting the Pencil code in SourceForge can turn the notification process relatively painless. Should we do it now or leave it for the pencil code meeting 2009? Another option is Google Code. They even have a space for blogs, so this blog could be moved there as well.

Postprocessing with Python (Boris)

2008-08-19T04:58:00.000-07:00

Overview:
-----------------

- What is Python and why use it?
- Installation procedure

What is Python?

Python was created in 1991 by Guido van Rossum (Amsterdam). Named after Monty Python's flying circus.

Goals:
----------

- an easy and intuitive language
- opens source
- plain English
- suitable for every tasks, short development times

Python v. 3.0 in 2009.

Why Python?
---------------------

- it's free
- quite easy to use; object-oriented; highly modular, etc.
- faster than IDL, and even parallel (MPI)

Python makes extensive use of modules, and they all have to be loaded.

The main Python website is www.python.org. Another important website is www.scipy.org. scipy is a collection of scientific modules for Python.

How to install Python?
--------------------------------------

Required:
- Python 2.5: the engine
- numpy: the scientific computing package
- scipy: modules for integrating ODEs, optimizing functions, etc.
- matplotlib: MATLAB-inspired mostly-2D plotting modules

numpy may be moving into scipy.

Optional:
- ipython: interactive shell
- basemap: map projection
- Pypar: parallel Python (interface with MPI libraries)
- PyQt: to do Qt-like widgets very easily
- MayaVi: 3D plotting

Download the code at http://www.python.org/download.

Should you download sources or binary packages? For Linux, Windows & Mac binaries are provided (e.g. using apt-get for Linux or Fink for Mac). For Mac there is the website www.pythonmac.org/packages/py25-fat. Can also install Scipy Superpack for OS X.

Jeff Oishi was the first to code Python scripts for postprocessing with the Pencil Code.

Need to set $PYTHONPATH.

if ($?PYTHONPATH) then
setenv PYTHONPATH "${PYTHONPATH}:${PENCIL_HOME}/numpy"
else
setenv PYTHONPATH "${PENCIL_HOME}/numpy"
endif

Postprocessing with Python:
------------------------------------------------

Start with 'import pencil as pc'

Can either do pc.read_ts() or a=pc.read_ts(). Then 'import matplotlib as P'. Then 'plot(a.t,a.urms)'

Most of the functionality of the idl pc_read modules (including magic) is implemented in the python counterparts.

Be careful: Python's arrays are ordered like f[nvar,mz,my,mz].
REVERSED COMPARED TO FORTRAN OR IDL

For 2D plots, use imshow.
For 3D plots one can use MayaVi. http://code.enthought.com/projects/mayavi

Migrating from IDL to Python:
---------------------------------------------------

There are good webpages with overview of IDL commands and the corresponding Python commands. Space Telescope has a nice page on that.

Other useful page

https://mathesaurus.sourceforge.net/idl-numpy.html

For lazy people there is the i2py converter (http://code.google.com/p/i2py).

Some tricks when using Python:
------------------------------------------------------

Load modules by default in ~/.ipython/ipythonrc.

Important only what you need, otherwise it takes too long to load Python.

Launch ipython with the '-pylab' option to call directly plotting commands.

Accelerate reading of snapshots by passing arguments.

Accelerate graphics by using a handle.

Take advantage of class and objects.

Parallel Python:
---------------------------

Can run Python with mpirun and do postprocessing faster.

Widgets using Q Designer + PyQt:
----------------------------------------------------------

Write script in Qt and convert to Python. Simpler to do widgets than in IDL.

Conclusion:
--------------------

- Python can do a very good job in the Pencil Code postprocessing
- Its use is rapidly increasing in astrophysics
- More in the Pencil Code philosophy (i.e. under GPL)

Weaknesses:
- the actual Python subroutines must be rewritten in a more object oriented form.
- need to reorganize the Python tree, in something like f90/pencil-code/python.
- calling Fortran or C subroutines to increase the speed?

Tutorial: "Changing from CVS to SVN" (Tobi)

2008-08-19T01:57:00.000-07:00

CVS is outdated. Move to SVN.

SVN itself is not so modern either, it is being replaced by Distributed Version Control.

SVN -> very similar to CVS. There are actually several "SVN for CVS users" guides in the web.

Important:

CVS repository moved to f90/pencil-inactive

and will not allow commits anymore.

Check out the code with svn:

svn co https://ajohan@svn.nordita.org:/svn/pencil-code/trunk pencil-code

Like cvs, update with

svn update

the latest revision is revision 9758.

For a short list of svn commands, do "svn help". For a comprehensive guide, google "svn book".

svn help also works for a command. For instance, type

svn update --help

for a list of subcommands of svn update. An important difference as compared to cvs is that the svn update is recursive. Use

svn update -N

for updating only the current directory and not the subtrees.

Commiting files. Same as with cvs:

svn ci -m "comment" file

We could show more ten commands or so, but it is better to start using it. The questions come with the experience.