5.4.7 Hydrodynamics Control (HD and MHD)

Many parameters in the parameter file that are valid for the HD module also apply for the MHD module. Only a few parameters for the hydrodynamics control (character hdscheme: Sect. 5.4.7, character reconstruction: Sect. 5.4.7, real c_slopered: Sect. 5.4.7) are slightly modified or extended for the MHD module.

Table 11: List of reasonable combinations of options for HD and MHD solver. The columns are the type of the (M)HD scheme, see Sect. 5.4.7, the time-integration scheme, see Sect. 5.4.7, the (directional) splitting method, see Sect. 5.4.7, the maximum Courant number (where $D$ is the dimension of the computational box), see Sect. 5.4.15, the maximum sound-Courant number (where $D$ is the dimension of the computational box) only specified for the case where it is actually a new constraint in addition to the normal Courant condition, see Sect. 5.4.15, and a comment.

hdScheme	hdTimeIntegrationScheme	hdSplit	$C_\mathrm{Cour}$	$C_\mathrm{SndCour}$	comment
Roe	Single	123	$1$		old default
Roe	Single	unsplit	$1/D$
Roe	Single	CTU	$1$	$1/2$ for $D$=$3$	recommended
Roe	RungeKutta2	unsplit	$1/D$
Roe	RungeKutta3	unsplit	$1/D$
HLLMHD	Hancock	123
HLLMHD	Hancock	unsplit
HLLMHD	RungeKutta2	unsplit
HLLMHD	RungeKutta3	unsplit

• character hdscheme:
  With this parameter the type of the hydrodynamics scheme can be specified as in
  character hdscheme f=A80 b=80 n='Hydrodynamics scheme' &
    c0='Roe (approximate Riemann solver of Roe type)' &
    c1='RoeMagKin (Roe solver + kinetic magnetic field transport)' &
    c2='RoeMHD/HLLMHD (MHD solver)' &
    c3='None (skip hydrodynamics step entirely)'
  Roe
  
  Possible values are
  • None: The hydrodynamics step is skipped entirely (for test purposes). Note that in this case some initializations necessary for the generation of the mean file are omitted, too.
  • Roe: (default) The standard approximate Riemann solver of Roe-type is activated. This value will be chosen in most cases.
  • RoeMagKin: The standard Roe solver is extended to transport passively a magnetic field. This is a test implementation to check, if the general magnetic field handling works. For proper MHD simulations use HLLMHD instead.
  • RoeMHD: Use the MHD Roe solver.
  • HLLMHD: Use the MHD HLL solver (the current recommend value for MHD simulations).

• real c_hydlowmachcsfactor:
  For very low Mach numbers the over-sensitive sound-wave detection of the Roe solver (for 2D and 3D cases) can be suppressed by choosing a non-zero value for real c_hydlowmachcsfactor. This is in fact a tiny modification of the Roe solver but changes its behavior significantly and requires (for C_hydLowMachcsFactor>0.0) hdsplit=unsplit and hdtimeintegrationscheme=RungeKutta2 or hdtimeintegrationscheme=RungeKutta3. The value can be set e.g. with
  real c_hydlowmachcsfactor f=E15.8 b=4 &
      n='Low-Mach-Number sound-speed reduction parameter' &
      c0='<=0.0: off (old settings, recommended default)' &
      c1='1.0: small value; 2.0: reasonable value; 5.0: large value'
  2.0
  
  Typical choices are
  • 0.0: The extra reduction factor is not applied. This is the standard Roe scheme. It is still the recommended default, excepd for very-low-Mach-number flows.
  • 1.0: Smallest meaningful value to switch on the reduction.
  • 2.0: Reasonable value
  • 5.0: Rather large value
  The (recommeded) default is not to use this reduction and leave the parameter at 0.0. This parameter is not recognized by the MHD module.

• character hdtimeintegrationscheme:
  With this parameter the type of the time-integration scheme can be specified as in
  character hdtimeintegrationscheme f=A80 b=80 n='Time-integration scheme' &
    c0='Single/Hancock/Euler1/RungeKutta2/RungeKutta3'
  Single
  
  Possible values are
  • Single: (default for HD Roe) Compute all updates in one single step. Terms higher order in time are provided by the Roe solver depending on the order of the reconstruction scheme.
  • Hancock: (default for MHD HLLE) Terms of second order in time are approximated by a Hancock predictor step.
  • Euler1: (only for MHD) A single Euler step is used, which is only first order in time.
  • RungeKutta2: A two-step (predictor-corrector) Runge-Kutta step is used.
  • RungeKutta3: A three-step Runge-Kutta step is used.
  For possible combinations see 11.

