Input parameters (full list)

Below you will find the full list of input parameters divided into modules, i.e., the origin of its definition.

Command line support

Each experiment executable offers command line support. Execute with -h argument to see available options. E.g.:

> ./photo_tr.x -h

Bifrost command-line options:

  -i, --input  input_file  // provide alternative input file (default: mhd.in)
  -v, --version            // print git version of this executable and exit
  -b, --bc.                // print standard boundary condition options and exit
  -h, --help               // print usage information and exit

mhd.in - the input file

The mhd.in input file must be set up correctly in order to run the code. Below follows comments to a specific example Note also that the mhd.in file is IDL compliant and can be run as such (this is also true for the tabparam.in file used by the square_gas EOS table).

params

 ; ************************* From   params *************************
              mx =   100
              my =   150
              mz =   100
              mb =   5
           nstep =   10
      nstepstart =   0
           debug =   0
        time_lim =  -1.000E+00
           tstop =  -1.000E+00

• mx,my,mz - number of grid points in x, y and z directions, respectively.
• mb - number of ghost zones. By default mb=5 because Bifrost uses a 6^th order differential operator (don't touch!).
• nstep - number of steps to run
• debug - 0 off, 1 on; note module specific debug switches in the below
• nstepstart - steps counted from this value, can be reset to any value
• time_lim - run-time limit in seconds (disable if time_lim < 0)
• tstop - stelar time (simulation time) limit in code units (value * u_t) (disable if tstop < 0)

parallel

 ; ************************* From parallel *************************
      periodic_x =    1
      periodic_y =    1
      periodic_z =    0
            ndim =    3
         reorder =    1
     manual_cuts =    0
            dims =    [ 1, 0, 0 ]

• periodic_[x,y,z] - sets periodicity for boundary box. Most runs will be done with 1,1,0 (Bifrost is not prepared to run periodicity in the z-direction).
• ndim - defines the dimensions of your problem (For 2.5D runs, use ndim=3 and my=1).
• reorder - always set to .true.. Used in MPI communication (don't touch)
• manual_cuts - 1 will enable manual domain decomposition which is setup from dims variable
• dims - cx cy cz, number of cuts (number of MPI processes) in x,y,z direction. If cx*cy*cz != nproc, automatic domain decomposition will be use istead

units

 ; ************************* From    units *************************
             u_l =  1.000E+08
             u_t =  1.000E+02
             u_r =  1.000E-07
             u_p =  1.000E+05
             u_u =  1.000E+06
            u_kr =  1.000E-01
            u_ee =  1.000E+12
             u_e =  1.000E+05
            u_te =  1.000E+11
            u_tg =  9.692E+03
             u_B =  1.121E+03

• u_[xx] - Multiply model values by these values to get cgs units. You can set scalling by modifying units module files in $BIFROST/UNITS/*.f90

stagger

 ; ************************* From  stagger *************************
        meshfile = "nomesh.dat"
              dx =  1.0
              dy =  1.0
              dz =  1.0

• meshfile - if ="nomesh.dat" regular mesh will be created using m[x,y,z] from params and d[x,y,z].
• dx,dy,dz - Grid cell size, disregarded if meshfile exists.

timestep

 ; ************************* From timestep *************************
             Cdt =  0.3
              dt =  1e-6 ! set in params init
         constdt =  0
     msperupdate =  0
               t =  0.0
   timestepdebug =    0

• Cdt - Courant factor. Use 0.3 or smaller for Hyman (Runge-Kutta not working).
• dt - Code will attempt to run first timestep with this dt.
• constdt - If 1, initial value of dt will not change. Simulation will be executed with constat time stepping.
• msperupdate - If 1, prints core-microseconds-per-cell-update instead of Mz/s along with immbalance i.e.,: time_it_max_in_s/(time_it_max_in_s - time_it_min_in_s)*100.0.
• t - Code starts at time t
• timestepdebug - Detailed information on which variables/terms are constraining the timestep. Set with >0.

mhd

; ************************* From      mhd *************************
             nu1 =  0.03
             nu2 =  0.05
             nu3 =  0.3
            nu_r =  0.05
         nu_r_mz =  0.0
          nu_r_k =  0
          nu_r_z =  0.0
          nu_r_w =  0.068
           nu_ee =  0.3
        nu_ee_mz =  0.0
         nu_ee_k =  0
         nu_ee_z =  0.0
         nu_ee_w =  0.068
     symmetric_e =  1
     symmetric_b =  0
            grav = -2.74
            eta3 =  0.3
          ca_max =  0.0
        mhddebug =  0
          do_mhd =  1
        mhdclean = -1
     mhdclean_ub =  0
     mhdclean_lb =  0

• nu1,nu2,nu3 - wave, advection and shock hyperdiffusion, use (of order) 0.1,0.3,0.5 or smaller.
• nu_r - The value of the mass diffusion in the upper atmosphere (bottom of grid).
• nu_r_mz - The value of the mass diffusion in the lower atmosphere (top of the grid).
• nu_r_k - The central grid-point point of the hyperbolic function that transitions between the two mass diffusion parameters nu_r and nu_r_mz. If set to a number greater than or equal to zero, this will overwrite nu_r_z.
• nu_r_z - The height, in Mm, of the mid-point of the hyperbolic function that transitions between the two mass diffusion parameters nu_r and nu_r_mz. Can be overwritten by nu_r_k.
• nu_r_w - Width of the transition between the two mass diffusion values, nu_r and nu_r_mz.
• nu_ee - The value of the energy diffusion in the upper atmosphere (bottom of grid). Use of order 0.05.
• nu_ee_mz - The value of the energy diffusion in the lower atmosphere (top of the grid).
• nu_ee_k - The central grid-point point of the hyperbolic function that transitions between the two energy diffusion parameters nu_ee and nu_ee_mz. If set to a number greater than or equal to zero, this will overwrite nu_ee_z.
• nu_ee_z - The height, in Mm, of the mid-point of the hyperbolic function that transitions between the two energy diffusion parameters nu_ee and nu_ee_mz. Can be overwritten by nu_ee_k.
• nu_ee_w - Width of the transition between the two energy diffusion values nu_r and nu_r_mz.
• symmetric_e - =1 turns on a symmetric method of calculating E. Not turned on if =0. Need to be careful with this if using different sized xy and z cells.
• symmetric_b - =1 turns on a symmetric method of calculating B. Not turned on if =0. Need to be careful with this if using different sized xy and z cells.
• grav - gravity in code units (-2.740, doesn’t change over the height of simulation)
• eta3 - scales joule heating, use 0.3 or less.
• ca_max - maximum value of Alfvén speed, will change Lorentz force to accomplish this, 0 means no maximum value.
• mhddebug - print values of many terms in MHD equations
• do_mhd - solve magnetic hydrodynamics =1 or just hydrodynamics =0.
• mhdclean - force ∇⋅B=0 every mhdclean timesteps, no cleaning if set to -1.
• mhdclean_ub - Off if =0, if =1 then Div B cleaning occurs on the boundary/ghost cells every nsteps at the lower atmospheric (upper grid) boundary. This is problematic if there is flux emergence in the simulation.
• mhdclean_lb - Off if =0, if =1 then Div B cleaning occurs on the boundary/ghost cells every nsteps at the upper atmospheric (lower grid) boundary. This is problematic if there is flux emergence in the simulation.

io

 ; ************************* From       io *************************
        snapname = "snapshot"
           isnap = 0
    large_memory = 1
           nsnap = 1000000
            nscr = 1000000
           dtaux = 0.0
             aux = "  "
          dtsnap = 100.0
          newaux = 0
         raw_aux = 0
      rereadpars = 1000000
           dtscr = 100.0
           tsnap = 0.000
            tscr = 0.000
     boundarychk = 0
     print_stats = 1
         time_io = 0
  dtsnap_precise = 0

• snapname - name root of idl, snapshot and aux files, note that an '_' will be appended to files that include a snapshot number
• isnap - current snapshot, if negative then starting from a scratch file and the absolute value is the current snapshot
• large_memory - how much of the cube that is written per write operation. Larger value == larger memorry on process ID==0 is needed.
• = -1 - Parallel IO (all MPI processes try to save snap at the same time)
• = 0 - 0 MPI Master does all the I/O in chunks of one x column
• = 1 - 1 MPI Master does all the I/O in chunks of one x-y plane
• > 1 - MPI Master does all the I/O in chunks of the whole grid
• nsnap - the number of time steps between the writing of snapshots
• nscr - the number of time steps between the writing of scratch files
• aux - which aux file variables to write to the aux file
• dtsnap - time lapse between the writing of snap shots
• dtaux - if != 0.0, aux variables have been save dtaux later then variables in snap (though newaux mechanism). If dtaux != 0.0 & newaux > 0, dtaux value will be used for time step size while progressing with 1 step of simulation to fill missing aux variables. In case of auxdt == 0.0 & newaux > 0, dt value for timestep will be use.
• newaux - if equal to snap number which already exists, performs 1 iteration to define and save auxiliary variable from aux list.
• raw_aux - =1 to turn off filering for negative values from numerical noise in qxxx terms for aux variables.
• rereadpars - number operation between re-reading params files.
• dtscr - the number of time steps between the writing of scratch files
• tsnap - time of last written snap shot, can be overwritten
• tscr - time of last written scratch file
• boundarychk - =1 if ghost zones on upper and lower boundary are to be written
• print_stats = =1 to print min/max statistics for selected MHD variables during I/O operations.
• time_io - =1 to print min/max ∆t in seconds for I/O operations
• dtsnap_precise = =1 adjusts dynamically dt to make sure that I/O will be triggered preciselly on t == tsnap+dtsnap.

math

 ; ************************* From     math *************************
           max_r = 5
        smooth_r = 3
      div_solver = 1
       div_niter = 128
         div_tol = 1.0E-05
        div_norm = 'max  '
   div_sor_niter = 1024
     div_sor_tol = 1.0E-10

• max_r - size of domain for picking maximum in diffusion quantities (don't touch)
• smooth_r - size of domain for picking smoothing in diffusion quantities (don't touch)
• div_solver - Switch to select divB cleaner: 1) for one using FFTs for Posson solver (might be slow for big parallel runs). 2) solver which uses SOR for Poisson solver (might need tunning of sor_niter and sor_tol to achive convergence)
• div_niter - max number of for divB cleaner (and/or max number of outer interation if div_solver == 2)
• div_tol - stop tolerance for divB cleaner (quantity subjects for test: divB / |B| + tiny())
• div_norm - max for abs(max), norm2 for euclidian norm (a.k.a norm2)
• div_sor_niter - max number of interation for Gauss-Seidel's method (active IFF div_solver == 2)
• div_sor_tol - global stop cryterion for Gauss-Seidel's solver (active IFF div_solver == 2)

quench

 ; ************************* From   quench *************************
            qmax =  8.000

• qmax - amplification factor for large gradients (don't touch)

eos

 ; ************************* From eos      *************************
           noneq =    0
         expieos =    0
         do_hion =    0
       do_helium =    0
   cs_from_table = 0
           gamma =  1.667
    debug_square =    0
       maxlnelim = 32.300
       minlnrlim = -34.000
    tabinputfile =    " tabparam.in  "

• noneq - if 1 use radiation table in chromosphere, if 0 then use LTE based radiation table in the chromosphere. Can only be set to 1 if the experiment provides electron densities and gas temperatures.
• expieos - if turned on, exrapolates EOS quantities for e/r above EOS table bound and/or r below table bounds. For opacity table bounds will be used. If expieos .eq. 0 EOS will stop the code if e/r or r are out of bounds.
• do_hion - turn Hion off =0 or on =1. Only set to 1 if you run an HION or HELIUM experiment
• do_helium - turn Helium off =0 or on =1. Only set to 1 if you run an HELIUM experiment
• cs_from_table - if on, speed of sound is calulated from table (if possible) though format definition: cs^s = {Dp/Dr} | s = cnst, otheriwse (default) cs=sqrt(gamma * p/r).
• gamma - ratio of specific heats used in computing speed of sound
• debug_square - =1 turns debug printouts
• maxlnelim - if e values are above this value, the ideal gas approximation will be used instead of EOS table
• minlnrlim - if r values are below this value, the ideal gas approximation will be used instead of EOS table
• tabinputfile - inputfile that contains information on EOS tables, must be there for square_gas_mpi

radiation

 ; ************************* From      rad *************************
          do_rad =  1
           dtrad =  0.000
        showinfo =  0
         quadalg = 'legacy'
      quadrature =  3
         zrefine =  1
         maxiter =  100
       taustream =  1.000E-05
        accuracy =  1.000E-03
       strictint =  1
          linear =  1
       monotonic =  1
          minbin =  1
          maxbin =  4
       dualsweep =  0
            teff =  5.854E+03

• do_rad - switch for radiative transfer
• dtrad - timestep for the radiation: If dtrad > 0.0, RT will be computed at t>=trad=tprev+dtrad only.
• showinfo - if 1, info about ray geometry and dJ will be printed.
• guadalg - Selector for different ways of seting quadrature. Default legacy, use quadrature parameter and it has only one setup for 2d and 3d case. If set to gauleg (affective only for 2D case), number of rays per octant can be set via quadrature parameter and it can be set to any positive integer. For 3d case use set_XY to select from predefine B.G. Carlson angle quadratures where X in ['a','b'] and Y in [ 2,4,6,8 ] (Y is number of rays per octant).
• quadrature - If quadalg is legacy this is a selector for the quadrature method for the integral over solid angle, =1 - 2 vertical rays, =2 - 8 rays, no vertical, Gauss-Legendre/Carlson quadrature with 1 ray/octant, =3 - 24 rays, no vertical Carlson quadrature with 3 rays/octant. For guadalg = legacy, this is a parameter to set numer of rays per PI/4.
• zrefine - The degree of z-axis refinement (1: no refinement, >1: insert up to zrefine layers)
• maxiter - Maximum number of iterations to be performed before the code gives in
• taustream - Mean optical depth threshold above which radiative transfer is computed
• accuracy - The convergence criterion for the solver
• strictint - Strictness of interpolation for the source function integral: =1 - Prevent overshooting interpolation, =0 - Only prevent negative source functions
• linear - Switches for linear interpolation onto the characteristics grid
• monotonic - Switches for monotonic interpolation onto the characteristics grid
• minbin - The minimum opacity bin to compute
• maxbin - The maximum opacity bin to compute
• dualsweep - Switch for dual sweep mode (compute S corrections during upsweeps and downsweeps). Only useful for strong scattering as it is more expensive.
• teff - efective temperature

genrad

This module computes sources of heating and cooling from input data that is read from a file. Processes that can be treated include:

• optically thin emission from the corona
• absorption of coronal radiation in the chromosphere
• various recipes for cooling and heating due to neutral H I, Ca II and Mg II
• localised non-physical heating
• heating below a threshold temperature
• if Hion is available, extra heating through absorption of Lyman alpha radiation

 ; ************************* From      rad *************************
       do_genrad = 1
      genradfile = 'radtab.dat'
    debug_genrad = 0
   incrad_detail = 0
     incrad_quad = 3
        dtincrad = 0.001
    dtincrad_lya = 0.0001
    debug_incrad = 0
       anglequad = 'legacy'

• do_genrad - switch for this module.
• genradfile - name of input file
• debug_genrad - if 1, prints info about ray geometry and convergence.
• incrad_detail - sets method for computing incrad radiation. Set to 1D flux based (=0), 1D intensity based (=1, not recommended) or 3D (=3).
• incrad_quad - selects the angle quadrature in case incrad_detail=3. It works the same as quadrature in radiation module.
• dtincrad - sets the time interval in internal units for recomputation of the radiation field. It is never bigger than dtrad from radiation.
• dtincrad_lya - sets the time interval in internal units for recomputation of the radiation field for Lyman Alpha.
• debug_incrad - turn debugging of incrad on (=1) or off (=0)
• anglequad = Selector for different ways of seting quadrature if incrad_detail=3. It works the same as guadalg in radiaton. To change number or rays for gauleg use incrad_quad.

timing

 ; ************************* From   timing *************************
          timing =    0

• timing - turn timing on (=1) or off (=0)

hion

 ; ************************* From     hion *************************
  hionhasoldpops =    1
      debug_hion =    0
        ltesnaps =    0
       firstsnap =    0
        lastsnap =    0
           vdamp =  7.000E+00
  adaptive_vdamp =    0
           icrsw =    0
         cswitch =  1.000E+00
           tcrsw =  1.000E+10
    HionAtomFile =                    "h_6.atom"
  HionNeeTabFile =           "opctab_helium.dat"
      HionSolTol = 1.000E-04
    HionSolMaxIt = 512
  HionFastSolver = F
  HionRateInterp = F

• hionhasoldpops - Start up a new run from a standard snapshot only (=0) or from a previously run HION experiment with an additional HION snapshot.
• debug_hion - Write error information to a file called debug.hion
• ltesnaps - Does nothing (=0), computes LTE hydrogen populations using e and ρ (=1) or computes LTE hydrogen populations using ne and Tg (=2) from previously computed HION snap files.
• firstsnap - starting snapshot for which to compute LTE population, use together with non-zero ltesnaps.
• lastsnap - last snapshot for which to compute LTE population, use together with non-zero ltesnaps.
• vdamp - set damping of Newton-Raphson corrections. Larger value improves stability, but increases running time.
• adaptive_vdamp - if set to =1, vdamp will set to 0.0 at first try and in case of solver failure will swtich back to default value.
• HionSolTol - desired accuracy before stopping Hion N-R solver
• HionSolMaxIt - Max number of iterations for Hion N-R solver
• HionFastSolver - if =1, switches-off pivoting in Hion LU solver (reduce computation cost but might cause problem with stability. Switch on, only for relaxed simulations)
• HionRateInterp - if =1, rate coeficients will be calculated from interpolation array created at the beginning of the run.
• icrsw -
• cswitch -
• tcrsw -
• HionAtomFile - input atom file
• HionNeeTabFile - input file for table containing electron densities and internal energy for elements other than hydrogen.

helium

 ; ************************* From     helium *************************
   HeHasOldPops = 0
 HeliumAtomFile = he_3.atom
  HeliumEUVFile = cemisstab_helium.dat
        nsolver =  4
       HeFadeIn = -1.000E+00
  HeSolverMaxIt =  512
    HeSolverTol =  1.000E-04
 HeSolverPopThr =  1.000E-12
        HeDebug =  0

• HeHasOldPops - Start up a new run from a standard snapshot (=0) or from a previously run HION experiment with an additional HION snapshot.
• HeliumAtomFile - Input atom file (file in $BIFROST/INPUT/HELIUM_TAB/).
• HeliumEUVFile - Input opacity file (file in $BIFROST/INPUT/HELIUM_TAB/).
• nsolver - Number of different solution in case of convergence error. Max number is 4 but forcing 1 might help in some cases.
• HeFadeIn - Positive value will slowly turn on photoionization and radiative recombination. This is done by multiplying the photoionization cross section the value HeFadeIn. Set it to 1e-3 and it will increase to unity over 20 solar seconds. At this point the value of HeFadeIn will automatically switch to -1, which means that no fading is taking place.
• HeSolverMaxIt - Max number of iterations for Helium N-R solver.
• HeSolverTol - Desired accuracy before stopping Helium N-R solver.
• HeSolverPopThr - Threshold for remove equation from HeSolver matrix for population wchich density number is < HeSolverPopThr * TotalHeDensityNumber.
• HeDebug - For values >0, print debug informations.

spitzer

 ; ************************* From  spitzer *************************
         spitzer = "impl"
   debug_spitzer =    0
    info_spitzer =    1
     spitzer_amp =  1.000
        theta_mg =  0.900
          dtgerr =  1.000E-06
        ntest_mg =        10
            tgb0 =  5.000E+05
            tgb1 =  1.000E+04
          tau_tg =  5.000E-03
     fix_grad_tg =    1
        niter_mg = [   2,   20,   10,   10,   30 ]
            bmin = 1.000E-05
          lgtmin = alog(8.0e3)

• spitzer - impl implicit conduction, expl explicit conduction with varying coefficient, wave new explicit implementation, none turns conduction off
• debug_spitzer - turn debugging on (don't use at present)
• info_spitzer - on gives information on convergence in finest grid interation
• spitzer_amp - mplification factor for Spitzer value κ0=1.1×10−6
• theta_mg - degree of implicitness θ, use > 0.55. For "impl", theta_mg=0.95, for "wave", theta_mg=0.60.
• dtgerr - limit to achieve in iteration before stopping iterations at finest grid
• ntest_mg - do iteration test only evey ntest_mg timesteps
• tgb0 - temperature at upper boundary if fix_grad_tg=0
• tgb1 - not used at present
• tau_tg - timescale of bringing upper boundary temperature to tgb0
• fix_grad_tg - dT/dz=0 at upper boundary
• niter_mg - number of V cycles (negative number, no muligrid), number of iterations at upper level, at mid levels going "down" the V, at mid levels going "up" the V, and at the level at the bottom of the V, i.e. the coarsest level.
• bmin - Minimum B field value for conduction to occur, if B is below this value at a point then the heat will diffuse isotropically from that location.
• lgtmin - Log10 of Tg above which qspitz [in expl or wave solver] is added to dedt.

genohm

 ; ************************* From   genohm *************************
   debug_genohm =    0

  do_genohm_amb = "onfly"
      tstep_amb = "stsf"
       eta_ambo =  1.000E+00
       eta4_amb = [ 0.100,  0.100,  0.100,  0.000 ]

     sts_method = "yy"
          sts_n =    3
      sts_max_n =   10
         sts_nu =  1.850E-03
          numin =  1.850E-03
          numax =  2.500E-01

  mts_max_n_amb =   10

      gol_tabin = "gol_coltab.in"
     gol_exptab =    1

 do_genohm_hall = "onfly"
     tstep_hall = "ntsf"
      eta_hallo =  1.000E-01
      eta4_hall = [ 0.100,  0.100,  0.100 ]

 mts_max_n_hall =   10

• debug_genohm - switch for debugging information for the Generalized Ohm's Law module.
• do_genohm_amb - '.....' chooses the method to compute the ambipolar diffusion term.
• onfly calculates the whole ambipolar diffusion term inside the code.
• table uses tables for the frequency of collisions to calculate eta_ambb=eta_amb/b^2.
• funct is used to try specific functions that you want for eta_ambb.
• const uses a constant value for eta_ambb.
• false switches the ambipolar diffusion off.
• tstep_amb - '....'
• tstep_amb(1:3) can be nts, mts or sts
  • nts No-TimeStep improvement.
  • mts MiniTimeStep method.
  • sts SuperTimeStep method.
• tstep_amb(4:4) - fix or not the method to speed up the ambipolar diffusion term.
  • f fixes the previous method for each iteration.
  • v allows the code to choice automatically the most adequate method in each iteration.
• eta_ambo - amplification factor for the ambipolar diffusion coefficient (by default 1).
• eta4_amb - [ . , . , . , . ] hyperdiffusion for the ambipolar diffusion.
• Hyperdiffusion for the fast mode (by default 0.1).
• Hyperdiffusion for the advection of the magnetic field (by default 0.1).
• Hyperdiffusion for the magnetic shocks (by default 0.1).
• Hyperdiffusion in places where the current is close to zero (see Martinez-Sykora et al. 2017, ApJ 847).
• sts_method - '..' different methods to solve the STS. 'n' means "non-variable" and 'y' means "variable"
• yy - D. Nobrega-Siverio recommends to use yy (default value).
• yn - for testing purposes.
• ny - for testing purposes.
• sts_n - Value for n for the STS method.
• sts_max_n - Max. value for n for the STS method.
• sts_nu - Value for nu for the STS method.
• numin - Min. value for nu for the STS method.
• numax - Max. value for nu for the STS method.
• mts_max_n_amb - Max. value for n for the MTS method in the ambipolar diffusion calculation.
• gol_tabin - Name of the input file for the ambipolar diffusion (by default gol_coltab.in).
• gol_exptab - Expands ambipolar diffusion table (by default .false.).
• do_genohm_hall - '.....' chooses the method to compute the Hall effect.
• onfly calculates the whole Hall term inside the code.
• funct is used to try specific functions that you want for eta_hallb.
• const uses a constant value for eta_hallb.
• false switches the Hall term off.
• tstep_hall - "...."
• tstep_hall(1:3) can be nts, mts or sts
  • nts No-TimeStep improvement.
  • mts MiniTimeStep method.
• tstep_hall(4:4) - fix or not the method to speed up the ambipolar diffusion term.
  • f fixes the previous method for each iteration.
  • v allows the code to choice automatically the most adequate method in each iteration.
• eta_hallo - amplification factor for the Hall coefficient (by default 1).
• eta4_hall - [ . , . , . ] hyperdiffusion for the Hall term.
• Hyperdiffusion for the fast mode (by default 0.1).
• Hyperdiffusion for the advection of the magnetic field (by default 0.1).
• Hyperdiffusion for the magnetic shocks (by default 0.1).
• mts_max_n_hall - Max. value for n for the MTS method in the Hall term calculation.

dampuz

; ************************* From   dampuz *************************
       damptype = "none"
       tau_damp =  1.000E-04

• damptype - Controls which type of velcity damping is enabled
• velnull - if tau_damp > 0, it will set p[x,y,z] to 0 at each iteration.
• dampuz - damps contribution in the momentum and energy equation in z direction
• dampvel - works like dampuz, but damps all momentum components
• tau_damp - Strengh of damping.

pmode_damp

 ; ************************* From           pmode_damp *************
          t_damp =  0.000
          k_damp =  0

• t_damp - Damp vertical motion if t_damp is not 0.0
• k_damp - Apply moementum damping term to pde for k gt k_damp. Ramp up damping over 20 points.

out_of_equilibrium

; ************************* From out_of_e *************************
   do_out_of_eq =    1
           steq =    0
          runse =    0
       write_se =    1
            iel =   26
         atomfl =     " fe9-15-sim-12_redlvl  "

• do_out_of_eq - This will turn on (1) or off the module (0).
• steq - If starting from initial statistical equilibrium snapshot, this should be 1. If the initial conditions already has ooe output this should be 0. So far, this has to be done manually.
• runse - This does not work. It will be removed (use write_se).
• write_se - Saves the statistical equilibrium solution (1). The name of the snapshot looks like this: snapname.seq_isnap.snap. If not interested, set to 0.
• iel - Number of the element in the periodic table, e.g., Si = 14, O = 8, Fe = 26, etc.
• atomfl - Name of the input file. Make sure that the files has lower letters and the input name also has lower letters.

ebeam

; ************************* From    ebeam *************************
       do_ebeam =    1
        verbose =    0
   no_transport =    0
 max_qbeam_diff =  1.000E-05
         ds_out =  0.010
     z_rec_ulim = -1.000E+03
     z_rec_llim =  0.000E+00
       krec_lim =  1.000E-04
   cut_en_floor =  0.000
   cut_en_roof  =  1.000E+06
    min_beam_en =  1.000E-02
 out_dep_thresh =  0.000
   min_residual =  1.000E-05
     min_dep_en =  1.000E+07
  min_stop_dist =  0.500
       max_dist =  1.000E+02
power_law_index =  4.000
qjoule_acc_frac =  0.200
         buflen = -30.0

• do_ebeam - Whether to include electron beams.
• verbose - Whether to print relevant information when simulating electron beams.
• no_transport - Whether to return before tracing any beams (only compute initial beam quantities).
• max_qbeam_diff - Maximum estimated change in qbeam before it is recomputed [same units as dedt].
• ds_out - Distance between output points along the field line [Mm].
• z_rec_ulim - Smallest depth at which to detect reconnection events [Mm].
• z_rec_llim - Largest depth at which to detect reconnection events [Mm].
• krec_lim - Sites with reconnection factors smaller than this are excluded.
• cut_en_floor - Lower cut-off energies are truncated at this lower limit [keV].
• cut_en_roof - Lower cut-off energies are truncated at this upper limit [keV].
• min_beam_en - Beams with initial power densities smaller than this are excluded [erg/s/cm^3].
• out_dep_thresh - Maximum distance outside the acceleration region energy deposition starts for the beam [Mm].
• min_residual - Beams may not be stopped until heating reaches this fraction of the initial heating.
• min_dep_en - Beams may not be stopped until deposited power per distance reaches this value [erg/s/cm].
• min_stop_dist - Beams travelling shorter than this are excluded [Mm].
• max_dist - Maximum distance a beam can propagate before propagation is terminated [Mm].
• power_law_index - If positive: fixed power law index of all electron distributions, if negative: the magnitude represents r_0/(L*R) in the acceleration model of Arnold et al. (2021), where r_0 is the characteristic island radius, L is the current sheet half-width and R is the normalized reconnection rate (reasonable values are r_0/L = 0.2 and R = 0.1)
• qjoule_acc_frac - Fraction of Joule heating allotted to electron acceleration.
• buflen - Determines the size to allocate for the MPI sending buffer. Increase the size of the buffer if you get segfaults relating to allocation in MPI_BSEND. If buflen == 0.0 no buffer is allocated. If buflen > 0.0, the buffer will be allocated to fit nint(buflen) values of type BF_REAL. If buflen < 0.0, the buffer will be allocated to fit nint(buflen*nghost) values, where nghost is the total number of ghost cells in the subdomain. (Note that buflen doesn't have to be an integer.)

Boudaries setup

standard_boundary

Open boundaries in z direction are controlled by bctypelower & bctypeupper strings. Each of them consist 9 characters which handle boundaries for different variables.

Position	Control
1	Reserved position to control type of boundaries (e.g., m for characteristic boundary conditions)
2	r - density
3	ux or px - x component of velocity of momentum (depends on experiment)
4	uy or py
5	uz or pz
6	e - energy or ee - energy/density
7	bx - x component of magnetic field
8	by
9	bz

Possible values for possitions 2-9:

 # : [!] do nothing
 % : [!] just update
 a : antisymm
 c : constant
 d : [!] div_free(Bx,By,INPUT)
 e : extrap
 f : extrap_constrained
 g : symm_gradient
 h : [!] hydrostatic
 l : extrap_log
 p : [!] potential_field(Bx,By,INPUT)
 s : symmetric
 v : vanishing
 w : vanishing_gradient

You can extract that list of possible options using command line option -b, e.g., ./photo_tr.x -b

For example bctypelower = ncggvcccd in photo_tr experiment will translate to:

 Var     std::bc routine
 ---     ---------------
 rho 'c' constant                  # fullz
 Ux  'g' symm_gradient             # fullz
 Uy  'g' symm_gradient             # fullz
 Uz  'v' vanishing                 # halfz
 ee  'c' constant                  # fullz
 Bx  'c' constant                  # fullz
 By  'c' constant                  # fullz
 Bz  'd' [!] div_free(Bx,By,INPUT) # halfz

h - hydrostatic boundary conditions

h is a special type of standard boundary condition which can be set for rho or ee [e/r] in photo_tr experiment. If h is set for density, the density in BC will be calculated though inverted EOS lookup from ee values and pressure values exponentially extrapolated from last layer pressure values with scale height.

If h is set for ee, then ee will be interpolated though inverted EOS lookup via rho and exponential pressure with scale height.

x - constrained boundary conditions

Constrained BC can be activated by changing the first character in the string bctype{lower,upper} to x. When the boundary type for a given variable is set to c (constant extrapolation), the user can specify a specific value to be applied at the boundary. This is achieved by setting the mask string to y and providing the desired value in the bc{lower,upper}_val list.

p - Potential Field Extrapolation (for Boundary Conditions)

This boundary condition (BC) type performs a potential (force-free) field extrapolation from a selected height within the computational domain. The height is controlled by the parameters pfe_{lower,upper}_s (*p*otential *f*ield *e*xtrapolate lower/upper *s*tart at). By default, this is set to -1, which means the extrapolation is performed from z index = 1 for the lower BC and from z index = ze for the upper BC.

The p type can only be set for the Bz component of the field, i.e., the 9th character in the bctype{lower/upper} string. However, like the d (div_free) boundary condition, this type of boundary will modify all three vector components (Bx, By, and Bz).

Warning: This boundary condition is computationally intensive and does not scale well with many MPI ranks due to the non-parallelized FFTs involved in the computation.

Tip

You can use this BC to perform field extrapolation throughout the entire computational domain. This can be useful for overwriting an existing field with a force-free version or to start a new simulation with a field provided only at the bottom of the domain.

This is not a standard use case, so it requires a specific approach. First, find an MPI domain decomposition without any cuts in the z-direction (see the manual_cuts option), or run with a single MPI rank if your problem fits into the RAM of one machine.

Next, set nscr to 1 to save *.scr at every iteration. Then, set the 9^th character in bctypelower to p and set pfe_lower_s to a value equal to mz, which is the number of grid points in the z-direction. This will overwrite the field in the entire box with the force-free extrapolated field from the bottom layer of the computational domain. Note that the extrapolation will be performed at each iteration, so one step will suffice.

After this step, modify *.scr.idl to remove the MPI manual cuts, set nscr to a larger value, and restore pfe_lower_s to the default -1. Then, restart your run from the *.scr.idl file to continue with the extrapolated field.

Key Variables:

• bctype{lower,upper}: A character string that controls the boundary condition types for different variables. The first character determines whether the constrained BC is enabled (x) or not.

• bc{lower,upper}_msh: A character string specifying mesh centering for boundary conditions.

• bc{lower,upper}_msk: A character string acting as a mask for the constrained BC. y enables the constrained BC for a variable, while n ignores it.

• bc{lower,upper}_val: An array of real values representing the specific values to be used for constrained BC when the setup type is c (constant).

• pfe_{lower,upper}_s: If positive integer, then it describle layer from which the potential field extrapolation is performed (only in use if bctype{lower,upper} for Bz is set to p). Default value is -1, which means the extrapolation is performed from z index = 1 (first layer at the top of the domain).

Example of Use:

To use this code, follow these steps:

1. Modify the bctype{lower,upper} string to enable the constrained BC by changing its first character to 'x'.

2. Set the boundary type for the desired variable(s) to 'c' (constant).

3. If you want to apply a specific value for the constrained BC:

4. Set the mask argument to 'y'.

5. Provide the desired value in the const_value argument.

6. Call the standard_{lower,upper}_boundary subroutine with appropriate arguments to apply the constrained BC to your data.

When all is set correctly Bifrost will inform you during the runtime about constrained values. Here is an example of such std output:

******************************* bc_lower *******************************
 $Id: bc_lower_empty (24538657) Mon Oct 2 11:51:54 2023 by M1kol4j $
    bctypelower =    xcssscsss
    bclower_msk =    -----y---
    bclower_val = [ 2.00, 0.00, 0.00, 0.00, 0.20, 0.00, 0.00, 0.00 ]
    pfe_lower_s =  -1

 -[!]- bc_lower::bctypelower setup:

    Var     std::bc routine
    ---     ---------------
    rho 'c' constant                  # fullz
    Px  's' symmetric                 # fullz
    Py  's' symmetric                 # fullz
    Pz  's' symmetric                 # halfz
    ee  'c' constant                  # fullz | set to -> 2.0000E-01
    Bx  's' symmetric                 # fullz
    By  's' symmetric                 # fullz
    Bz  's' symmetric                 # halfz

bc_lower_magnetic

Note : Since 29/10/2021 bctypelower & bctypeupper has new format but legacy (shorter) format described below will still work. To read about new format check section std_bc

 ; ************************* From    bc_lower_magnetic *************
     bctypelower = "mfgc"
        bcldebug =    0
       bcllegacy =  1
     pfe_lower_s = -1
         tau_bcl =  1.00000000E+00
      tau_ee_bcl =  1.00000000E+02
      tau_d2_bcl =  1.000E+02
      tau_d3_bcl =  1.000E+02
      tau_d4_bcl =  1.000E+02
  qjoule_amp_bcl =  1.000E+00
              s0 =  1.55399361E+01
              e0 =  7.00427964E-02
              r0 =  1.86682008E-02
             cs0 =  6.54462993E-01
              p0 =  4.81020007E-03
     nsmooth_bcl =        50
      nclean_bcl =        50
      nclean_lbl =        -4
      nclean_ubl =         5
     nextrap_bcl =         1
       naver_bcl =        -1

• bctypelower - In legacy format up to 4 characters.
• If any part of the 4 characters is absent, it will default to mdgc.
• The First character is m or not m. (Using n is standard as an alternative)
• Adjusts the equations in the ghost zones at the edges of the experiment domain using the following scheme:
  1. The ghost zones at the edges of the experiment are filled by extrapolating values from the points inside the experiment to update the ghost zone values.
  2. The time derivative equations are applied throughout the whole domain including the ghost zones.
  3. If m is switched the first character, then special characteristic boundary conditions are applied to these zones, which overwrites the standard time-derivative equations in the ghost zones. If any other letter is present then the specific ghost zone characteristic equations are not applied.
• Note this is not true of ghost zones for internal subdomains, which take data from neighbouring subdomains.
• 2^nd - 4^th letters:
  • 2^nd defines extrapolation for u_z
  • 3^rd defines extrapolation for u_x and u_y
  • 4^th defines extrapolation for Bx, By, and Bz

• 2^nd - 4^th letter options. How to extrapolate into the ghost zones:

  • f → ramps values to zero
  • e → simply extrapolates
  • c → sets a constant value using the values at the edge of the experiment
  • d → constrained extrapolation - i.e. a damped extrapolation, that doesn’t allow the value of the variable to go too wild (somewhere between values given by e and c)
  • a → antisymmetric values, antisymmetric values compared between border of experiment real part and ghost zone part (see Gaalsgaard and Nordlund)
  • g → symmetric gradient of variables inside and in ghost zones (See Nordlund Gaalsgaard)

• bcldebug - controls verbosity of debug printous, =0 to turn off any prinouts set to ={1,2,3} increase volume of information printed during execution.

• bcllegacy - =0 to run new implementation.
• tau_bcl - end incoming entropy at boundary to given vaule 0 to 1; to current value at boundary one, to s0 zero.
• tau_ee_bcl - not in use
• tau_d[2,5,6,7,8]_bcl - damping of characteristics, leave large (>>1!!)
• s0 - The parameter for boundary incoming entropy talked about for tau_bcl
• e0 - specifies a proxy of the Temperature that the code drives towards in the corona
• r0 - specifies a proxy of the Temperature that the code drives towards in the corona
• cs0 - speed of sound at upper boundary (there for reference)
• p0 - pressure at upper boundary (there for reference)
• nsmooth_bcl - specifies number of steps between hydrodynamic terms being smoothed (Magnetic terms are smoothed elsewhere)
• naver_bcl - =0 off, =1 smooths and averages hydrodynamic terms - avoid if possible!
• nclean_bcl - Cleans B on the boundary every nclean_bcl timesteps, match with nsmoothbxy_bcl
• nclean_lbl - bottom of ncleaning region (-4 = set to bottom of ghost cells)
• nclean_ubl - top of ncleaning region

bc_upper_damp

 ; ************************* From        bc_upper_damp *************
     bctypeupper = "dcs"
       bcu_udamp = "ouz"
          t_bdry =  1.00000000E+00
            rbot =  2.64239990E+02
            ebot =  2.90000000E+03
      nclean_bcu =   50
  nsmoothbxy_bcu =    0
         pb_spot =  0.00000000E+00
      spot_upper =  0.00000000E+00
        beta_lim =  1.00000000E+02

• bctypeupper - '...'

• first letter

  • d -> enter a damping regime (see routine bc_upper_late) that changes the derivatives after d/dt’s have been calculated (fixing), somewhat similarly to m in bctypelower. Damping velocities to zero and trying to put in a pressure node or not to “d”, i.e. put in a placeholder.
    • No need to have not d, even if introducing a magnetic field at the bottom.
  • 2^nd letter is for magnetic field B[x,y,z]
    • same letters as bctypelower
  • 3^rd letter is for velocity u[x,y,z]
    • same letters as bctypeupper
  • c -> constant extrapolation • s -> symmetric etc

• bcu_udamp -

• t_bdry - Timescale for driving r and e towards the parameters rbot and ebot at the convection zone boundary, and also the timescale for damping at the boundary.
• rbot - Density contributions to constant entropy at the bottom of the model, if not the granules may become unstable. If nevegative, it will calculate value from inital snapshot.
• ebot - Energy contributions to constant entropy at the bottom of the model, if not the granules may become unstable. If nevegative, it will calculate value from inital snapshot.
• nclean_bcu - cleaning every nclean_bcu time steps, currently not in use.
• nsmoothbxy_bcu -
• pb_spot -
• spot_upper -
• beta_lim -

bc_upper_emer

 ; ************************* From        bc_upper_emer *************
             bx0 =  0.00000000E+00
             by0 =  0.00000000E+00
             bz0 =  0.00000000E+00
          x0_bcu =  0.00000000E+00
          x1_bcu =  1.60000000E+01
          y0_bcu =  0.00000000E+00
          y1_bcu =  1.60000000E+01
          uz_bcu =  2.00000009E-03
           strtb =          "n"

• bx0 - x-component of magnetic field entering boundary (form depends on value of strtb).
• by0 - y-component of magnetic field entering boundary
• bz0 - z-component of magnetic field entering boundary
• x0_bcu - x-coordinate for lower bound of flux emergence, use with x1_bcu to define a rectangular strip region of the lower boundary.
• x1_bcu - x-coordinate for upper bound of flux emergence, use with x1_bcu to define a rectangular strip region of the lower boundary.
• y0_bcu - x-coordinate for lower bound of flux emergence, use with x1_bcu to define a rectangular strip region of the lower boundary.
• y1_bcu - x-coordinate for upper bound of flux emergence, use with x1_bcu to define a rectangular strip region of the lower boundary.
• uz_bcu - smoothing width for border between field entering and leaving computational domain (only active if strb=c).
• strtb - Specifies the type of field entering the domain. Possible values are:
  • n - No flux enters.
  • b - When set to 'b' (b for Boris), boundary field is driven by changing the electric field in ghost zones, as opposed to the magnetic field as for case c. The boundary is implemented in the routine efield_boundary.
  • f - Functions similarly to b, but reads the electric field from external binary files (separate file for each component).
    • elbc_fbname - Base name of the input field file. For example, if elbc_fbname = test, the field should be provided in files test_{ex,ey,ez}.dat. Input fields are expected to be in correct code units, vector orientation, and staggered position.
    • elbc_fdlen - Number of z layers in the provided field (default ==6, 5 guard zones + 1 bottom layer).
    • elbc_ftlen - Number of time slots for the external field. ==1 indicates that the field is constant in time.
    • elbc_ffadein - If positive, the external field will be faded in over the course of time equal to the set value. Time should be provided in code units, i.e., expected time in seconds divided by u_t.
  • a - Random (Alleatory) of "nseed" number of tubes.
  • c - constant flux in region
  • r - constant tube in region
  • t - horizontal tube with twist
    • rtb - Tube radius (Mm)
    • frtb - Use 2.3, which ensures 99.9% of the magnetic flux of the tube entering the domain.
    • qtb - Twist of the tube (Mm^{-1})
    • zotb - Vertical distance with respect to the boundary of the convection zone where the center of the tube is going to be located initially. Use a value greater than frtb*rtb (Mm).
    • xotb - Initial horizontal position of the center of the tube (Mm).
  • o - -loop with twist
  • z - set field to value Bz0