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.
Contents
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
isRadial=1
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.
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.
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
isRadial=1
derivative=divergence
</Updater>
<Updater computeFieldCurrent>
kind=vector2d
onGrid=domain
in=[magneticField]
out=[magneticFieldCurrent]
numberOfInterpolationPoints=8
isRadial=1
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.
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
isRadial=0
writeGeom=0
writeConn=0
writeHalos=0
<Vertices vertices>
kind=funcVertCalc
<Function myGrid>
kind=exprFunc
# Step 1: Add Variables
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
# Step 2: Add Pre-Expressions
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) ]
# Step 3: Add Expressions
exprs=[ xp yp ]
</Function>
</Vertices>
</Grid>
Notice how each of the three steps are added to the function block that defines the grid condition:
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
isRadial=0
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.
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:
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\):
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>)
addEntityMaskExpression(<entityName>,<logicalExpression>)
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
<Function <entityName>Mask>
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
<Function inflowEntityMask>
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
<Function wallEntityMask>
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
<Function outflowEntityMask>
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>,)
addBoundaryConditionVariable(<boundaryCondition>,<entityName>,<variableName>,<variableValue>)
addBoundaryConditionPreExpression(<boundaryCondition>,<entityName>,<preExpression>)
addBoundaryConditionExpression(<boundaryCondition>,<entityName>,<Expression>)
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 1: Add Variables
addBoundaryConditionVariable(userSpecified,inflowEntity,gasGamma,GAS_GAMMA)
# Step 2: Add Pre-Expressions
addBoundaryConditionPreExpression(userSpecified,inflowEntity,rho = gasGamma)
addBoundaryConditionPreExpression(userSpecified,inflowEntity,vx = 3.0)
addBoundaryConditionPreExpression(userSpecified,inflowEntity,vy = 0.0)
addBoundaryConditionPreExpression(userSpecified,inflowEntity,vz = 0.0)
addBoundaryConditionPreExpression(userSpecified,inflowEntity,pr = 1.0)
# Step 3: Add expressions specifying boundary condition on density,
# momentum, total energy
addBoundaryConditionExpression(userSpecified,inflowEntity,rho)
addBoundaryConditionExpression(userSpecified,inflowEntity,rho*vx)
addBoundaryConditionExpression(userSpecified,inflowEntity,rho*vy)
addBoundaryConditionExpression(userSpecified,inflowEntity,rho*vz)
addBoundaryConditionExpression(userSpecified,inflowEntity,(pr/(gasGamma-1.0))+0.5*rho*(vx*vx))
These macros expand to yield the block:
<Updater userSpecifiedBoundaryOnEntityinflowEntity>
kind=functionBc2d
onGrid=forwardFacingStep
in=[q]
out=[q]
entity=inflowEntity
<Function userSpecifiedBoundaryOnEntityinflowEntityFunc>
kind=exprFunc
# Step 1: Add Variables
gasGamma=1.4
# Step 2: Add preExpressions
preExprs=[ rho=gasGamma vx=3.0 vy=0.0 vz=0.0 pr=1.0 ]
# Step 3: Add expressions
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:
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=[]
</UpdateStep>
<UpdateStep startStep>
updaters=[setVar wallBoundaryOnEntitywallEntity userSpecifiedBoundaryOnEntityinflowEntity copyBoundaryOnEntityoutflowEntity]
syncVars=[q]
</UpdateStep>
<UpdateStep restoreStep>
updaters=[]
syncVars=[]
</UpdateStep>
<UpdateStep bcStep>
updaters=[wallBoundaryOnEntitywallEntity userSpecifiedBoundaryOnEntityinflowEntity copyBoundaryOnEntityoutflowEntity]
syncVars=[q]
</UpdateStep>
<UpdateStep hyperStep>
updaters=[mainIntegrator copier]
syncVars=[]
</UpdateStep>
<UpdateStep writeStep>
updaters=[computePressure computeVelocity computeDensity computemachNumber]
syncVars=[]
</UpdateStep>
<UpdateSequence simulation>
startOnly=[generateStep startStep]
restoreOnly=[generateStep restoreStep bcStep]
writeOnly=[bcStep writeStep]
loop=[hyperStep]
</UpdateSequence>
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
<Function <Entity_Name1>Mask>
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
<Function <Entity_NameN>Mask>
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
# 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 generateStep>
updaters=[generate<Entity_Name1> ... generate<Entity_NameN>]
</UpdateStep>
<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>