• character hdsplit:
  With this parameter the type of the hydrodynamics operator (directional) splitting scheme can be specified as in
  character hdsplit f=A80 b=80 n='Hydrodynamics directional splitting scheme' &
    c0='123: directional splitting (default)' &
    c1='unsplit: unsplit (direct) operator' &
    c2='CTU: Colellas Corner Transport Upwind scheme (1990)'
  123
  
  Possible values are
  • 123: (default) The standard directional splitting is activated, where the 1D operators for the individual directions are applied in the given order. So far, only the order 123 is possible.
  • unsplit: The standard solver is applied. However, the changes from the individual steps are computed from (and applied to) the same model configuration. The result is an unsplit scheme.
  • CTU: Use Colella's "Corner Transport Upwind" scheme ('CTU', 1990). It first computesand stores the changes from hydrodynamics in x1, x2, and x3 direction. Then (half of) the contributions from the steps in x1 and x2 direction are used as sort of a predictor step for the full ('corrector') step in x3 direction (and accordingly for the two other directions). At the end, the three corrector-step contributions are added. That means, that two hydrodynamics steps have to be performed in each direction with some overhead for copying and adding the contributions. Still, together with e.g., reconstruction=FRweno, there is a dramatic improvement for convection at small Mach numbers compared to 123, vanLeer.
  This parameter is recognized by hdscheme=Roe and hdscheme=HLLMHD, although hdscheme=HLLMHD does not know hdsplit=CTU.

• character hdcheckflux:
  With this parameter the CheckFlux routine of the HLLMHD solver can be (des)activated, as e.g. in
  character hdcheckflux f=A80 b=80 n='Switch to activate checking of MHD fluxes' &
    c0='on/off'
  off
  
  This parameter is only recognized by the HLLMHD scheme. Leaving it on might be safer but is definitely (slightly) slower.

• character reconstruction:
  This parameter determines the order and ``aggressiveness'' of the reconstruction scheme with e.g.
  character reconstruction f=A80 b=80 n='Reconstruction method' &
    c0=Constant c1=Minmod/VanLeer/Superbee c2=PP/PPmimo &
    c3=FRmimo,FRmono(2),FRweno(2),FRcont
  FRweno
  
  Possible values are
  • Constant: The run of the partial waves inside the cells is assumed to be constant. A highly dissipative first order scheme results. This values will usually only be used for test (or comparison) purposes.
  • Minmod: Chooses the smallest slope which still results in a second order scheme. It is the most diffusive (and most stable) one in this class.
  • VanLeer: (default) The recommended first-order reconstruction method.
  • Superbee: The ``most aggressive'' first-order reconstruction method. It results in the steepest shocks, which works well in some test cases but tends to turn smooth gradients into steps and might be to difficult for the radiation transport module to handle.
  • PPmimo: MinMod reconstruction based on boundary values from PP
  • PP: Chooses the piecewise parabolic reconstruction of the PPM scheme (``Piecewise Parabolic Method'', [Colella and Woodward (1984)]). Results in 3rd order accuracy for the advection.
  • FRmimo: MinMod reconstruction based on boundary values from FRmono
  • FRmono: 2nd-order monotonic reconstruction derived from PP, WENO, and MinMod: ``Frankenstein's Method''
  • FRmono2: 2nd-order monotonic reconstruction derived from PP, WENO, and MinMod. It is identical to FRmono if real c_slopered=0.0 (the usually adopted default choice, see Section. 5.4.7). However, for real c_slopered$>$0.0, it uses Constant reconstruction as fallback, instead of
  • Minmod as
  • FRmono does. This is the recommended version for low-resolution red-supergiant models.
  • FRweno: 2nd-order non-monotonic version of FR
  • FRweno2: 2nd-order non-monotonic version of FR. It is identical to FRweno if real c_slopered=0.0 (the usually adopted default choice, see Section. 5.4.7). However, for real c_slopered$>$0.0, it uses Constant reconstruction as fallback, instead of
  • Minmod as
  • FRweno does. This is the recommended version for low-resolution red-supergiant models.
  • FRcont: 2nd-order continuous version of FR, unstable except for low Mach numbers.
  Usually, the VanLeer reconstruction has been a good choice. If a more stable (and diffusive) scheme is needed, take Minmod. However, the variants of the new 2nd-order ``Frankenstein's Method'' are a better choice, if the additional computational cost is acceptable. The PP reconstruction gives the highest accuracy. However, there is the possibility that it produces somewhat ``noisy'' models with small wiggles e.g. in the velocity for extreme cases - solar models look fine with PP.
  The options Constant, SuperBee, and FRcont are just for testing purposes: they have either too much, sometimes the wrong sign of, or too little diffusion, respectively.
  The order of decreasing diffusivity and increasing accuracy is roughly (schemes on a line are similar; schemes in brackets are for testing only):
  (Constant)
  MinMod
