# 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

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

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  ]

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
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
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>

operation=boundary
updaters=[copyBoundaryOnEntityleft copyBoundaryOnEntityright periodicBoundaryOnEntityghost]
syncVars=[q]

operation=integrate
updaters=[hyper]
syncVars=[qNew]

startOnly=[]
restoreOnly=[]
writeOnly=[]
loop=[boundaryStep   integrationStep]

</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]

updaters=[]
syncVars=[]

updaters=[copyBoundaryOnEntityleft copyBoundaryOnEntityright periodicBoundaryOnEntityghost]
syncVars=[q]

updaters=[getMaxWaveSpeed mainIntegrator copier]
syncVars=[]

updaters=[computePressure computeMagneticField computeVelocity computeDensity computeFieldEnergy computeFieldDivergence computeFieldCurrent]
syncVars=[]

startOnly=[startStep]
restoreOnly=[restoreStep   bcStep]
writeOnly=[bcStep   writeStep]
loop=[hyperStep]


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.

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

VARIABLE_1 = <float>
VARIABLE_2 = <float>
VARIABLE_N = <float>

preExprs = [ \
"<PreExpression_1>", \
"<PreExpression_2>", \
"<PreExpression_N>"
]

# Step 3: Add expressions specifying initial condition on density,
# momentum, total energy

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
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>

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

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

loop=[boundaryStep   integrationStep]

</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>

VARIABLE_1 = <float>
VARIABLE_2 = <float>
VARIABLE_N = <float>

preExprs = [ \
"<PreExpression_1>", \
"<PreExpression_2>", \
"<PreExpression_N>"
]

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

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

updaters=<String List of Updaters>
syncVars=<String List of DataStructs>

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

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

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