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. Then, in Advanced Methods for Solving for Solving Problems in Multi-Dimensions with Usim, we demonstrated how to extend these concepts to multi-dimensions for both the Euler and MHD without the use of macros.

In this tutorial, we show how USim can be used to solve problems on more advanced meshes, building on the concepts presented in Solving Problems on Advanced Structured Meshes in USim and Solving Problems on Unstructured Meshes in USim to demonstrate how to solve problems in USim on a range of advanced mesh types.

## Solving Problems in Axi-Symmetric Curvilinear Coordinates¶

An example of using USim to solve the MHD equations in axi-symmetric curvilinear coordinates is found in Unstable Plasma zPinch (zPinch.pre).

In Solving Problems on Advanced Structured Meshes in USim, the use of axisymmetric cylindrical coordinates in this problem is specified through the use of the macro:

addCylindricalGrid(lowerBounds, upperBounds, numCells, periodicDirections)


For the Unstable Plasma zPinch (zPinch.pre) example, this macro expans to yield:

<Grid domain>
kind=cart2d
lower=[0.001  -0.1]
upper=[0.1  0.1]
cells=[32  64]
periodicDirs=[1]
ghostLayers=2
writeGeom=0
writeConn=0
writeHalos=0
</Grid>


The isRadial attribute in this Grid block specifies that the grid is utilizing curvilinear coordinates. When we work directly with the input file, we have to manually specify which Grids and updaters we want to use curvilinear coordinates, unlike for the macros described in Solving Problems on Advanced Structured Meshes in USim.

Note

Cylindrical grids in USim are designed with axisymmetric coordinates in mind, so a one-dimensional mesh simulates the $$R$$ coordinate, a two-dimensional mesh simulates the $$(R,Z)$$ coordinates and a three-dimensional mesh simulates the $$(R,Z,\phi)$$ coordinates.

### Evolving the Fluid¶

The MUSCL scheme implemented in USim requires additional source blocks to correctly integrate a hyperbolic conservation law in curvilinear coordinates. In Solving Problems on Advanced Structured Meshes in USim, these terms were added automatically when we used the addCylindricalGrid macro. When we work directly with the input file, we have to add these by hand. This is done through the addition of a Source block to the classicMusclUpdater (1d, 2d, 3d) Updater in a similar fashion to the gravitational acceleration block discussed in Advanced Methods for Solving for Solving Problems in Multi-Dimensions with Usim:

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

<Equation idealMhd>
kind=mhdDednerEqn
basementPressure=20047.379936313715
basementDensity=8.000000000000001e-06
gasGamma=1.6667
mu0=1.0
correctNans=1
correct=1
</Equation>

equations=[idealMhd]

<Source axiSymmetricSource>
kind=mhdSym
gasGamma=1.6667
mu0=1.0
symmetryType=cylindrical
model=mhdDednerEqn
</Source>

sources=[axiSymmetricSource]
</Updater>


The axiSymmetricSource applied here is documented at mhdSym. This is the only changed needed for the classicMusclUpdater (1d, 2d, 3d) updater to work in axisymmetric cylindrical coordinates.

### 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, 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=vector2d
onGrid=domain
in=[magneticField]
out=[magneticFieldDivergence]
numberOfInterpolationPoints=8
derivative=divergence
</Updater>

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


Compared to the Advanced Methods for Solving the Magnetohydrodynamics Equations with USim, these blocks have an additional attribute, isRadial=1, which tells USim to calculate these derivatives in a form consistent with axisymmetric cylindrical coordinates.

## Advanced Methods for Solving Problems on Body-Fitted Meshes in USim¶

An example of applying USim to a problem on a body-fitted mesh is found in Magnetized Ramp Flow (rampFlow.pre). This simulation uses a bodyFitted (1d, 2d, 3d) grid, which is added to the simulation through the addBodyFittedGrid macro:

addBodyFittedGrid(lowerBounds, upperBounds, numCells, periodicDirections)


This macro expands to yield:

<Grid domain>
kind=bodyFitted2d
lower=[0.0  0.0  0.0]
upper=[1.0  1.0  1.0]
cells=[37  25  25]
periodicDirs=[]
ghostLayers=2
writeGeom=0
writeConn=0
writeHalos=0

<Vertices vertices>
kind=funcVertCalc

<Function myGrid>
kind=exprFunc

xmin=0.1
xmax=0.8
ymax=0.45
slope=0.17632698070846498
inletCells=12.0
rampCells=25.0
yCells=25.0
dxc=0.02702702702702703
dyc=0.04

preExprs=[ ix=rint(x/dxc)   iy=rint(y/dyc)   dxi=xmin/inletCells   dxr=(xmax-xmin)/rampCells   xp=if(inletCells>=ix,dxi*ix,xmin+dxr*(ix-inletCells))   dyi=ymax/yCells   b=-slope*xmin   y0=slope*xp+b   dyr=(ymax-y0)/yCells   yp=if(inletCells>=ix,dyi*iy,y0+dyr*iy)  ]

