# Advanced Methods for Solving for Solving Problems in Multi-Dimensions 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. Next, in Advanced Methods for Solving the Magnetohydrodynamics Equations with USim, we extended these ideas to demonstrate how to solve the MHD equations in USim, again without the use of macros.

This tutorial is based on the Rayleigh-Taylor Instability (rtInstability.pre) example in USimBase, which demonstrates the well-known Rayleigh-Taylor instability problem described by:

Jun, Norman, & Stone, ApJ 453, 332 (1995).


Using this example, we will examine how to solve the equations of inviscid compressible hydrodynamics and ideal compressible magnetohydrodynamics directly without the use of macros.

## Defining the Simulation Grid¶

In order to execute the two-dimensional Rayleigh-Taylor, we must define a two-dimensional grid simulation grid. In USim, this is done by replacing the Grid input block defined in Defining the Simulation Grid with:

<Grid domain>
kind=cart2d
lower=[-0.25  -0.75]
upper=[0.25  0.75]
cells=[64  192]
periodicDirs=[0]
ghostLayers=2
writeGeom=0
writeConn=0
writeHalos=0
</Grid>


Notice that the lower, upper and cells input blocks now take two entries, one for each dimension. A three-dimensional Cartesian grid would therefore be defined by:

<Grid domain>
kind=cart3d
lower=[-0.25  -0.75  -0.25]
upper=[0.25  0.75  0.25]
cells=[64  192  64]
periodicDirs=[0 2]
ghostLayers=2
writeGeom=0
writeConn=0
writeHalos=0
</Grid>


Note the usage of the periodicDirs attribute in these blocks. This attribute species the directions in the computational mesh that are periodic. For this example in two-dimensions, the 0 (or x) direction is periodic; similarly, in three-dimensions, both the 0 and 2 (z) directions are periodic. Note that this attribute is necessary to setup information within the simulation grid regarding periodicity; however, periodic boundary conditions must still be applied in the simulation. This is discussed below.

## Allocating Simulation Memory¶

USim data structures automatically account for the chosen dimensionality of the simulation. As a result, there is no need to change how data structures are added for the Euler equations compared to the discussion in Advanced Methods for Solving the Euler Equations with USim or the MHD equations compared to the discussion in Advanced Methods for Solving the Magnetohydrodynamics Equations with USim.

## Initializing the Fluid¶

The initial condition for a multi-dimensional simulation follows the same basic three-step pattern outlined previously; however, the kind attribute for the initialization updater has to match the dimension of the mesh:

<Updater setVar>
kind=initialize2d
onGrid=domain
in=none
out=[q   qNew]

<Function initFunc>
kind=exprFunc
gas_gamma=1.4
# Upper fluid density
rhoTop=2.0
# Lower fluid density
rhoBottom=1.0
# Perturbation amplitude
perturb=0.01
# Location of upper y boundary
ytop=0.75
# acceleration due to gravity
gravity=0.1
# X-extent of domain
lx=0.5
# Y-extent of domain
ly=1.5
preExprs=[ p0=0.01   pert=0.01   pi=3.14159   rho=if(y<0.0,rhoBottom,rhoTop)   pr=(1.0/gas_gamma)-(gravity*rho*y)   vx=0.0   vy=(0.25*perturb)*(1.0+cos(2.0*pi*x/lx))*(1.0+cos(2.0*pi*y/ly))   vz=0.0  ]
exprs=[ rho   rho*vx   rho*vy   rho*vz   (pr/(gas_gamma-1.0))+0.5*rho*vy*vy  ]
</Function>

</Updater>


As in Solving Multi-Dimensional Problems in USim, this initialization block creates a hydrostatic equilibrium consisting of a heavier fluid supported by a lighter fluid, which is perturbed with a single mode. Compared to the initialization block for the Sod shock tube problem discussed in Advanced Methods for Solving the Euler Equations with USim, the key difference as far as problem dimensionality is concerned is that the kind block is specified as initialize2d for this problem. This pattern is repeated for all updaters within a USim input file, e.g. for the generic updater update:

<Updater Updater_Name (type=string)>
kind = updateKind<NDIM>d
onGrid = <Grid_Name>
</Updater>


Here, updateKind is the kind of the USim updater that is being utilized (e.g. for the initialization input block, this is initialize) and <NDIM> is the dimensionality of the problem, e.g. 1,2,3 for one-, two- and three-dimensions.

## 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 Solving Multi-Dimensional Problems in USim, we added this algorithm through:

#
# Add the spatial discretization of the fluxes
finiteVolumeScheme(DIFFUSIVE)