FRmimo, PPmimo, vanLeer
FRmono
FRweno, (Superbee,) PP
  (FRcont).
  However, there are differences how the increase in accuracy is achieved: vanLeer uses (potentially dangerous) anti-diffusion (in contrast to MinMod) but needs only a 3-point stencil, whereas PPmimo and FRmimo don't have anti-diffusion but need a (more expensive) 5-point stencil. FRweno achieves actually a 2nd-order-accurate reconstruction also near extrema for reasonably resolved structures (e.g., sufficiently long wavelengths of sine waves). I.e., it gives up strict monotonicity, whereas PP is strictly monotonic and achieves low dissipation from avoiding jumps at cell boundaries wherever possible. However, it chops off (flattens) extrema like all TVD schemes and it shows Gibbs wiggles near a jump that is not surrounded by constant plateaus but lies in a slope. The new 2nd-order scheme FRmono is designed to suppress these wiggles, too. The highest-resolution scheme is PP. It is implemented in the MHD module, too.
  By specifying e.g., reconstruction=FRweno PP it is possible to use the FRweno scheme for the hydrodynamics scheme as such and the monotonic PP only for the advection of additional quc quantities.

• real c_slopered:
  A new extra stabilization mechanism can be activated that reduces the slope and flattens the reconstruction function in case of a strong density contrast ( $\rho_\mathrm{one side}/\rho_\mathrm{other side} > 2$ between neighboring cells). This value can be set e.g. with
  real c_slopered f=E15.8 b=4 &
      n='Slope reduction parameter in case of strong density contrast' u=1 &
      c0='0.00: off (default), 0.02: reasonable value, 0.10: large value'
  0.02
  
  Typical choices are
  • 0.0: Slope reduction switched off. Original reconstruction is used.
  • 0.02: Moderate slope reduction in case of large density jumps.
  • 0.10: More pronounced slope reduction in case of strong density contrast.
  The default is not to use this trick and leave the parameter at 0.0. However, if numerical problems occur close to strong density jumps one can try to cure them by setting it e.g., to 0.02 or so. One could even think off other criteria (e.g., low temperatures, very high Mach numbers,...) to activate this slope reduction. That could be implemented if if would deem beneficial.

• real c_reccontshift:
  The second-order reconstruction schemes except FRcont have the option to make the reconstruction even more continuous in the case of very low Mach numbers, controlled with C_RecContShift and C_RecContSteep. The value can be set e.g. with
  real c_reccontshift f=E15.8 b=4 &
      n='Enhanced reconstruction continuation parameter 1 for shifting' &
      c0='0.0: off; 1.0: reasonable value'
  0.2
  
  Typical choices are
  • 0.0: No attempt is made to make the reconstruction more continuous. The original reconstruction is used.
  • 0.1: Make the transition to a smoother reconstruction at lower Mach numbers.
  • 1.0: Reasonable value.
  • 10.0: Make the transition to a smoother reconstruction at higher Mach numbers.
  While the resulting scheme works for solar simulations, the parameters are really of benefit only in the case of brown dwarfs or exoplanets or the deep layers of M dwarfs.

• real c_reccontsteep:
  The second-order reconstruction schemes except FRcont have the option to make the reconstruction even more continuous in the case of very low Mach numbers, controlled with C_RecContShift and C_RecContSteep. The value can be set e.g. with
  real c_reccontsteep f=E15.8 b=4 &
      n='Enhanced reconstruction continuation parameter 1 for steepening' &
      c0='0.0: off; 0.1: reasonable value'
  0.5
  
  Typical choices are
  • 0.0: No attempt is made to make the reconstruction more continuous. The original reconstruction is used.
  • 0.01: Make the transition to a smoother reconstruction steeper (at even lower Mach numbers)
  • 0.1: Reasonable value.
  • 1.0: Make the transition to a smoother reconstruction smoother
  While the resulting scheme works for solar simulations, the parameters are really of benefit only in the case of brown dwarfs or exoplanets or the deep layers of M dwarfs.