exprs=[ xp   yp  ]
</Function>

</Vertices>

</Grid>


Notice how each of the three steps are added to the function block that defines the grid 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 addGridPreExpression calls described in Solving Problems on Advanced Structured Meshes in USim adds one entry to this vector, which contains the function specified in the addGridPreExpression call.
3. Step 3 creates a string vector called exprs (short for expressions). As with Step 2, each of the addGridExpression call described in Solving Problems on Advanced Structured Meshes in USim adds one entry to this vector, which contains the function specified in the addExpression call. Note that the order of the addGridExpression calls is preserved in the string vector; this is critical as these entries specify each of the real space coordinates for the computational grid.

## Advanced Methods for Solving Problems on Unstructured Meshes with USim¶

An example of applying USim to a problem on an unstructured mesh is found in Flow over a Forward-Facing Step (forwardFacingStep.pre). In Solving Problems on Unstructured Meshes in USim, we added an unstructured mesh to a USim simulation through the use of one of the two macros:

addExodusGrid(GRIDFILE)


or:

addGmshGrid(GRIDFILE)


The choice of which block to use corresponds to the format of the mesh; either GMSH or ExodusII format. The GRIDFILE is the name of the file containing the mesh without the file extension.

These macros expand to yield a grid block that follows the pattern:

<Grid GRIDFILE>
kind=unstructured
ghostLayers=2
writeGeom=0
writeConn=0
writeHalos=0
decomposeHalos=0

<Creator ctor>
kind=<Mesh_Format>
ndim=<NDIM>
file=GRIDFILE.<Mesh_Extension>
</Creator>

</Grid>


Here, <Mesh_Format> = gmsh,exodus and <Mesh_Extension> = msh,g correspondigly, while <NDIM> specifies the dimension of the simulation (as usual); here, it additionally corresponds to the highest dimensionality element in the mesh, i.e. for quads <NDIM> = 2, while <NDIM> = 3 for hexes.

Note

Unstructured meshes in USim can only have dimension 2,3. One-dimensional meshes should make use of a cartesian or a body fitted mesh.

### Advanced Methods for Boundary Conditions on Unstructured Meshes¶

As was discussed in Solving Problems on Unstructured Meshes in USim, the complexity of performing calculations on an unstructured mesh in USim is associated with application of boundary conditions. USim’s approach to applying boundary conditions on an unstructured mesh is a two-step process:

1. The user defines the regions of the mesh (entity) that will be used to apply boundary conditions.
2. The user specifies the boundary conditions to apply on each region of the mesh.

For the specific case of Flow over a Forward-Facing Step (forwardFacingStep.pre), there are three boundaries that we need to define, which are delimited by position in the streamwise direction, $$x$$:

1. Inflow boundary: defined for $$x < 0.0$$
2. Wall boundary: defined for $$0.0 < x < 3.0$$
3. Outflow boundary: defined for $$x > 3.0$$

In Solving Problems on Unstructured Meshes in USim, we generated boundary entities that exist on the exterior of the mesh using a combination of the createNewEntityFromMask and addEntityMaskExpression using the macros:

createNewEntityFromMask(<entityName>)


These expand to yield a block that follows the pattern:

<Updater generate<entityName> >
kind=entityGenerator<NDIM>d
onGrid=<gridName>
in=none
out=none
newEntityName=<entityName>
onEntity=ghost

kind=exprFunc
exprs=[ <logicalExpression> ]
</Function>

</Updater>


Of particular note here is the attribute onEntity = ghost. This attribute specifies an existing entity in the grid which is will be subdivided by <logicalExpression>. All USim grids include an entity named ghost which corresponds to the exterior surface of the mesh; as a result, our new entities only include mesh elements on the exterior surface of the mesh that obey <logicalExpression>. In this fashion, it is possible to construct complex boundary regions for unstructured meshes in USim using only simple logical conditions.

For the three entities needed for Flow over a Forward-Facing Step (forwardFacingStep.pre), these operations take the form:

# Inflow for x < 0.0
<Updater generateinflowEntity>
kind=entityGenerator2d
onGrid=forwardFacingStep
in=none
out=none
newEntityName=inflowEntity
onEntity=ghost

kind=exprFunc
exprs=[ if(x<0.0,1.0,-1.0)  ]
</Function>

</Updater>

# Wall for 0.0 < x < 3.0
<Updater generatewallEntity>
kind=entityGenerator2d
onGrid=forwardFacingStep
in=none
out=none
newEntityName=wallEntity
onEntity=ghost

kind=exprFunc
exprs=[ if((x>0.0)and(x<3.0),1.0,-1.0)  ]
</Function>

</Updater>

# Outflow for x > 3.0
<Updater generateoutflowEntity>
kind=entityGenerator2d
onGrid=forwardFacingStep
in=none
out=none
newEntityName=outflowEntity
onEntity=ghost