#
# Add source term for gravitational acceleration to the equations


For the Rayleigh-Taylor Instability (rtInstability.pre) example in two-dimensions, this macro expands to give:

<Updater hyper>
kind=classicMuscl2d
onGrid=domain
in=[q]
out=[qNew]
cfl=0.4
timeIntegrationScheme=none
numericalFlux=hllcEulerFlux
limiter=[muscl]
variableForm=primitive
preservePositivity=0

<Equation euler>
kind=eulerEqn
gasGamma=1.4
</Equation>

equations=[euler]

<Source gravity>
kind=exprHyperSrc
gravity=0.1
rhoSrc=0.0
mxSrc=0.0
mzSrc=0.0
inpRange=[0  1  2  3  4]
outRange=[0  1  2  3  4]
indVars=[rho rhou rhov rhow Er]
exprs=[ rhoSrc   mxSrc   -rho*gravity   mzSrc   -gravity*rhov  ]
</Source>

sources=[gravity]
</Updater>


Compared to the Shock Tube (shockTube.pre) example discussed in Advanced Methods for Solving the Euler Equations with USim, we see that this is following the pattern discussed above, i.e. kind = classicMuscl1d has been replaced by kind = classicMuscl2d. We also have added an additional Source block to this updater, corresponding to the use of the macro:

# Add source term for gravitational acceleration to the equations


The use of this macro resulted in the addition of exprHyperSrc to the classicMusclUpdater (1d, 2d, 3d) updater:

<Source gravity>
kind=exprHyperSrc
gravity=0.1
rhoSrc=0.0
mxSrc=0.0
mzSrc=0.0
inpRange=[0  1  2  3  4]
outRange=[0  1  2  3  4]
indVars=[rho rhou rhov rhow Er]
exprs=[ rhoSrc   mxSrc   -rho*gravity   mzSrc   -gravity*rhov  ]
</Source>


Note that this block is following the classic three step pattern that used throughout USim (in this case, no preExpressions have been defined). The additional attributes inpRange and outRange tell USim which components of the input and output variables for the classicMusclUpdater (1d, 2d, 3d) updater to work with.

## Applying Boundary Conditions¶

The Rayleigh-Taylor instability problem requires more advanced boundary conditions that were used for the Sod Shock tube in the previous example. In two-dimensions, appropriate boundary conditions for the Rayleigh-Taylor instability are periodic boundaries in the direction transverse to the gravitational field and reflecting boundaries in the parallel direction, which are the x and y directions respectively.

In the Rayleigh-Taylor Instability (rtInstability.pre) example, we applied these boundary conditions through the use of the macros:

#
# Boundary conditions
boundaryCondition(wall,top)
boundaryCondition(wall,bottom)
boundaryCondition(periodic)


The first two of these macros expand to yield the blocks:

<Updater wallBoundaryOnEntitytop>
kind=eulerBc2d
onGrid=domain
in=[q]
out=[q]
model=eulerEqn
bcType=wall
entity=top
gasGamma=1.4
</Updater>

<Updater wallBoundaryOnEntitybottom>
kind=eulerBc2d
onGrid=domain
in=[q]
out=[q]
model=eulerEqn
bcType=wall
entity=bottom
gasGamma=1.4
</Updater>


These reflecting wall boundary conditions are described below. The third macro expands to yield a periodic boundary condition:

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


We discuss this in the next section.

### Periodic Boundary Conditions¶

Periodic boundary conditions in USim are specified in two ways, both of which must be used simultaneously. First, an input block peroidicDirs is added to the Grid block:

<Grid domain>
kind = cart2d
ghostLayers = 2
lower = [-0.25, -0.75]
upper = [  0.25,  0.75]
cells = [64, 192]
periodicDirs = [0]
</Grid>


The dimensions in which periodic boundary conditions are applied are specified by a vector. In this case, we apply periodic boundary conditions in the x direction, so we set:

periodicDirs = [0]


If we wished (for example) to extend the Rayleigh-Taylor instability problem to three-dimensions, then it might be appropriate to apply periodic boundary conditions to the x and z directions. In this case, we would set:

periodicDirs = [0,2]


In addition to the periodicDirs input block that is added to Grid, we must also add an updater block to apply the periodic boundary conditions within the simulation loop. This is done through the updater:

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


where the meanings of the various input blocks are described as in Evolving the Fluid. Note that if we wished to apply periodic boundary conditions in a three-dimensional problem, then we would set:

kind = periodicCartBc3d


This follows the pattern for changing the dimension of a USim simulation discussed previously.

