Advanced Methods for Solving the Magnetohydrodynamics Equations with USim

In Advanced USim Simulation Concepts we examined the basic ingredients of a USim input file: the simulation grid (see Defining the Simulation Grid); data structures (see Allocating Memory); how to assign initial conditions (see Setting Initial Conditions) and how to write out additional data (see Writing Out Data). In Advanced Methods for Solving the Euler Equations with USim, we built on these concepts and demonstrated the basic methods used by USim to solve the Euler equations without the use of macros.

This tutorial is based on the Shock Tube (shockTube.pre) example. The Shock Tube simulation is designed to set up a variety of tube simulations including those by Einfeldt, Sod, Liska & Wendroff, Brio & Wu, and Ryu & Jones. In this tutorial, we will look at the Brio & Wu shock Tube, described by:

Brio, M., & Wu, C.~C. (1988), Journal of Computational Physics, 75, 400

and the use of equations for ideal compressible magnetohydrodynamics (the MHD equations) in one-dimension, which were described in Using USim to solve the Magnetohydrodynamic Equations. In this tutorial, we discuss how to setup USim algorithms for solving the MHD equations directly without using macros.

Allocating Simulation Memory

The MHD equations consist of partial differential equations that describe the evolution of density, momentum. total energy and the magnetic field. mhdDednerEqn documents the number of components that each data structure needs to contain to solve this set of equations. In addition to the quantities defined by the MHD equations, the mhdDednerEqn evolves an additional partial differential equation that controls the divergence of the magnetic field (the correction potential). As a result, our data structures need nine components:

<DataStruct q>
  kind = nodalArray
  onGrid = domain
  numComponents = 9
</DataStruct>

i.e. modifying the numComponents input block to contain 9 entries. This must be applied to all data structures that contain the conserved state, fluxes or the primitive variables in the input file.

Initializing the Fluid

The initial condition for a MHD simulation follows the same basic three-step pattern outlined previously, but with the need to specify the additional componets associated with the magnetic field and the correction potential:

<Updater init>
   kind = initialize1d
   onGrid = domain
   out = [q,qnew]

     <Function initFunc>
     kind=exprFunc

     # Step 1: Add Variables
     pi_value=3.141592653589793
     gas_gamma=1.6666666666666667
     densityL=1.0
     densityR=0.125
     pressureL=1.0
     pressureR=0.1
     normalVelocityL=0.0
     normalVelocityR=0.0
     perpendicularVelocityL=0.0
     perpendicularVelocityR=0.0
     tangentialVelocityL=0.0
     tangentialVelocityR=0.0
     mu0=1.0
     normalFieldL=0.75
     normalFieldR=0.75
     perpendicularFieldL=1.0
     perpendicularFieldR=-1.0
     tangentialFieldL=0.0
     tangentialFieldR=0.0

     # Step 2: Add Pre-Expressions
     preExprs=[ rho=if(x>0.0,densityR,densityL)   pr=if(x>0.0,pressureR,pressureL)   vx=if(x>0.0,normalVelocityR,normalVelocityL)   vy=if(x>0.0,perpendicularVelocityR,perpendicularVelocityL)   vz=if(x>0.0,tangentialVelocityR,tangentialVelocityL)   bx=if(x>0.0,normalFieldR,normalFieldL)   by=if(x>0.0,perpendicularFieldR,perpendicularFieldL)   bz=if(x>0.0,tangentialFieldR,tangentialFieldL)   psi=0.0  ]

     # Step 3: Add Expressions
     exprs=[ rho   rho*vx   rho*vy   rho*vz   pr/(gas_gamma-1.0)+0.5*rho*(vx*vx+vy*vy+vz*vz)+0.5*((bx*bx+by*by+bz*bz)/mu0)   bx   by   bz   psi  ]
   </Function>

 </Updater>