kind=exprFunc
exprs=[ if(x>3.0,1.0,-1.0)  ]
</Function>

</Updater>


Now that we have created the inflowEntity, wallEntity and outflowEntity, we can specify boundary conditions on them. For the wallEntity and the outflowEntity, the boundary conditions are familiar from our previous tutorials:

boundaryCondition(wall,wallEntity)
boundaryCondition(copy,outflowEntity)


The first of these boundary condition macros expands to yield:

<Updater wallBoundaryOnEntitywallEntity>
kind=eulerBc2d
onGrid=forwardFacingStep
in=[q]
out=[q]
model=eulerEqn
bcType=wall
entity=wallEntity
gasGamma=1.4
</Updater>


while the second of these expands to yield:

<Updater copyBoundaryOnEntityoutflowEntity>
kind=copy2d
onGrid=forwardFacingStep
in=[q]
out=[q]
entity=outflowEntity
</Updater>


These boundary conditions should be familiar from previous tutorials; the difference here being that the entityName parameter is set to match the specific entity that the boundary condition is applied to.

For the inflowEntity, we specify a new type of boundary condition userSpecified to determine the inflow properties using the macros:

boundaryCondition(<boundaryCondition>,<entityName>,)


For Flow over a Forward-Facing Step (forwardFacingStep.pre), the inflow boundary was specified according to:

# Add a user specified boundary condition
boundaryCondition(userSpecified,inflowEntity)

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


These macros expand to yield the block:

<Updater userSpecifiedBoundaryOnEntityinflowEntity>
kind=functionBc2d
onGrid=forwardFacingStep
in=[q]
out=[q]
entity=inflowEntity

<Function userSpecifiedBoundaryOnEntityinflowEntityFunc>
kind=exprFunc
gasGamma=1.4
preExprs=[ rho=gasGamma   vx=3.0   vy=0.0   vz=0.0   pr=1.0  ]
exprs=[ rho   rho*vx   rho*vy   rho*vz   (pr/(gasGamma-1.0))+0.5*rho*(vx*vx)  ]
</Function>

</Updater>


This block follows the standard three step pattern:

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 addBoundaryConditionPreExpression calls described in Solving Problems on Unstructured Meshes in USim adds one entry to this vector, which contains the function specified in the addBoundaryConditionPreExpression call.
3. Step 3 creates a string vector called exprs (short for expressions). As with Step 2, each of the addBoundaryConditionExpression call described in Solving Problems on Unstructured Meshes in USim adds one entry to this vector, which contains the function specified in the addBoundaryConditionExpression call. Note that the order of the addBoundaryConditionExpression calls is preserved in the string vector; this is critical as these entries specify each of the entries in the output data structure.

### Putting it all Together¶

In Solving Problems on Unstructured Meshes 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 Flow over a Forward-Facing Step (forwardFacingStep.pre) example, the computations performed by this macro results in the following UpdateSteps and UpdateSequence being generated:

<UpdateStep generateStep>
updaters=[generatewallEntity generateinflowEntity generateoutflowEntity]
syncVars=[]

updaters=[setVar wallBoundaryOnEntitywallEntity userSpecifiedBoundaryOnEntityinflowEntity copyBoundaryOnEntityoutflowEntity]
syncVars=[q]

updaters=[]
syncVars=[]

updaters=[wallBoundaryOnEntitywallEntity userSpecifiedBoundaryOnEntityinflowEntity copyBoundaryOnEntityoutflowEntity]
syncVars=[q]

updaters=[mainIntegrator copier]
syncVars=[]

updaters=[computePressure computeVelocity computeDensity computemachNumber]
syncVars=[]

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


Notice that the Updaters that generate the wall, inflow and outflow entities are added to the the generateStep. This step is only called during the startOnly and restoreOnly operations in the UpdateSequence: an entity only needs to be created once per simulation startup.

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. For unstructured meshes, 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)>
kind=unstructured
<Grid Parameters>

<Creator ctor>
kind=<Mesh_Format>
ndim=<NDIM>
file=GRIDFILE.<Mesh_Extension>
</Creator>
</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>

<Updater generate<Entity_Name1> >
kind=entityGenerator<NDIM>d
onGrid=<Grid_Name>
in=none
out=none
newEntityName=<Entity_Name1>
onEntity=ghost

kind=exprFunc
exprs=[ <logicalExpression> ]
</Function>
</Updater>

...

<Updater generate<Entity_NameN> >
kind=entityGenerator<NDIM>d
onGrid=<Grid_Name>
in=none
out=none
newEntityName=<Entity_NameN>
onEntity=ghost

kind=exprFunc
exprs=[ <logicalExpression> ]
</Function>
</Updater>

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

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

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>

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

updaters=[generate<Entity_Name1> ... generate<Entity_NameN>]

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>