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.
Contents
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
isRadial=0
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
isRadial=0
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.
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.
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
# 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
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.
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
addGravitationalAcceleration(GRAVITY_ACCEL)
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
addGravitationalAcceleration(GRAVITY_ACCEL)
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.
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 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.
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
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>
<UpdateStep boundaryStep>
operation=boundary
updaters=[wallBoundaryOnEntitytop wallBoundaryOnEntitybottom periodicBoundaryOnEntityghost]
syncVars=[q]
</UpdateStep>
<UpdateStep integrationStep>
operation=integrate
updaters=[hyper]
syncVars=[qNew]
</UpdateStep>
<UpdateSequence sequence>
startOnly=[]
restoreOnly=[]
writeOnly=[]
loop=[boundaryStep integrationStep]
</UpdateSequence>
</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]
</UpdateStep>
to:
<UpdateStep boundaryStep>
operation=boundary
updaters=[wallBoundaryOnEntitytop wallBoundaryOnEntitybottom periodicBoundaryOnEntityghost]
syncVars=[q]
</UpdateStep>
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=[]
</UpdateStep>
<UpdateStep startStep>
updaters=[setVar wallBoundaryOnEntitytop wallBoundaryOnEntitybottom periodicBoundaryOnEntityghost]
syncVars=[q]
</UpdateStep>
<UpdateStep restoreStep>
updaters=[]
syncVars=[]
</UpdateStep>
<UpdateStep bcStep>
updaters=[wallBoundaryOnEntitytop wallBoundaryOnEntitybottom periodicBoundaryOnEntityghost]
syncVars=[q]
</UpdateStep>
<UpdateStep hyperStep>
updaters=[mainIntegrator copier]
syncVars=[]
</UpdateStep>
<UpdateStep writeStep>
updaters=[computePressure computeVelocity computeDensity]
syncVars=[]
</UpdateStep>
<UpdateSequence simulation>
startOnly=[generateStep startStep]
restoreOnly=[generateStep restoreStep bcStep]
writeOnly=[bcStep writeStep]
loop=[hyperStep]
</UpdateSequence>
</UpdateSequence>
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]
</UpdateStep>
becomes:
<UpdateStep bcStep>
updaters=[wallBoundaryOnEntitytop wallBoundaryOnEntitybottom periodicBoundaryOnEntityghost]
syncVars=[q]
</UpdateStep>
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.
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.
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
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]
<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>
<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.