• character hdtransvelomode:
  This parameter controls the treatment of advection of transverse velocities. It can be set e.g. with:
  character hdtransvelomode  f=A80 b=80 &
    n='Mode for transverse velocity advection' &
    c0='Normal(default)/CA1/CA2/CA3/CA4'
  CA4
  
  The choices are
  • Normal, ' ': Compute flux of transversal momentum with standard Roe treatment (old, default)
  • CA1: Compute flux of transversal momentum with consistent advection, version 1. Use mass flux and upwind density to determine flow speed. Actually recompute transversal velocities.
  • CA2: Compute flux of transversal momentum with consistent advection, version 2. Use mass flux and Roe-mean density ($\tilde{\rho}$) to determine flow speed. Actually recompute transversal velocities.
  • CA3: Compute flux of transversal momentum with consistent advection, version 3. Use mass flux and upwind density to determine flow speed. Compute transversal momentum fluxes from mass flux and diffusive momentum flux.
  • CA4: Compute flux of transversal momentum with consistent advection, version 4. Use mass flux and Roe-mean density ($\tilde{\rho}$) to determine flow speed. Compute transversal momentum fluxes from mass flux and diffusive momentum flux, (recommended).
  For instance for a 2D Sedov blast wave problem, the new modes give less artifacts close to the axes than the old mode Normal. The four new advection types are made of to choices for the way to compute the effective flow speed across a cell boundary (different density values can be used) and two ways to compose the final transversal momentum fluxes. Version CA4 is recommeded. This parameter is not recognized by the MHD module.

• character hdentropywavemode:
  This parameter controls the computation of the Roe correction flux for the entropy wave. It can be set e.g. with:
  character hdentropywavemode f=A80 b=80 &
    n='Mode for entropy wave type' &
    c0='Normal/Entropy1/Entropy2/Entropy3M/Entropy3Ms/Entropy3gMps(recommended)'
  Entropy3gMps
  
  The choices are
  • Normal, ' ': Compute entropy-wave flux with standard Roe treatment from pressure and density, (old, default).
  • Entropy1: Compute entropy-wave flux from entropy. This only works well for small Mach numbers.
  • Entropy2: Compute entropy-wave flux from entropy and pressure.
  • Entropy3M: Compute entropy-wave flux from entropy, pressure, and density. In fact, an cases Normal and Entropy1 are interpolated according to the Mach number.
  • Entropy3Ms: Compute entropy-wave flux from entropy, pressure, and density. In fact, an cases Normal and Entropy1 are interpolated according to the Mach number and entropy fluctuations.
  • Entropy3gMps: Compute entropy-wave flux from entropy, pressure, and density. In fact, an cases Normal and Entropy1 are interpolated according to the Mach number and fluctuations in pressure (directly and after subtracting the hydrostatic pressure correction) and entropy. This case is recommended.
  The Normal computation with the standard Roe procedure is optimized for shocks and works also for small fluctuations in density and pressure. However, in the case of a nearly adiabatic stratification, the density and pressure fluctuations between adjacent depth points induce some spurious fluxes. These are avoided by reconstructing the entropy directly (Entropy1). However, that does not work well for shocks. Therefore, an interpolation between the standard reconstruction based on pressure and density and another one based on entropy is recommended (Entropy3gMps). This parameter is not recognized by the MHD module.