Notice how each of the three steps are added to the function block that defines the initial condition:

  1. Step 1 creates a list of <variablesName> = <value> pairs; one per line.
  2. Step 2 creates a string vector called preExprs (short for preExpressions). Each of the addPreExpression calls described in Using USim to solve the Magnetohydrodynamic Equations adds one entry to this vector, which contains the function specified in the addPreExpression call.
  3. Step 3 creates a string vector called exprs (short for expressions). As with Step 2, each of the addExpression call described in Using USim to solve the Magnetohydrodynamic Equations adds one entry to this vector, which contains the function specified in the addExpression call. Note that the order of the addExpression calls is preserved in the string vector; this is critical as these entries specify each component of the nodalArray specified in the [out] attribute.

We will see this pattern repeated whenever we used this three step process in the simple input files.

Evolving the Fluid

USim implements the well-known MUSCL scheme to compute the spatial discretization of a hyperbolic conservation law (see Using USim to solve the Euler Equations for a more detailed description). In Using USim to solve the Magnetohydrodynamic Equations, we added this algorithm through:

finiteVolumeScheme(DIFFUSIVE)

For the MHD equations, this macro expands to produce a :ref: refmanual-classicMuscl, shown below:

<Updater hyper>
  kind=classicMuscl1d
  onGrid=domain
  # input data-structures
  in=[q]
  # output data-structures
  out=[qNew]
  # CFL number to use
  cfl=0.5
  # legacy time integration scheme attribute; should be set to "none"
  timeIntegrationScheme=none
  # Riemann solver
  numericalFlux=hlldFlux
  # Limiter to use
  limiter=[muscl]
  # Form of variable to limit
  variableForm=primitive
  # Whether to check variables for physical validitiy
  preservePositivity=0
  # Maximum wave speed in the grid at this time step
  waveSpeeds=[maxWaveSpeed]

  # Hyperbolic equation system
  <Equation idealMhd>
    kind=mhdDednerEqn
    # Adiabatic index
    gasGamma=1.6666666666666667
    # Permeability of free space; should always be set to unity
    mu0=1.0
  </Equation>

  # Hyperbolic equation to solve
  equations=[idealMhd]
</Updater>

Compare to the case for the Euler equations described in Advanced Methods for Solving the Euler Equations with USim, this example of the classicMusclUpdater (1d, 2d, 3d) updater has two key differences:

  1. The use of the mhdDednerEqn to describe the hyperbolic equation system.
  2. The use of the waveSpeeds = <string vector> to specify a list of dynVectors to provide the fastest wave speed in the grid at the current time step.

The waveSpeeds = <string vector> attribute is required for the mhdDednerEqn in order to control the divergence of the magnetic field. These wave speeds are computed dynamically by USim using a timeStepRestrictionUpdater (1d, 2d, 3d), which in the shock tube example looks like:

<Updater getMaxWaveSpeed>
  kind=timeStepRestrictionUpdater1d
  onGrid=domain
  # input data-structures
  in=[q]
  # no output data-structures
  out=[]
  # CFL number to use
  courantCondition=0.5
  # DynVector to hold the maximum wave speed
  waveSpeeds=[maxWaveSpeed]

  # Time step restriction to compute
  <TimeStepRestriction idealMhdRestriction>
    # Compute a time step restriction for a hyperbolic equation
    kind=hyperbolic1d
    # Adiabatic index
    gasGamma=1.6666666666666667
    # Permeability of free-space, should always be 1.0
    mu0=1.0
    # Hyperbolic equation to use
    model=mhdDednerEqn
  </TimeStepRestriction>

  # List of time step restrictions to compute
  restrictions=[idealMhdRestriction]
</Updater>

This updater computes the global maximum wave speed over the entire simulation based on the nodalArray specified in the in string vector. The wave speed is stored in the dynVector specified in the waveSpeeds string vector.

Applying Boundary Conditions

At each time step, we have to apply boundary conditions at the left and right of the domain to ensure that at the next time step, physically-valid data is used to update the conserved state. Without this, the simulation will fail. It is possible to specify arbitrary boundary conditions in USim.

The boundary conditions used for the Brio & Wu shock tube are the same as those used in Using USim to solve the Euler Equations, Using USim to solve the Magnetohydrodynamic Equations and Advanced Methods for Solving the Euler Equations with USim: outflow (“open”) boundary conditions at both ends of the domain. This was done in Using USim to solve the Magnetohydrodynamic Equations using the macros:

boundaryCondition(copy,left)
boundaryCondition(copy,right)

