English versionEN


(Ce texte n'est pas disponible en français, nous nous en excusons.)


This page provides a detailed description of the changes being made from S. Fromang's version of Zeus3d to Zeus3d+ only. It is not ment to be a user manual for Zeus. And remember: the best manual for a code is its current source code :) !

The last implementation is described in Benchmarks.



The current version can be downloaded here.


Zeus3d+ is being developed by P. Lesaffre (speaking, or rather writing...), starting from S. Fromang's version which itself had contributions from numerous authors: J. Sone, J. Hawley, C. Evans, T. Fleming, and probably many others who did not bother to sign their names in the code...



I make extensive use of CVS so that you can follow the evolution of the modifications of the code in the source code itself as CVS comments. I will provide any version you may want to get on request.


The README file provided in the main root of the archive is my personal log file and as such is probably not very readable... It consists of a day by day account of the modifications to the code.


The NOTES file is a file I wrote while first reading through the code. It may help you getting started with the organisation of the code and which routine does what, but it is now probably obsolete as I brought many strong modifications to the code.


As the code has just been benchmarked against new test cases, not all keywords are currently available.


I list here the bugs I found in the original version I got from S. Fromang (He has corrected them since). Of course you won't find here all the bugs I have introduced in this new version of Zeus since I don't know them yet ! (You will however find in the README file and in the CVS comments a few temporary bugs I introduced which are now fixed).


I reckon the call tcons should placed at the very beginning of the time step such that the shearing box boundary conditions are performed at the same time for the whole time step. (In Seb's original version it was called only at the beginning of the transport step.)

Note: This might not be a bug, but rather a matter of preference. Only benchmarks will be able to tell.

MPI shearing box boundary conditions

I noticed the mpi version and the sequential versions did not yield identical results. Neither did different ways of splitting the domain. I found two bugs (wrong loop boundaries in compute_bval_in and an error in the azymuthal shift computed in compute_y_inout_all. But fixing them was not all and I could never find the remaining bug(s).

I hence rewrote completely the shearing box boundary conditions and it now works perfectly : sequential=mpi=mpi for different splits.

tranx2 and tranx3

In the non isothermal (adiabatic) mode, calls to boundary conditions for internal energy were missing from routines tranx2 and tranx3.


This routine had two flaws: inverse of density was missing and wrong centering in one term (as noted by Seb himself). I rewrote visphys.F, compared it to the old visphys.F and benchmarked it against my own test cases so I trust it now (see VISPHYS and VISNU below).

New features

I sort the new features by their corresponding keywords which have to be "defined" in the dat/zeus3d.def file at compilation.


I implemented a new version of the boundary conditions which allows to split the computational domain in all 3 directions (azimuthal, radial and vertical). This is useful if you want to run on a large number of processors.

Note that the calls to boundary conditions are now embedded in an easy to call src/bc.F routine. For example, if you want to apply boundary conditions to variables density, velocity in x,y,z directions and energy, just use

call bc(d,v1,v2,v3,e,vmode=3)

The optional keyword vmode tells the position in the argument list of the variable which has an azimuthal velocity type of shearing box boundary condition. Other keywords are bmode to skip the ie+1 update (for b1), smode for azimuthal momentum and emode for total energy. This syntax is the same whether you compile in WITHMPI mode or if you compile the sequential version.

Note that the WITHMPI and sequential version have not been checked to be rigourously identical since the last boundary condition (keyword FBC) changes.


MPI communications are performed in serveral binary interactions rather than in two annuli (Seb's like). I did not yet notice any difference in performance between the two versions on the machine jxb.


Routines src/etot.F, src/doetot.F, src/initetot.F

This makes use of a conservative version of the energy equation. It solves for the internal energy as well as the total energy, and computes the rate of change of internal energy that needs to be applied after the MHD step. This rate is then added in the cooling step (routine src/cool.F when COOLING is on).


Routine src/transprt.F

This includes tidal potential energy in the total energy (as it should !). The complete conservative scheme has to use both ETOT and ETOTpot. Note: this is currently bugged (only radial flux gradient is applied).


This will make use of a conservative version of the momentum equations, but is currently bugged and not fully implemented...


This integrates isochore cooling with a user provided cooling function (routine src/cool.F). A linear cooling function is used if icool=1 in the pgen namelist which generates the parameters of the problem. Otherwise, the cooling function has to be supplied in routine src/smallsteps.f90.


Routine src/visphys.F

Uses a physical viscosity (advised, as the artificial viscosity is not a tensor). By default, the viscosity mu=rho nu is assumed to be a uniform constant. But the routine is written in a way such that a non-uniform viscosity can be included straightforwardly.


Routine visphys.F

The viscosity nu=mu/rho is a constant. This is the most physical choice and this keyword should be on for compressible applications.


Routine mhd/resist.F

Uses an exponential (Saha-like) temperature dependence for the resistivity.


Routine mhd/resist.F.

Uses a 3/2 power (fully ionised gas) temperature dependence for the resistivity.


Forces the mass density to be equal to 1.


Discards the pressure gradient term and the call to mhd/magpres.F (which is not the magnetic pressure gradient term... caution !).


I recently found new sets of benchmarks for shearing box MHD flows that include resistivity, viscosity and cooling (Lesaffre & Balbus, 2007 MNRAS submitted). The code has been benchmarked against these test cases.

In particular, these solutions include modes with non radial wavenumbers. The code is unable to model them accurately unless the transport step is modified to advect separately the mean shear and the perturbed part of the flow. This is the last implementation I did, for which ETOT and ETOTpot need to be updated.

Even more recently, I found analytical solutions for non zero azimuthal wavenumbers and am about to test the implementation of the shearing box boundary conditions in an MHD case.


An example of a snapshot of a simulation ran on jxb with Zeus3d+. It uses the temperature-dependent resistivity of a fully ionised gas, linear cooling and a cubic box of 4 scale heights aside (no stratification). The plotting routine is pro/zload.pro.


Voir aussi dans la même rubrique