• character hdenthalpyavgmode:
  This parameter controls the type of averaging used to get from the cell-centered enthalphy (flux) the boundary-centered enthalphy (flux). To this flux, the Roe corrections are added. It can be set e.g. with:
  character hdenthalpyavgmode f=A80 b=80 &
    n='Mode for enthalpy averaging' &
    c0='Normal(old default)/CA1/CA1g/CA1gMps(recommeded)/CA1M/CA1Ms/CA1p/CA1s' &
    c1='CA2/CA2g/CA2M/CA2p/CA2s/CA3/CA4'
  CA1gMps
  
  The choices are
  • Normal, ' ': Use the standard Roe-averaging of the enthalpy + kinetic-energy fluxes, (old, default).
  • CA1: Average the enthalpy and multiply with the average mass flux. Average the kinetic energy fluxes as in the standard Roe case.
  • CA1g: Use the consistent-advection centering for the average enthalpy flux as in CA1 for small fluctuations. Switch smoothly to the Normal Roe averaging for large non-hydrostatic pressure fluctuations.
  • CA1gMps: Use the consistent-advection centering for the average enthalpy flux as in CA1 for small Mach numbers. Switch smoothly to the Normal Roe averaging for large Mach numbers or large fluctuations fluctuations in pressure (directly and after subtracting the hydrostatic pressure correction) and entropy. This case is recommended.
  • CA1M: Use the consistent-advection centering for the average enthalpy flux as in CA1 for small Mach numbers. Switch smoothly to the Normal Roe averaging for large Mach numbers.
  • CA1Ms: Use the consistent-advection centering for the average enthalpy flux as in CA1 for small Mach numbers. Switch smoothly to the Normal Roe averaging for large Mach numbers or large entropy fluctuations.
  • CA1p: Use the consistent-advection centering for the average enthalpy flux as in CA1 for small pressure fluctuations. Switch smoothly to the Normal Roe averaging for large pressure fluctuations.
  • CA1s: Use the consistent-advection centering for the average enthalpy flux as in CA1 for small entropy fluctuations. Switch smoothly to the Normal Roe averaging for large entropy fluctuations.
  • CA2: Average the enthalpy plus kinetic energy and multiply with the average mass flux.
  • CA2g: Use the consistent-advection centering for the average enthalpy plus kinetic energy flux as in CA2 for small fluctuations. Switch smoothly to the Normal Roe averaging for large non-hydrostatic pressure fluctuations.
  • CA2M: Use the consistent-advection centering for the average enthalpy plus kinetic energy flux as in CA2 for Mach numbers. Switch smoothly to the Normal Roe averaging for large Mach numbers.
  • CA2p: Use the consistent-advection centering for the average enthalpy plus kinetic energy flux as in CA2 for small pressure fluctuations. Switch smoothly to the Normal Roe averaging for large pressure fluctuations.
  • CA2s: Use the consistent-advection centering for the average enthalpy plus kinetic energy flux as in CA2 for small entropy fluctuations. Switch smoothly to the Normal Roe averaging for large entropy fluctuations.
  • CA3: Take the Roe tilde average enthalpy plus kinetic energy and multiply with the computed mass flux. This is a more experimental version and not generally recommended.
  • CA4: Averaging the enthalpy and multiply with the computed mass flux. This is a more experimental version and not generally recommended.
  Here, 'CA1' and 'CA2' behave in a similar way, whereas'CA3' and 'CA4' are essentially just for tests. I prefer the CA1 variant for the time being. However, both CA1 and CA2 work only well for small fluctuations and Mach numbers. Therefore, one of the modified schemes should be used. Again, they are rather similar even if the implemented creteria are rather different. While most of the versions work for ``ordinary'' main-sequence models, the comparably pooly resolved red-supergiant models (with large entropy step and Mach numbers) require a combination of the criteria above. Therefore, CA1gMps is recommended (for all type of models in general). This parameter is not recognized by the MHD module.

• real c_hydpredfactor:
  The "hydrostatic pressure correction" terms in the two acoustic waves are always present. However, for the entropy wave and the "ionization wave" it is not quite clear if the terms should be there (the classical default) or if they should be set to zero. The value of a factor in front of these terms can be set e.g. with
  real c_hydpredfactor f=E15.8 b=4 &
      n='hydrostatic pressure reduction in waves 3 and 6' u=1 &
      c0='0.0: Deactivation of pressure reduction terms' &
      c1='1.0: Activation of pressure reduction terms (default)'
  0.0
  
  Possible choices are
  • 0.0: Deactivation of pressure reduction terms in waves 3 and "6" in Roe solver
  • 1.0: Activation of pressure reduction terms in waves 3 and "6" in Roe solver (default)
  This parameter is not recognized by the MHD module.

• real c_hydsdiffvelo:
  With this parameter an extra energy diffusion down an entropy step (after reconstructing the entropy) can be activated. For the diffusion velocity all three velocity components are used, i.e., there is a diffusion even transversal to the flow speed. This makes the diffusion more isotropic and lets numerical solutions to the Sedov-blast-wave problem look smoother. The value can be set e.g. with
  real c_hydsdiffvelo f=E15.8 b=4 &
      n='Parameter for energy diffusion down an entropy step' u=1 &
      c0='typically 1.0'
  1.0
  
  Possible choices are for example
  • 0.0: No extra energy diffusion down an entropy gradient (default).
  • 0.2: Small value
  • 1.0: Reasonable value
  • 3.0: Large value
  It acts on the entropy but only on the fraction of the jump that is still there after the reconstruction - unlike C_hydTdiffLin and C_hydTdiffMach that act on local temperature minima (see below). This parameter is not recognized by the MHD module.