These expand to yield:

<Updater copyBoundaryOnEntityleft>
  kind=copy1d
  onGrid=domain
  in=[q]
  out=[q]
  entity=left
</Updater>

<Updater copyBoundaryOnEntityright>
  kind=copy1d
  onGrid=domain
  in=[q]
  out=[q]
  entity=right
</Updater>

The copy boundary condition block is described in copy (1d, 2d, 3d). This boundary condition updater copies the values on the layer next to the ghost cells into the ghost cells - this is equivalent to a zero derivative boundary condition. Note that the blocks that describe these boundary conditions are identical to those used for the Euler equations (described in Advanced Methods for Solving the Euler Equations with USim). This is because these boundary conditions simply copy the values in each component of the nodalArray into the ghost cells; USim does this for all components of the supplied input nodalArray.

As in Using USim to solve the Euler Equations, Using USim to solve the Magnetohydrodynamic Equations and Advanced Methods for Solving the Euler Equations with USim, If we are evolving in more than one-dimension, we have to specify boundary conditions on the rest of the domain boundaries. This was done in Using USim to solve the Magnetohydrodynamic Equations using the macro:

boundaryCondition(periodic)

This expands to yield:

<Updater periodicBoundaryOnEntityghost>
  kind=periodicCartBc1d
  onGrid=domain
  in=[q]
  out=[q]
</Updater>

The periodic boundary condition updater (documented at periodicCartBc (1d, 2d, 3d)) is again identical to that used for the Euler equations (described in Advanced Methods for Solving the Euler Equations with USim), for the same reason as discussed above for the copy boundary condition. The copy and periodicCartBc boundary conditions are the two easiest boundary conditions to apply in USim.

Advancing By A Time Step

In order to solve the MHD equations, we have to advance the conserved quantities from time \(t\) to \(t+\Delta t\). This is done by applying a time integration scheme. In Using USim to solve the Magnetohydrodynamic Equations, we did this using the macro:

timeAdvance(TIME_ORDER)

This expands to yield:

<Updater mainIntegrator>
  kind=multiUpdater1d
  onGrid=domain
  in=[q]
  out=[qNew]

  <TimeIntegrator timeStepper>
    kind=rungeKutta1d
    onGrid=domain
    scheme=second
  </TimeIntegrator>


  <UpdateStep boundaryStep>
    operation=boundary
    updaters=[copyBoundaryOnEntityleft copyBoundaryOnEntityright periodicBoundaryOnEntityghost]
    syncVars=[q]
  </UpdateStep>


  <UpdateStep integrationStep>
    operation=integrate
    updaters=[hyper]
    syncVars=[qNew]
  </UpdateStep>


  <UpdateSequence sequence>
    startOnly=[]
    restoreOnly=[]
    writeOnly=[]
    loop=[boundaryStep   integrationStep]
  </UpdateSequence>

</Updater>

Note that this time integration scheme is identical to that described in Advanced Methods for Solving the Euler Equations with USim; in USim, we can evolve many different systems of hyperbolic equations using the same basic algorithmic ingredients.

Simulation Diagnostics

The macro-based simulations described in Basic USim Simulations compute a range of quantities that are output from the simulation to aid the user in understanding simulation behavior. For the MHD equations, these quantities include the primitive variables:

primitiveVariables = [density Velocity Pressure magneticField magneticFieldEnergy]

These quantities are computed in the simulation through use of a combiner (1d, 2d, 3d):

<Updater computeDensity>
  kind=combiner1d
  onGrid=domain
  in=[q]
  out=[density]
  indVars_q=[rho rhou rhov rhow Er bx by bz psi]
  exprs=[ rho  ]
</Updater>

<Updater computeVelocity>
  kind=combiner1d
  onGrid=domain
  in=[q]
  out=[velocity]
  indVars_q=[rho rhou rhov rhow Er bx by bz psi]
  preExprs=[ vx=rhou/rho   vy=rhov/rho   vz=rhow/rho  ]
  exprs=[ vx   vy   vz  ]
</Updater>

