User Guide: Setting up a generic MHD and/or Euler multiphysics simulation

Table of Contents

• Basics
• Euler equations
• Magnetohydrodynamics
  • Resistivity
  • Viscosity
• Keeping things positive
  • Combining positivity-enforcing limiters with other features
  • Flux floors
  • Non-conservative density and pressure floors
  • Artificial dissipation
  • Slope limiters
• Gradient calculations
• Cylindrical geometry
  • Extra source terms

This is a guide for users of WARPXM who want to use the physics and numerical applications that are already implemented, specificially targeted at MHD and neutral Euler fluids. It covers the basics of setting up and running a multiphysics simulation using the `warpy` toolkit. It does not cover how to implement new physics in the WARPXM framework, nor the numerical methods that WARPXM uses.

Basics

A WARPXM simulation has two top-level parameters:

• The temporal order, which indicates which Runge-Kutta integrator to use. Options include RK1, TVDRK2, SSPRK3, and others.
• The spatial order, which indicates the order of approximation of the Discontinuous Galerkin solver. This is a string such as firstOrder, secondOrder, thirdOrder, etc. The spatial order must be specified several times during the setup script, so it should be defined at the top.

When setting up a simulation, you will follow this general outline:

• Construct variables, which have one or more "components". An example would be an MHD fluid, which has eight components: density, momentum x 3, energy, and B field x 3.
• If diffusive terms are desired:
  • Construct corresponding gradient variables to hold gradients of primary variables.
  • Construct variable_adjusters for each gradient variable you need to calculate.
• Construct physics apps which act on your variables. For example, a resistive MHD simulation would have an imhd app and a resistive_diffusion app.
• Construct bc_apps to apply boundary conditions to your variables (both primary and gradient).
• Construct any variable_adjusters that you want to apply at every Runge-Kutta stage. For example, limiters adjust the solution values directly.
• Construct ic_variable_adjusters to apply initial conditions to your variables at the start of the simulation.

Once the lists of apps, bc_apps, variable_adjusters, and ic_variable_adjusters have been constructed, you will bring them together in a simulation:

• Construct a spatial_solver which will solve the apps in the DG formulation.
• Construct a boundary_conditions variable adjuster whose job is to execute the bc_apps, and add it to the variable_adjusters list.
• Construct a host_action to apply the ic_variable_adjusters
• Construct a time integrator from the spatial_solver and the variable_adjusters list
• Finally, construct the simulation from the time integrator and the initial conditions host_action.

These are the steps for a moderately complicated simulation. The absolute minimum to solve, say, advection with periodic boundary conditions are a variable, an advection app, an ic_variable_adjuster to set up the initial condition, a spatial_solver, and a time integrator.

Euler equations

Refer to the shock_1d_sod.py example.

Magnetohydrodynamics

To follow along with this section, refer to the shock_1d_brio-wu.py example script.

The convention in WARPXM is to use eight components for MHD fluids, no matter the dimensionality of the overall problem. This is typically baked into the C++ implementation of variables, which look for e.g. energy in the fifth index of the component vector. As such, the way to create an MHD variable is as follows:

plasma = warpy.variable(
    name='plasma',
    components=['rho', 'px', 'py', 'pz', 'e', 'Bx', 'By', 'Bz'],
    basis_array_set=sorder_str
    )

The components should appear in that order: density, momentum x 3, energy, and B field x 3.

The actual ideal MHD fluxes are implemented by the WmApplication_IMHD_Flux application, which is constructed via the imhd constructor in warpy:

apps.append(
    warpy.apps.mhd.imhd.imhd(
        name='imhd',
        gamma=gamma,
        variable=plasma,
        numerical_flux_type=num_flux
        )
    )

This will include the ideal MHD governing equations in the Runge-Kutta update. For more details on the meaning of the floor_density and floor_pressure parameters, see the section on flux floors.

Resistivity

Resistive MHD is implemented by the WmApplication_Resistive_Diffusion application. This is constructed in warpy via the resistive_diffusion constructor:

if use_resistivity:
    apps.append(
        warpy.apps.mhd.rmhd.resistive_diffusion(
            name='resistive_diffusion',
            ...

For information about the arguments, see the documentation of WmApplication_Resistive_Diffusion. Note that, since resistive diffusion is a second-order differential operator term, we need to set up calculation of gradient variables. See Gradient calculations for information on how to do that.

Viscosity

TODO

Keeping things positive

WARPXM implements a positivity-preserving Discontinuous Galerkin scheme, which is guaranteed to preserve positivity of solution cell averages for certain combinations of fluxes and source terms. A detailed writeup with references can be found in the Positivity-preserving DG writeup (PDF). When used correctly, this scheme should provide a robust solution that never crashes due to a negative density or pressure. This section describes the correct use of the positivity-preserving elements of WARPXM and the recommended combinations. It then describes some of the alternative limiters that are implemented.

From a user perspective, the main thing that is required to take advantage of positivity-preserving DG is to install one of the "positivity-enforcing" limiters. The role of the positivity-enforcing limiter is to categorically guarantee a positive cell average at the next RK stage, provided that a CFL condition is satisfied. To install this limiter:

• For an Euler fluid, the NDGEulerPositivityPreservingLimiter:

  var_adj.append(warpy.variable_adjusters.limiters.dg_euler_positivity_enforcing_limiter(name='pos_limiter_atoms',
                                                                                         priority=2,
                                                                                         spatial_order=sorder_str,
                                                                                         fluid_variables=[atoms],
                                                                                         rho_rhou_e_B_components=[atoms_component_names],
                                                                                         include_gaussian_quad_nodes=cyl_geometry,
                                                                                         include_interior_lgl_nodes=True))

• For an MHD fluid, the NDGMHDPositivityPreservingLimiter

  var_adj.append(warpy.variable_adjusters.limiters.dg_mhd_positivity_enforcing_limiter(name='pos_limiter_plasma',
                                                                                       priority=2,
                                                                                       spatial_order=sorder_str,
                                                                                       fluid_variables=[plasma],
                                                                                       rho_rhou_e_B_components=[plasma_component_names],
                                                                                       include_gaussian_quad_nodes=cyl_geometry),
                                                                                       include_interior_lgl_nodes=True)

  Note the priority=2. It is critical that these limiters run last out of all the variable adjusters that are acting on their respective variables. For example, if you have a moment-slope limiter installed like the Moe-Rossmanith limiter, its priority parameter must be lower (indicating that it will run before) the priority of the positivity-enforcing limiter.

If your simulation is running in cylindrical geometry and you are using the GQ (Gaussian quadrature) source quadrature style (see Cylindrical geometry, then you must pass include_gaussian_quad_nodes=True. Similarly, if you are using LGL source quadrature, and your simulation includes additional source terms, then you should pass the include_interior_lgl_nodes=True parameter. There is no issue with passing True for these parameters when they are not needed, merely a performance hit from evaluating the solution at more nodes than necessary.

Combining positivity-enforcing limiters with other features

The positivity-preserving DG scheme plays more or less nicely with other WARPXM physics features.

Recommended:

• Strong stability preserving RK method: It is highly recommended to use an SSP RK method. Examples include RK1, TVDRK2, and SSPRK3.
• Flux floors: It is recommended to use the flux floors in combination with positivity-enforcing limiters. Without the flux floors, it is possible to have floating point exceptions during flux calculations due to pointwise negative density or pressure.
• Slope limiters: In the absence of a slope limiter to smooth out discontinuities in the solution, the positivity-preserving DG method may result in the recommended dt "tanking" all the way to machine precision. A combination of a slope limiter with a positivity-enforcing limiter is very robust.
  • Note: recall the admonition above about variable adjuster priority. Any slope limiter must have its priority parameter set lower than the priority of the positivity-enforcing limiter.

Supported:

• Arbitrary source terms: as long as the source term has its recommended dt set by the positivity_preserving_source_dt estimator, source terms present no difficulty to the positivity preservation scheme. Note: this applies to source terms arising from cylindrical geometry as well.
• Viscosity and resistivity: these diffusive terms add additional fluxes which technically violate the guarantees of the positivity-preserving DG scheme. However, for moderate levels of diffusivity, the limiter appears to work out of the box.

Discouraged:

• The non-conservative floor apps are discouraged, as they inject additional mass and energy into the simulation, which may result in instability.
• Artificial dissipation is discouraged, since it is an additional flux that is not covered by the positivity-preserving DG analysis, and which may be quite large in particularly sensitive regions of the simulation.

Flux floors

Despite our best efforts, it is possible for a high-order DG solution to take on negative values of density and pressure at the collocation nodes. This is a nuisance, but in itself does not necessarily indicate that the solution is destroyed. WARPXM includes options that allow the simulation to continue even if unphysical values are predicted at a collocation node.

The Euler and Ideal MHD applications provide the option to set floors on the value of density and pressure that will be used when calculating the flux functions. These options do not affect conservation of mass or energy. Rather, they will modify the values of density and pressure in the internal and numerical flux calculations. The update of the conserved variables occurs in the same way as always, just with a slightly different amount of flux in or out of the element.

To use these options, pass the floor_density and floor_pressure options to the Euler or Ideal MHD application. For example:

apps.append(
    warpy.apps.mhd.imhd.imhd(
        ...
        floor_density=rho_min,
        floor_pressure=P_min,
        ...
        )
    )

Non-conservative density and pressure floors

The following applications apply a non-conservative floor to the density and internal energy of the solution:

• SetFloorsSingleSpecies for an Euler-like fluid
• SetFloors for MHD

Artificial dissipation

TODO

Slope limiters

Slope limiters work in concert with the positivity-enforcing limiter to avoid crashes. While the positivity-enforcing limiters can keep the density and pressure positive, they cannot guarantee that the solution wavespeeds remain reasonable. For this, a limiter specifically designed to smooth out discontinuities is very helpful.

Currently the only supported slope limiter is the Moe-Rossmanith limiter.

Gradient calculations

TODO

Cylindrical geometry

• TODO: document radial index idea

WARPXM supports simulations in cylindrical geometry in both 1D and 2D. Cylindrical geometry support is built in to applications, and is not a property of the mesh as such. To use cylindrical geometry, follow these steps:

• Pass the keyword argument

  coordinate_system="Cylindrical"

  to those applications which support it. Examples include
  • WmApplication_Resistive_Diffusion
  • intraspecies_collisions (MHD)
• Switch to using Gaussian quadrature nodes for source term evaluation.
• Add additional apps for the source terms arising from cylindrical geometry (see this section

Extra source terms

When using cylindrical geometry, the divergence and curl operators introduce additional source terms due to the coordinate system. These must be included explicitly by the user in the list of applications. Applications which require a cylindrical geometry source term, and the corresponding source term app, are as follows:

• Ideal MHD: IMHD CylSource
• Resistive MHD: RMHD CylSource
• Viscous MHD: IntraspeciesCyl