• real c_hydvdiffvelo:
  With this parameter an extra diffusion of the transversal velocities (after reconstructing them) can be activated. For the diffusion velocity all three velocity components are used, i.e., there is a diffusion even transversal to the flow speed. This makes the diffusion more isotropic and lets numerical solutions to the Sedov-blast-wave problem look smoother. The value can be set e.g. with
  real c_hydvdiffvelo f=E15.8 b=4 &
      n='Parameter for transversal  momentum diffusion down a velocity step' u=1 &
      c0='typically 1.0'
  1.0
  
  Possible choices are for example
  • 0.0: No extra energy diffusion down an entropy gradient (default).
  • 0.2: Small value
  • 0.5: Reasonable value
  • 1.0: Reasonable value
  • 3.0: Large value
  The type of dissipation that can be added with this option is similar to that provided by the tensor viscosity (Sect. 5.4.8). Therefore, I recommend to use only small values (0.5 to 1.0) to reduce the directional dependence of the dissipation applied by the Roe solver and to use an explicite tensor viscosity (Sect. 5.4.8) in case of emergencies. This parameter is not recognized by the MHD module (but in fact, the HLLE solver applies a somewhat similar - but stronger - dissipation by construction).

• real c_hydtdifflin:
  The Roe-solver correction does not have any dissipation for standing waves and very little for slowly moving ones - unlike the MHD-HLLE solver that always has quite some minimum dissipation (depending on the reconstruction, though). The small amount of dissipation is a positive feature when it comes e.g., to convection with smaller Mach numbers in the deeper layers of a solar model, or (the absence) of vertical net mass fluxes in a hydrostatic atmosphere. However, it can be a drawback, e.g., in the case of a transsonic rarefaction wave, where an entropy fix is applied (since the beginning, see Leveque ``Numerical Methods for Conservation Laws''), or in the (rather quiet) center of a Sedov blast wave where a cross-like structure with comparably low temperature is visible. This artefact can be removed by applying the standard tensor viscosity including the energy dissipation down an entropy gradient (controlled with C_visPrturb) of by the new, rather simple option to the Roe solver to apply additional energy dissipation down a temperature gradient but only in cells adjacent to a local temperature minimum (in 1D). It works nicely for the Sedov blast wave and does not have any visible adverse impact in any other cases, so far. I don't expect it to do any harm for convection simulations and recommend to activate it per default. Cautious people might not like it and see how well it works without, though. Still, unlike the tensor viscosity with turbulent Prandtl number there is no viscous flux down the subphotospheric temperature or entropy gradient (unless there is a local minimum at the bottom). The linear diffusion can be set e.g. with
  real c_hydtdifflin f=E15.8 b=4 &
      n='Parameter 1 for energy diffusion near local temperature minimum' u=1 &
      c0='typically 0.2, if activated; usually 0.0'
  0.0
  
  It just acts on the temperature gradient (in the neighborhood of a local minimum), whereas C_hydTdiffMach (see below) has an additional factor, that takes velocity gradients into account, which prevent this type of diffusion from acting in deeper convective layers. I recommend as a default C_hydTdiffLin=0.0 and non-zero values only in case of an emergency and C_hydTdiffMach=1.0 or 2.0 as a default. These parameters are not recognized by the MHD module.

• real c_hydtdiffmach:
  Analogously to c_hydtdifflin, an additional energy dissipation down a temperature gradient but only in cells adjacent to a local temperature minimum (in 1D) can be applied, that has a scaling factor derived from the Mach number. This Mach-number-dependent diffusion can be set e.g. with
  real c_hydtdiffmach f=E15.8 b=4 &
      n='Parameter 1 for energy diffusion near local temperature minimum' u=1 &
      c0='typically 0.2, if activated; usually 0.0'
  1.0
  
  I recommend as a default C_hydTdiffLin=0.0 and non-zero values only in case of an emergency and C_hydTdiffMach=1.0 or 2.0 as a default. These parameters are not recognized by the MHD module.

• integer n_hyditer:
  After each complete hydiative time step the recommendation for the next time step will be chosen so that n_hyditer iterations will (probably) needed. The parameter can be set e.g. with
  integer n_hyditer f=I4 b=4 &
      n='Number of hydrodynamics iterations' c0=10
  8
  
  For a simulation of a solar-type star it will typically be set to 1. E.g. for brown dwarfs with shorter hydrodynamical time scales values around 10 may be considered. Note, that the hydrodynamics iteration works somewhat differently than the radiation transport iteration: in the latter case the size of the actual time step can be determined after computing the fluxes, whereas the hydrodynamics step is (possibly) of at least second order in time and the time step has to be known in advance.

  In case of the MHD module, this can reduce the computation time significantly because the MHD time step is often much smaller than the other timestep limits. However, this parameter does not determine the nuber of substeps directly. Instead, the number of substeps is determined by the requested timestep and the timestep of each MHD substep which is determined by the standart Courant condition. So, the number of MHD substeps varies during the simulation. So, the exact value of this parameter has no direct influence on the number of substeps. If this parameter is set to zero, this features is turned of. (See also parameter va_max for speeding up MHD simulations.)