<Updater computePressure>
  kind=combiner1d
  onGrid=domain
  in=[q]
  out=[pressure]
  gas_gamma=1.6666666666666667
  mu0=1.0
  indVars_q=[rho rhou rhov rhow Er bx by bz psi]
  preExprs=[ pr=(Er-0.5*(rhou*rhou+rhov*rhov+rhow*rhow)/rho-0.5*(bx*bx+by*by+bz*bz))*(gas_gamma-1.0)  ]
  exprs=[ pr  ]
</Updater>

<Updater computeMagneticField>
  kind=combiner1d
  onGrid=domain
  in=[q]
  out=[magneticField]
  indVars_q=[rho rhou rhov rhow Er bx by bz psi]
  exprs=[ bx   by   bz  ]
</Updater>

<Updater computeFieldEnergy>
  kind=combiner1d
  onGrid=domain
  in=[magneticField]
  out=[magneticFieldEnergy]
  mu0=1.0
  indVars_magneticField=[bx by bz]
  exprs=[ (bx*bx+by*by+bz*bz)*(0.5/mu0)  ]
</Updater>

Each of these updaters follow the pattern described in Writing Out Data. For the MHD equations, USim also computes a set of variables that describe the properties of the magnetic field:

magneticFieldProperties = [magneticFieldDivergence magneticFieldCurrent]

These quantities are computed through the use of USim capabilties to compute derivatives of quantities defined on the simulation mesh (described at vector (1d, 2d, 3d)):

<Updater computeFieldDivergence>
  kind=vector1d
  onGrid=domain
  in=[magneticField]
  out=[magneticFieldDivergence]
  numberOfInterpolationPoints=8
  derivative=divergence
</Updater>

<Updater computeFieldCurrent>
  kind=vector1d
  onGrid=domain
  in=[magneticField]
  out=[magneticFieldCurrent]
  numberOfInterpolationPoints=8
  derivative=curl
</Updater>

USim capabilties for computing first- and second-order derivatives are discussed in later tutorials.

Putting it all Together

In Solving Multi-Dimensional Problems in USim, we told USim that we’re done specifying the simulation and that it can be run by calling the macro:

runFluidSimulation()

This macro collects up all of the updaters that we have added so far and, based on the sequence that we have added them, figures out how to call them to run a USim simulation. For the Brio & Wu shock case in ShockTube example, the computations performed by this macro results in the following UpdateSteps and UpdateSequence being generated:

<UpdateStep startStep>
  updaters=[setVar copyBoundaryOnEntityleft copyBoundaryOnEntityright periodicBoundaryOnEntityghost getMaxWaveSpeed]
  syncVars=[q]
</UpdateStep>


<UpdateStep restoreStep>
  updaters=[]
  syncVars=[]
</UpdateStep>


<UpdateStep bcStep>
  updaters=[copyBoundaryOnEntityleft copyBoundaryOnEntityright periodicBoundaryOnEntityghost]
  syncVars=[q]
</UpdateStep>


<UpdateStep hyperStep>
  updaters=[getMaxWaveSpeed mainIntegrator copier]
  syncVars=[]
</UpdateStep>


<UpdateStep writeStep>
  updaters=[computePressure computeMagneticField computeVelocity computeDensity computeFieldEnergy computeFieldDivergence computeFieldCurrent]
  syncVars=[]
</UpdateStep>


<UpdateSequence simulation>
  startOnly=[startStep]
  restoreOnly=[restoreStep   bcStep]
  writeOnly=[bcStep   writeStep]
  loop=[hyperStep]
</UpdateSequence>

Note the similarlity between these blocks and those discussed in Advanced Methods for Solving the Euler Equations with USim. There are two major differences:

  1. The updaters specified in the hyperStep now calls the getMaxWaveSpeed updater prior to the mainIntegrator; this computes the maximum wave speed for the simulation needed for evolving the MHD equations by \(\Delta t\).
  2. The updaters specified in the writeStep now includes all of the updaters needed to compute additional quantities for the MHD equations. These updaters are only run when USim writes out data; helping to minimize their impact on USim performance.

An Example Simulation File