### Reflecting Wall Boundary Conditions¶

In the transverse direction to the flow (here the y direction), we specify reflecting wall boundary conditions. In USim, this is done using the updater:

<Updater wallBoundaryOnEntitytop>
kind=eulerBc2d
onGrid=domain
in=[q]
out=[q]
model=eulerEqn
bcType=wall
entity=top
gasGamma=1.4
</Updater>

<Updater wallBoundaryOnEntitybottom>
kind=eulerBc2d
onGrid=domain
in=[q]
out=[q]
model=eulerEqn
bcType=wall
entity=bottom
gasGamma=1.4
</Updater>


The eulerBc (1d, 2d, 3d) applies boundary condition appropriate for systems of equations that follow the pattern of the Euler equations, i.e. a conserved state vector similar to that described at eulerEqn. The specific type of hydrodynamic equation that the boundary condition is applied to is specified by the attribute modelEqn. Valid options for this attribute include eulerEqn, realGasEqn and realGasEosEqn. Note that attributes required by these equations should be specified in the block associated with the boundary condition. In the examples above, this includes the gasGamma parameter required for eulerEqn.

A range of boundary conditions can be applied by eulerBc (1d, 2d, 3d); the specific choice is specified by the attribute bcType. Valid options for this attribute include wall, noInflow and noSlip.

Finally, we must specify the boundary that the boundary condition is to be applied on. For a structured mesh such as used in the Rayleigh-Taylor Instability (rtInstability.pre) example, a range of boundaries are predefined:

left The lower x-boundary.

right The upper x-boundary.

bottom The lower y-boundary.

top The upper y-boundary.

back The lower z-boundary.

front The upper z-boundary.

We specify the boundary to apply the boundary through the use of the entityType attribute, e.g. to apply the wall boundary condition on the upper y-boundary, we set:

entity=top


## 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 Solving Multi-Dimensional Problems in USim, we did this using the macro:

timeAdvance(TIME_ORDER)


This expands to yield:

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

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

operation=boundary
updaters=[wallBoundaryOnEntitytop wallBoundaryOnEntitybottom periodicBoundaryOnEntityghost]
syncVars=[q]

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

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

</Updater>


This time integration scheme is very similar to that described in Advanced Methods for Solving the Euler Equations with USim; we follow the pattern described above to switch the dimension of the mainIntegrator and the timeStepper blocks. To handle the changes to the boundary conditions compared to Advanced Methods for Solving the Euler Equations with USim, we change the boundaryStep block from:

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


to:

<UpdateStep boundaryStep>
operation=boundary
updaters=[wallBoundaryOnEntitytop wallBoundaryOnEntitybottom periodicBoundaryOnEntityghost]
syncVars=[q]


## 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 Rayleigh-Taylor Instability (rtInstability.pre) example, the computations performed by this macro results in the following UpdateSteps and UpdateSequence being generated:

<UpdateStep generateStep>
updaters=[]
syncVars=[]

updaters=[setVar wallBoundaryOnEntitytop wallBoundaryOnEntitybottom periodicBoundaryOnEntityghost]
syncVars=[q]

updaters=[]
syncVars=[]

updaters=[wallBoundaryOnEntitytop wallBoundaryOnEntitybottom periodicBoundaryOnEntityghost]
syncVars=[q]

updaters=[mainIntegrator copier]
syncVars=[]

updaters=[computePressure computeVelocity computeDensity]
syncVars=[]

