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 this tutorial, we build on these concepts and demonstrate the basic methods used by USim to solve the Euler equations without the use of macros.
This tutorial is based on the Shock Tube (shockTube.pre) example. The Shock Tube simulation is designed to set up a variety of tube simulations including those by Einfeldt, Sod, Liska & Wendroff, Brio & Wu, and Ryu & Jones. In this tutorial, we will look at the Sod Shock Tube based on the classic paper:
Sod, Gary A. "A survey of several finite difference methods for
systems of nonlinear hyperbolic conservation laws." ; Journal of Computational Physics 27.1 (1978): 1-31.
and the use of equations for inviscid compressible hydrodynamics (the Euler equations), which were described in Using USim to solve the Euler Equations. The addition of a basic grid and how to setup an initial condition was covered in Advanced USim Simulation Concepts. In this tutorial, we discuss how to setup USim algorithms for solving the Euler equations directly without using macros.
Contents
The Euler equations consist of partial differential equations that describe the evolution of density, momentum and total energy. eulerEqn documents the number of components that each data structure needs to contain to solve this set of equations; which is five:
<DataStruct q>
kind = nodalArray
onGrid = domain
numComponents = 5
</DataStruct>
USim implements the well-known MUSCL scheme to compute the spatial discretization of a hyperbolic conservation law (see Using USim to solve the Euler Equations for a more detailed description). In Using USim to solve the Euler Equations, we added this algorithm through:
finiteVolumeScheme(DIFFUSIVE)
For the Euler equations, this macro expands to produce a :ref: refmanual-classicMuscl, shown below:
<Updater hyper>
kind=classicMuscl1d
onGrid=domain
# input data-structures
in=[q]
# output data-structures
out=[qNew]
# CFL number to use
cfl=0.5
# legacy time integration scheme, attribute; should be set to "none"
timeIntegrationScheme=none
# Riemann solver
numericalFlux=hllcEulerFlux
# Limiter to use
limiter=[muscl]
# Form of variable to limit
variableForm=primitive
# Whether to check variables for physical validitiy
preservePositivity=0
# Hyperbolic equation system
<Equation euler>
kind=eulerEqn
# Adiabatic index
gasGamma=1.6666666666666667
</Equation>
# Hyperbolic equation to solve
equations=[euler]
</Updater>
Details of the block and the meaning of each component are given in classicMusclUpdater (1d, 2d, 3d) . The Updater computes the numerical flux for the hyperbolic system. This same Updater can be used for all hyperbolic equations available in USim.
At each time step, we have to apply boundary conditions at the left and right of the domain to ensure that at the next time step, physically-valid data is used to update the conserved state. Without this, the simulation will fail. It is possible to specify arbitrary boundary conditions in USim.
For the Sod shock tube example considered here, appropriate boundary conditions are outflow (“open”) boundary conditions at both ends of the domain. This was done in Using USim to solve the Euler Equations using the macros:
boundaryCondition(copy,left)
boundaryCondition(copy,right)
These expand to yield:
<Updater copyBoundaryOnEntityleft>
kind=copy1d
onGrid=domain
in=[q]
out=[q]
entity=left
</Updater>
<Updater copyBoundaryOnEntityright>
kind=copy1d
onGrid=domain
in=[q]
out=[q]
entity=right
</Updater>
The copy boundary condition block is described in copy (1d, 2d, 3d). This boundary condition updater copies the values on the layer next to the ghost cells into the ghost cells - this is equivalent to a zero derivative boundary condition.
If we are evolving in more than one-dimension, we have to specify boundary conditions on the rest of the domain boundaries. This was done in Using USim to solve the Euler Equations using the macro:
boundaryCondition(periodic)
This expands to yield:
<Updater periodicBoundaryOnEntityghost>
kind=periodicCartBc1d
onGrid=domain
in=[q]
out=[q]
</Updater>
The periodic boundary condition updater is documented at periodicCartBc (1d, 2d, 3d)).
In order to solve the Euler equations, we have to advance the conserved quantities from time \(t\) to \(t+\Delta t\). This is done by applying a time integration scheme. In Using USim to solve the Euler Equations, we did this using the macro:
timeAdvance(TIME_ORDER)
This expands to yield:
<Updater mainIntegrator>
kind=multiUpdater1d
onGrid=domain
in=[q]
out=[qNew]
<TimeIntegrator timeStepper>
kind=rungeKutta1d
onGrid=domain
scheme=second
</TimeIntegrator>
<UpdateStep boundaryStep>
operation=boundary
updaters=[copyBoundaryOnEntityleft copyBoundaryOnEntityright periodicBoundaryOnEntityghost]
syncVars=[q]
</UpdateStep>
<UpdateStep integrationStep>
operation=integrate
updaters=[hyper]
syncVars=[qNew]
</UpdateStep>
<UpdateSequence sequence>
startOnly=[]
restoreOnly=[]
writeOnly=[]
loop=[boundaryStep integrationStep]
</UpdateSequence>
</Updater>
This is the multiUpdater (1d, 2d, 3d), which is the most powerful updater available in USim. It enables the user to combine together almost any of the Updaters available in USim and apply explicit, super-time step or fully-implicit Time Integrators. Details of each block can be found in multiUpdater (1d, 2d, 3d), however here are a few key points. The multiUpdater integrates system of equations in time. The time integration method is specified in TimeIntegrator, after that a series of UpdateSteps are defined. Each UpdateStep contains an operation type (which is optional), a list of updaters that are evaluated during the step and an optional list of syncVars which specifies which variables should synchronized after the updateStep is called. Synchronization (syncVars) are important for parallel computation and tell USim to synchronize a DataStruct across a parallel domain. Default synchronization can be specified by setting defaultParallelSync=true at the top of the file - at that point all the syncVars can be left blank. This approach is less efficient than manual synchronization, but is less prone to user error.
The final element of advancing the conserved quantities from time \(t\) to \(t+\Delta t\) is to copy the updated data in qnew to q. In USim, we accomplish this through use of a Linear Combiner:
<Updater copier>
kind=linearCombiner1d
onGrid=domain
in=[qNew]
out=[q]
coeffs=[1.0]
</Updater>
This updater block is also generated when we expand the timeAdvance(TIME_ORDER) macros and is described in linearCombiner (1d, 2d, 3d) This updater solves the equation \(Out = Coeffs*In\), where in this case \(In = q_{New}\), \(Out = q\) and \(Coeffs = 1.0\).
In Using USim to solve the Euler Equations, 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 Sod shock case in ShockTube example, the computations performed by this macro results in the following UpdateSteps and UpdateSequence being generated:
<UpdateStep startStep>
updaters=[setVar copyBoundaryOnEntityleft copyBoundaryOnEntityright periodicBoundaryOnEntityghost]
syncVars=[q]
</UpdateStep>
<UpdateStep restoreStep>
updaters=[]
syncVars=[]
</UpdateStep>
<UpdateStep bcStep>
updaters=[copyBoundaryOnEntityleft copyBoundaryOnEntityright periodicBoundaryOnEntityghost]
syncVars=[q]
</UpdateStep>
<UpdateStep hyperStep>
updaters=[mainIntegrator copier]
syncVars=[]
</UpdateStep>
<UpdateStep writeStep>
updaters=[computePressure computeVelocity computeDensity]
syncVars=[]
</UpdateStep>
Each UpdateStep block calls a number of updaters in in the order given in the updaters string list and when USim is executed in parallel, handle synchronization of data across processors through the syncVar string list. Note that if these lists are empty and the Update Step is called, then it simply returns without doing anything.
As part of the expansion of the runSimulation macro, the UpdateStep blocks are added to appropriate string lists in the UpdateSequence:
<UpdateSequence simulation>
startOnly=[generateStep startStep]
restoreOnly=[generateStep restoreStep bcStep]
writeOnly=[bcStep writeStep]
loop=[hyperStep]
</UpdateSequence>
UpdateSteps listed in the startOnly string list are executed at the beginning of the simulation and generally are used to define generate entities in the grid and definine iniital conditions. UpdateSteps listed in the restoreOnly string list are executed when a simulation is restartedand generally are used to define generate entities in the grid, ensure that data structures contain valid data and apply boundary conditions. UpdateSteps in the writeOnly string list are executed only when data is output to a file. Finally, UpdateSteps listed in the loop string list are executed at each time step.
The input file for the problem Shock Tube in the USimBase package demonstrates each of the concepts described above to evolve the classic Sod Shock tube problem in one-dimensional hydrodynamics. You can view the actual input blocks to Ulixes in the Setup window by choosing Save And Process Setup and then clicking on the shockTube.in file. In the .in file all macros are expanded to produce input blocks.
In earlier tutorials, we developed a simple pattern that could be used to design USim simulations using macros. This pattern can be repeated when we don’t use macros; however, we now need to add the grids, data structures, updaters, boundary conditions and time integration schemes by hand and then tell USim how to run them. An initial pattern looks like the following:
# Specify parameters for the specific physics problem
$ PARAM_1 = <value>
$ PARAM_2 = <value>
$ PARAM_N = <value>
# Initialize a USim simulation
# Simulation start and end times
tStart = <float>
tEnd = <float>
# Number of data files to write
numFrames = <integer>
# Initial time-step to use
initDt = <float>
# Level of feedback to give user
verbosity = <info/debug>
<Component fluids>
kind = updaterComponent
# Setup the grid
<Grid Grid_Name (type=string)>
<grid parameters>
</Grid>
# Create data structures needed for the simulation
<DataStruct DataStruct_Name1 (type=string)>
kind = nodalArray
onGrid = <Grid_Name>
numComponents = <int>
</DataStruct>
<DataStruct DataStruct_Name2 (type=string)>
kind = nodalArray
onGrid = <Grid_Name>
numComponents = <int>
</DataStruct>
<DataStruct DataStruct_NameN (type=string)>
kind = nodalArray
onGrid = <Grid_Name>
numComponents = <int>
</DataStruct>
# Specify initial condition
<Updater Initialization_Updater_Name (type=string)>
kind = initialize<NDIM>d
onGrid = <Grid_Name>
out = <DataStruct_Name for t = 0>
# initial condition to use
<Function func>
kind = exprFunc
# 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]
</Updater>
# Boundary conditions
<Updater BoundaryCondition_Name1 (type=string)>
kind=<string><NDIM>d
onGrid=<Grid_Name>
in=<DataStruct_Name for t^n>
out=<DataStruct_Name for t^n>
entity=<Entity_Name (type=string)>
<Boundary Condition Parameters>
</Updater>
# Time integration
<Updater TimeIntegrationUpdater_Name (type=string)>
kind=multiUpdater<NDIM>d
onGrid=<Grid_Name>
in=<DataStruct_Name for t^n>
out=<DataStruct_Name for t^n+1>
<TimeIntegrator timeStepper>
kind=rungeKutta<NDIM>d
onGrid=<Grid_Name>
scheme=<string>
</TimeIntegrator>
<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.