• integer n_hydmaxiter:
  The absolute maximum number of hydro iterations can be specified e.g. with
  integer n_hydmaxiter f=I4 b=4 &
      n='Maximum number of hydro iterations' c0=14
  0
  
  If more iterations are needed the computation for the current time step is stopped and resumed with a smaller one. Usually, n_hydmaxiter will either be set to a value somewhat larger than the recommended number of iterations (n_hyditer) or to 0 which disables the check for too many iterations completely. This can be safely allowed in many cases. To disable the iteration of the hydrodynamics sub-step set n_hyditer=0.

  In case of the MHD module, this parameter sets the maximum number of MHD-substeps. If this parameter is set to zero, the number of MHD substeps is not limited.

• integer n_hydcellsperchunk:
  In every directional sub-step neighboring 1D columns are independent from each other. They can be grouped and computed in chunks of arbitrary size. The approximate number of grid cells per chunk can be specified e.g. with
  integer n_hydcellsperchunk f=I9 b=4 &
    n='Number of cells per hydro chunk' &
    c0='0 => one 2D slice at a time' &
    c1='1 => minimum chunk size (inefficient)' &
    c2='2500: reasonable value' &
    c3='1000000000: maximum chunk size (inefficient and memory intensive)'
  20000
  
  The exact number is determined at run time to get (approximately) equal sizes of the individual chunks. The choice of this parameter does not affect the result of the computation but the memory usage and performance: Smaller (and more) chunks may result in an optimum cache usage and need the smallest amount of memory, but result in additional overhead due to frequent subroutine calls. Bigger (and less) chunks are to be preferred for vector machines and processors with large caches. Very rough guide values may be
  • 2500: Pentium III, Core 2 Duo processor
  • 20000: RISC processor
  • 100000: Vector machine
  Note: For simulations with activated OpenMP on a parallel machine the chunk size has to be made small enough to allow at least as many chunks as processors available. This is particularly important for models with a small number of grid points (e.g. 2D models). An example is given for the Hitachi SR8000 in Sect. 3.7.8.

• real c_visbound:
  An additional drag force can be added locally in inflow cells in the outer layer when the transmitting boundary condition is chosen. The value can be set e.g. with
  real c_visbound f=E15.8 b=4 &
      n='Boundary drag viscosity parameter' u=1
  0.001
  
  This extra drag force is usually not necessary and should be switched off (with c_visbound=0.0).

The following parameters apply only to the MHD module (see also Sect. 2.2 for details).

• integer n_orderconstrainedtransport:
  Order of reconstruction of the electric field in the constrained-transport step. The parameter can be set e.g. with
  integer n_orderconstrainedtransport f=I4 b=4 &
    n='Number of hydrodynamics iterations' &
    c0='order of reconstruction of electric field in constrained-transport step' &
    c1='1: Simple arithmetic average (default), 2: Quadratic interpolation'
  2
  
  Possible values:
  • 1: Simple arithmetic average of fluxes (default).
  • 2: Quadratic interpolation of fluxes to the cell edges
  This parameter is only recognized by the hdscheme=HLLMHD.

• integer n_rescellsperchunk:
  The approximate number of grid cells per chunk in the 3D resistivity scheme can be specified e.g. with
  integer n_rescellsperchunk f=I9 b=4 &
    n='Number of cells per chunk in resistivity routine' &
    c0='0 => one 2D slice at a time' &
    c1='1 => minimum chunk size (inefficient)' &
    c2='10000: typical value'
  20000
  
  The exact number is determined at run time to get (approximately) equal sizes of the individual chunks. The choice of this parameter does not affect the result of the computation but the memory usage and performance: Smaller (and more) chunks may result in an optimum cache usage and need the smallest amount of memory, but result in additional overhead due to frequent subroutine calls. Bigger (and less) chunks are to be preferred for vector machines and processors with large caches.

• real c_rescourant:
  The tensor resistivity routines have their own time-step restriction. The recommended typical time step can be set e.g. with
  real c_rescourant f=E15.8 b=4 n='resistivity Courant factor'             u=1 &
    c0='range: 0.0 < C_resCourant, typically: 0.5'
   0.5
  

• real c_rescourantmax:
  The absolute upper limit for the resistivity time scale can be set with
  real c_rescourantmax f=E15.8 b=4 n='maximum resistivity Courant factor'  u=1 &
    c0='range: C_resCourant <= C_resCourantmax, typically 1.0'
   0.95
  
  Its value should be slightly above c_rescourant.