startOnly=[generateStep   startStep]
restoreOnly=[generateStep   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. In fact, the only difference is the change in the bcStep block, which is analogous to that discussed above, i.e. we have that:

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


becomes:

<UpdateStep bcStep>
updaters=[wallBoundaryOnEntitytop wallBoundaryOnEntitybottom periodicBoundaryOnEntityghost]
syncVars=[q]


## 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. In multi-dimensions, we can extend the pattern to:

# 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>"] </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 euler> kind=eulerEqn # Adiabatic index gasGamma=<float> </Equation> # Hyperbolic equation to solve equations=[euler] <Source Souce_Name_1> kind=<string> <Source_Name_1 params> </Source> ... <Source Souce_Name_N> kind=<string> <Source_Name_N params> </Source> sources=[Souce_Name_1, ..., Souce_Name_N] </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> ... <Updater BoundaryCondition_NameN (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> # 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> # 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> <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=<String List of Updaters> syncVars=<String List of DataStructs> </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. ## Advanced Methods for Solving the MHD Equations in Multi-Dimensions¶ We can extend the approach above to solving the MHD equations in multi-dimensions in a straightforward fashion. The steps for Creating a Fluid Simulation is unchanged from Using USim to solve the Magnetohydrodynamic Equations, while the steps for Adding a Simulation Grid are identical to that given above. The first differences come in Initializing the Fluid, where the initial condition now has to specify the magnetic field, magnetic energy and the scalar potential: <Updater setVar> kind=initialize2d onGrid=domain in=none out=[q qNew] <Function initFunc> kind=exprFunc # adiabatic index gas_gamma=1.4 # Upper fluid density rhoTop=2.0 # Lower fluid density rhoBottom=1.0 # Perturbation amplitude perturb=0.01 # Location of upper y boundary ytop=0.75 # acceleration due to gravity gravity=0.1 # X-extent of domain lx=0.5 # Y-extent of domain ly=1.5 # permeability of free space mu0=1.0 # ratio of gas to magnetic pressure beta=1000.0 preExprs=[ p0=0.01 pert=0.01 pi=3.14159 rho=if(y<0.0,rhoBottom,rhoTop) pr=(1.0/gas_gamma)-(gravity*rho*y) vx=0.0 vy=(0.25*perturb)*(1.0+cos(2.0*pi*x/lx))*(1.0+cos(2.0*pi*y/ly)) vz=0.0 bx=sqrt(2.0*pr/(beta*beta)) by=0.0 bz=0.0 psi=0.0 ] exprs=[ rho rho*vx rho*vy rho*vz pr/(gas_gamma-1.0)+0.5*rho*(vy*vy)+0.5*((bx*bx)/mu0) bx by bz psi ] </Function> </Updater>  The finite volume scheme used to compute the spatial discretization is becomes: <Updater hyper> kind=classicMuscl2d onGrid=domain in=[q] out=[qNew] cfl=0.4 timeIntegrationScheme=none numericalFlux=hlldFlux limiter=[muscl] variableForm=primitive preservePositivity=0 waveSpeeds=[maxWaveSpeed] <Equation idealMhd> kind=mhdDednerEqn gasGamma=1.4 mu0=1.0 </Equation> equations=[idealMhd] <Source gravity> kind=exprHyperSrc gravity=0.1 rhoSrc=0.0 mxSrc=0.0 mzSrc=0.0 bxSrc=0.0 bySrc=0.0 bzSrc=0.0 psiSrc=0.0 inpRange=[0 1 2 3 4 5 6 7 8] outRange=[0 1 2 3 4 5 6 7 8] indVars=[rho rhou rhov rhow Er Bx By Bz Psi] exprs=[ rhoSrc mxSrc -rho*gravity mzSrc -gravity*rhov bxSrc bySrc bzSrc psiSrc ] </Source> sources=[gravity] </Updater>  Compared to the exampled discussed in Advanced Methods for Solving the Magnetohydrodynamics Equations with USim, the differences follow the same pattern as discussed above; the switch of the dimension and the addition of the Source block. Note that for the MHD equations, we have to specify the source terms for the magnetic field and the potential, even though these are zero. Note The algorithms in USim automatically preserve $$\nabla \cdot B = 0$$ (the solenoidal constraint on the magnetic field), so there is no need for the user to make changes to the algorithm at the input file level when running in multi-dimensions. Finally, the wall boundary conditions become: <Updater conductingWallBoundaryOnEntitytop> kind=mhdBc2d onGrid=domain in=[q] out=[q] model=mhdDednerEqn bcType=conductingWall entity=top gasGamma=1.4 </Updater> <Updater conductingWallBoundaryOnEntitybottom> kind=mhdBc2d onGrid=domain in=[q] out=[q] model=mhdDednerEqn bcType=conductingWall entity=bottom gasGamma=1.4 </Updater>  The signficant difference compared to the hydrodynamic case discussed above, is that we have switch the kind attribute to mhdBc (1d, 2d, 3d). Valid bcType options for this attribute include conductingWall, noInflow and noSlip, while the options for the modelType attribute include mhdDednerEqn, gasDynamicMhdDednerEqn, simpleTwoTemperatureMhdDednerEqn and twoTemperatureMhdDednerEqn. ### Advanced USim Simulation Structure for MHD in Multi-Dimensions¶ 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. A multi-dimensional approach for the MHD equations, include sources resembles 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

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]

<Source Souce_Name_1>
kind=<string>
<Source_Name_1 params>
</Source>

...

<Source Souce_Name_N>
kind=<string>
<Source_Name_N params>
</Source>

sources=[Souce_Name_1, ..., Souce_Name_N]
</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>

...

<Updater BoundaryCondition_NameN (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>

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

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>