The input file for the problem Shock Tube in the USimBase package demonstrates each of the concepts described above to evolve the classic Brio & Shock tube problem in one-dimensional magnetohydrodynamics. You can view the actual input blocks to Ulixes in the Setup window by choosing Save And Process Setup and then clicking on the shockTube.in file. In the .in file all macros are expanded to produce input blocks.

Advanced USim Simulation Structure

In earlier tutorials, we developed a simple pattern that could be used to design USim simulations using macros. This pattern can be repeated when we don’t use macros; however, we now need to add the grids, data structures, updaters, boundary conditions and time integration schemes by hand and then tell USim how to run them. An initial for the MHD equations looks like the following:

# Specify parameters for the specific physics problem
$ PARAM_1 = <value>
$ PARAM_2 = <value>
$ PARAM_N = <value>

# Initialize a USim simulation
# Simulation start and end times
tStart = <float>
tEnd = <float>
# Number of data files to write
numFrames = <integer>
# Initial time-step to use
initDt = <float>
# Level of feedback to give user
verbosity = <info/debug>

<Component fluids>
  kind = updaterComponent

   # Setup the grid
  <Grid Grid_Name (type=string)>
    <grid parameters>
  </Grid>

   # Create data structures needed for the simulation
  <DataStruct DataStruct_Name1 (type=string)>
    kind = nodalArray
    onGrid = <Grid_Name>
    numComponents = <int>
  </DataStruct>

  <DataStruct DataStruct_Name2 (type=string)>
    kind = nodalArray
    onGrid = <Grid_Name>
    numComponents = <int>
  </DataStruct>

   <DataStruct DataStruct_NameN (type=string)>
    kind = nodalArray
    onGrid = <Grid_Name>
    numComponents = <int>
  </DataStruct>

  # Specify initial condition
  <Updater Initialization_Updater_Name (type=string)>
    kind = initialize<NDIM>d
    onGrid = <Grid_Name>
    out = <DataStruct_Name for t = 0>

    # initial condition to use
    <Function func>
      kind = exprFunc

      # Step 1: Add Variables
      VARIABLE_1 = <float>
      VARIABLE_2 = <float>
      VARIABLE_N = <float>


      # Step 2: Add Pre-Expressions
      preExprs = [ \
          "<PreExpression_1>", \
          "<PreExpression_2>", \
          "<PreExpression_N>"
          ]

     # Step 3: Add expressions specifying initial condition on density,
     # momentum, total energy
     addExpression(<expression>)

      exprs = ["<density_expression>",  \
                    "<xMomentum_expression>",  \
                    "<yMomentum_expression>",  \
                    "<zMomentum_expression>", \
                    "<totalEnergy_expression>", \
                    "<xMagneticField_expression>",  \
                    "<yMagneticField_expression>",  \
                    "<zMagneticField_expression>", \
                    "<scalarPotential_expression>"
                    ]

    </Function>

  </Updater>

  # Add the spatial discretization of the fluxes
  <Updater FiniteVolume_Updater_Name (type=string)>
    kind=classicMuscl<NDIM>d
    onGrid=<Grid_Name>

    # input data-structures
    in=<DataStruct_Name for t^n>
     # output data-structures
    out=<DataStruct_Name for \nabla \cdot F(<DataStruct_Name for t^n>)>
    # CFL number to use
    cfl=<float>
    # legacy time integration scheme, attribute; should be set to "none"
    timeIntegrationScheme=none
    # Riemann solver
    numericalFlux=<string>
     # Limiter to use
    limiter=[<string>]
    # Form of variable to limit
    variableForm=<string>
    # Whether to check variables for physical validitiy
    preservePositivity=<int>

    # Hyperbolic equation system
    <Equation idealMhd>
      kind=mhdDednerEqn
      # Adiabatic index
      gasGamma=<float>
      # Permeability
      mu0=<float>
    </Equation>

     # Hyperbolic equation to solve
    equations=[idealMhd]
  </Updater>

  # Boundary conditions
  <Updater BoundaryCondition_Name1 (type=string)>
    kind=<string><NDIM>d
    onGrid=<Grid_Name>
    in=<DataStruct_Name for t^n>
    out=<DataStruct_Name for t^n>
    entity=<Entity_Name (type=string)>
    <Boundary Condition Parameters>
  </Updater>

  # Time integration
  <Updater TimeIntegrationUpdater_Name (type=string)>
    kind=multiUpdater<NDIM>d
    onGrid=<Grid_Name>
    in=<DataStruct_Name for t^n>
    out=<DataStruct_Name for t^n+1>

    <TimeIntegrator timeStepper>
      kind=rungeKutta<NDIM>d
      onGrid=<Grid_Name>
      scheme=<string>
    </TimeIntegrator>

    <UpdateStep boundaryStep>
      operation=boundary
      updaters=<string list of boundary conditions>
      syncVars=<DataStruct_Name at t^n>
    </UpdateStep>

    <UpdateStep integrationStep>
      operation=integrate
      updaters=<string list of integrators>
      syncVars=<DataStruct_Name for \nabla \cdot F(<DataStruct_Name for t^n>)>
    </UpdateStep>

    <UpdateSequence sequence>
      loop=[boundaryStep   integrationStep]
    </UpdateSequence>

  </Updater>

  <Updater CopierUpdater_Name (type=string)>
    kind=linearCombiner<NDIM>d
    onGrid=<Grid_Name>
    in=<DataStruct_Name for t^n+1>
    out=<DataStruct_Name for t^n>
    coeffs=[1.0]
  </Updater>

  # Output Diagnostics
  <Updater OutputDiagnosticUpdater_Name1 (type=string)>
    kind=<string><NDIM>d
    onGrid=<Grid_Name>
    in=<DataStruct_Name1, DataStruct_Name2, ..., DataStruct_NameN>
    out=<OutputDiagnosticDataStruct_Name1>

    # Step 0: Specify components of input data structures
    indVars_DataStruct_Name1 = <DataStruct_Name1_Component1, DataStruct_Name1_Component2, ..., DataStruct_Name1_ComponentN>
    indVars_DataStruct_Name2 = <DataStruct_Name2_Component1, DataStruct_Name2_Component2, ..., DataStruct_Name2_ComponentN>
    indVars_DataStruct_NameN = <DataStruct_NameN_Component1, DataStruct_NameN_Component2, ..., DataStruct_NameN_ComponentN>

    # Step 1: Add Variables
    VARIABLE_1 = <float>
    VARIABLE_2 = <float>
    VARIABLE_N = <float>

    # Step 2: Add Pre-Expressions
    preExprs = [ \
          "<PreExpression_1>", \
          "<PreExpression_2>", \
          "<PreExpression_N>"
          ]

     # Step 3: Add expressions to compute output diagnostic
     exprs = [ \
          "<expression_1>", \
          "<expression_2>", \
          "<expression_N>"
          ]
  </Updater>

  <UpdateStep startStep>
    updaters=[Initialization_Updater_Name BoundaryCondition_Name1 ... BoundaryCondition_NameN]
    syncVars=[<DataStruct_Name for t^n>]
  </UpdateStep>

  <UpdateStep restoreStep>
    updaters=<String List of Updaters>
    syncVars=<String List of DataStructs>
  </UpdateStep>

  <UpdateStep bcStep>
    updaters=[BoundaryCondition_Name1 ... BoundaryCondition_NameN]
    syncVars=[<DataStruct_Name for t^n>]
  </UpdateStep>

  <UpdateStep hyperStep>
    updaters=[FiniteVolume_Updater_Name (type=string)]
    syncVars=<String List of DataStructs>
  </UpdateStep>

  <UpdateStep writeStep>
    updaters=< OutputDiagnosticUpdater_Name1,  OutputDiagnosticUpdater_Name2,  OutputDiagnosticUpdater_NameN>
    syncVars=<OutputDiagnosticDataStruct_Name1, OutputDiagnosticDataStruct_Name2, ..., OutputDiagnosticDataStruct_NameN>
  </UpdateStep>

  <UpdateSequence simulation>
    startOnly=[generateStep   startStep]
    restoreOnly=[generateStep   restoreStep   bcStep]
    writeOnly=[bcStep   writeStep]
    loop=[hyperStep]
  </UpdateSequence>

</Component>

Note that we have not filled in some of the entries in the UpdateSteps above. How to use these UpdateSteps is discussed in later tutorials.