• real c_resB:
  This parameter specifies the electric resistivity. Values between 0.0 and 1.0 may be reasonable. Higher values are possible but drastically reduce the time step. The default value is 0.0. Example:
  real c_resb f=E15.8 b=4 &
      n='Parameter for numerical resistivity'                       u=1
  0.0
  
  The artificial electric resistivity of [Stone & Pringle (2001)] is used. It is given by $\eta_{B} = \verb\vert c_resB\vert \cdot \vert j\vert\cdot (\Delta x)^2 / \sqrt{\rho}$, where $\eta_{B}$ is the magnetic diffusivity, $j$ the current density, $\Delta x$ a mean width of the grid cell, and $\rho$ is the density. Values $\le 0$ disactivate this feature.

• real c_resBconst:
  This parameter specifies a constant magnetic diffusivity $\eta_{\rm m}$. The default value is 0.0. Example:
  real c_resbconstant f=E15.8 b=4 &
      n='Parameter for constant magnetic diffusivity'                       u=cm^2/s
  0.1
  
  A reasonable value can be found by choosing a reasonable value for the magnetic Reynolds number ${\rm Re}_{\rm mag} = (\ell_0 v_0)/\eta_{\rm m}$. Values $\le 0$ disactivate this feature.

• real c_resepsilon:
  This parameter controls an additional numerical energy diffusion. Typical values are between 0.0 and 1.0. The default value is 0.0. Example:
  real c_resepsilon f=E15.8 b=4 &
      n='Parameter for additional energy diffusion'                 u=1
  0.5
  
  The diffusion coefficient of the additional energy diffusion is given by $\eta_{E} = \verb\vert c_resepsilon\vert \cdot \vert j\vert\cdot (\Delta x)^2 / \sqrt{\rho}$, analog to the magnetic diffusivity c_resB. Values $\le 0$ disactivate this feature.

• real beta_inv:
  This parameter is used for the dual energy method. It determines which cells are updated with the thermal energy equation. Example:
  real beta_inv f=E15.8 b=4 &
      n='1/beta dual energy parameter'                               u=1
  -1.0
  
  The criterion is as follows:
  if $\beta < 1/\verb\vert beta_inv\vert$ the equation for the thermal energy is used.
  if $\beta > 1/\verb\vert beta_inv\vert$ the equation for the total energy is used.
  $\beta$ is the plasma $\beta$, i.e., the ratio of the gas pressure to the magnetic pressure. If beta_inv is set to zero, the thermal energy equation is used for all cells. A reasonable value may be 100.0. To deactivate this feature, set beta_inv to a negative value. If not specified, this parameter is set to -1.0. This parameter works only with hdscheme=HLLMHD. When $\verb\vert va_max\vert > 0$ then a reasonable value is given by $\verb\vert beta_inv\vert \lessapprox v_{A {\rm max}}2/c_s2$, where $c_s$ is the speed of sound in the region where $\beta < 1/\verb\vert beta_inv\vert$.

• real va_max:
  This parameter limits the Alfvén speed to va_max by an arbitrary reduction of the Lorentz force. Example:
  real va_max f=E15.8 b=4 &
      n='maximum Alfven speed' u=cm/s &
  -1.0
  
  This parameter is only recognized in conjunction with hdscheme=HLLMHD. va_max=5.0E+06 might be a reasonable value for some specific applications to the solar photosphere. If not specified this parameter is set -1.0. Values $\le 0$ disactivate this feature. See also Eq. 19.

• integer N_magDiffSide:
  Increases the diffusivity of the scheme near the side boundaries within a thickness of the diffusion layer of N_magDiffSide computational cells. Example:
  integer n_magdiffside f=I9 b=4 &
    n='Number of cells of diffusive side boundary layer' &
    c0='0 => deactivates this feature' &
    c1='10: reasonable value'
  0
  

• integer N_magDiffTop:
  Increases the diffusivity of the scheme near the top boundary within a thickness of the diffusion layer of N_magDiffTop computational cells. Example:
  integer n_magdifftop f=I9 b=4 &
    n='Number of cells of diffusive top boundary layer' &
    c0='0 => deactivates this feature' &
    c1='10: reasonable value'
  0
  

• integer N_magDiffBottom:
  Increases the diffusivity of the scheme near the bottom boundary within a thickness of the diffusion layer of N_magDiffBottom computational cells. Example:
  integer n_magdiffbottom f=I9 b=4 &
    n='Number of cells of diffusive bottom boundary layer' &
    c0='0 => deactivates this feature' &
    c1='10: reasonable value'
  0