In the USimBase tutorials the basic concepts of USim were described. In this HEDP tutorial we show how to use the ideal MHD equations with general equation of state and a simple radiation model.
Contents
The plasma jets problem uses an ideal MHD model with general equation of state idealMhdEosEqn and a simple radiation loss model bremsPowerSrc. More sophisticated radiation loss can be included with purchase of the PROPACEOS tables from prism computational sciences or use of the SESAME tables from Los Alamos National Laboratory. The key to this whole development is the use of the idealMhdEosEqn in the hyperbolic updater given below:
<Updater hyper>
kind = classicMuscl$DIM$d
onGrid = domain
numericalFlux = NUMERICAL_FLUX
timeIntegrationScheme = none
variableForm = conservative
preservePositivity = true
limiter = [LIMITER, LIMITER, none]
in = [q, pressure, soundSpeed]
out = [qnew]
cfl = CFL
equations = [mhd]
<Equation mhd>
kind = idealMhdEosEqn
basementDensity = $0.1*BASEMENT_DENSITY$
basementPressure = $0.1*BASEMENT_PRESSURE$
mu0 = MU0
correctionSpeed = 0.0
</Equation>
</Updater>
A few key differences for this classicMusclUpdater (1d, 2d, 3d) updater are that it takes in multiple input variables. In this case the conserved variables q, the fluid pressure pressure and an estimate of the speed of sound soundSpeed. Since the pressure and sound speed are inputs the user can define these however they wish, with the caveat that the sound speed should be a good estimate in order to produce accurate (and stable) results.
For every variable input into the classicMusclUpdater (1d, 2d, 3d) updater a limiter must be provided limiter=[LIMITER, LIMITER, none]. In this case we limit the conserved variable and the pressure, but don’t limit the sound speed since it’s only used to add diffusion to the system and compute a time step. In this case the `correctionSpeed is being set to 0.0 as we will use a separate updater to evolve the correction potential.
The key issue then becomes how to compute the pressure and sound speed outside the hyperbolic update? As an example we use an ideal equation of state and use a combiner (1d, 2d, 3d) to compute the pressure from the energy:
<Updater computePressure>
kind = combiner$DIM$d
onGrid = domain
in = [q]
out = [pressure]
mu0 = MU0
gamma = GAMMA
bp = BASEMENT_PRESSURE
indVars_q = ["rho","mx","my","mz","en","bx","by","bz","phi"]
exprs = ["max(1.1*bp,(gamma-1)*(en-(0.5/mu0)*(bx*bx+by*by+bz*bz)-0.5*(mx*mx+my*my+mz*mz)/rho))"]
</Updater>
In this case a basement pressure is used so the maximum of the basement pressure and the computed pressure is taken. Alternatively, the propaceosVariables or sesameVariables updaters may be used to compute an equation of state for the pressure as a function of the temperature and density. Similarly, we use a combiner to compute the sound speed based on the pressure and the density:
<Updater computeSoundSpeed>
kind = combiner$DIM$d
onGrid = domain
in = [q, pressure]
out = [soundSpeed]
gamma = GAMMA
indVars_q = ["rho","mx","my","mz","en","bx","by","bz","phi"]
indVars_pressure = ["pressure"]
exprs = ["sqrt(gamma*pressure/rho)"]
</Updater>
These results can be easily applied in a multiUpdater to produce an MHD model with a user specified equation of state. The multiUpdater (1d, 2d, 3d) for this system follows where the pressure and density are computed before hyper updater is called:
<Updater rkUpdater>
kind = multiUpdater$DIM$d
onGrid = domain
timeIntegrationScheme = RKSCHEME
updaters = [bc, computePressure, computeSoundSpeed, hyper]
syncVars_bcBack =[q]
syncVars_hyper = [qnew]
integrationVariablesIn = [q]
integrationVariablesOut = [qnew]
dummyVariables_q = [dummy1, dummy2]
syncAfterSubStep = [qnew]
</Updater>
Adding an energy loss term is quite simple, except the user needs to make sure that the energy loss does not lead negative energies in the MHD model. In this case we use radiated power as estimated by analytical estimates of Bremsstrahlung radiation. The bremsPowerSrc computes radiated power density (\(W/m^{3}\)) given plasma density, temperature and the effective Z. The block below shows how to use the bremsPowerSrc:
<Updater computeRadPower>
kind = equation$DIM$d
onGrid = domain
in = [density, temperature, zEffective]
out = [radiationPower]
<Equation Bremsstrahlung>
kind = bremsPowerSrc
</Equation>
</Updater>
The dataStruct radiationPower has a single component and this value must then be subtracted from the fluid energy equation as shown below using a combiner (1d, 2d, 3d):
<Updater addRadiation>
kind = combiner$DIM$d
onGrid = domain
in = [qnew, radiationPower]
out = [qnew]
indVars_qnew = ["rho","mx","my","mz","en","bx","by","bz","phi"]
indVars_radiationPower = ["radiation"]
exprs = ["rho","mx","my","mz","en-radiation","bx","by","bz","phi"]
</Updater>
In certain cases the radiated power an be huge, so much so that if the time step is not reduced all the energy in the system will dissappear. To get around this problem we use a positiveValue time step restriction inside a timeStepRestrictionUpdater (1d, 2d, 3d). A timeStepRestrictionUpdater (1d, 2d, 3d) takes a list of input variables, and a list of restriction‘s. The restriction in this case is the positiveValue (1d, 2d, 3d) restriction. The variable positiveIndex refers to the index of the first input variable (in this case pressure). The variable sourceIndex refers to the index of the second variable sourceIndex. In this case the maximum time step is computed using
And what the positiveValue (1d, 2d, 3d) restriction does is it choses a time step so that you only subtract off the fraction alpha of the remaining energy in the input pressure. This time step is then compared with the other time steps computed by other restrictions and updaters and the smallest time allowable time step is used for the update:
<Updater timeStepRestriction>
kind = timeStepRestrictionUpdater$DIM$d
in = [pressure, radiationPower]
onGrid = domain
restrictions = [radiation]
<TimeStepRestriction radiation>
kind = positiveValue
positiveIndex = 0
sourceIndex = 0
basementValue = BASEMENT_PRESSURE
alpha = 0.25
coefficient = -1.0
</TimeStepRestriction>
</Updater>
In this case radiatedPower is always positive. Since radiatedPower is actually subtracted from pressure (not added) we set coefficient=-1.0. If radiatedPower was defined to be negative then coefficient=1.0 would be the correct value. With radiation added the multiUpdater (1d, 2d, 3d) looks like this:
<Updater rkUpdater>
kind = multiUpdater$DIM$d
onGrid = domain
timeIntegrationScheme = RKSCHEME
updaters = [correctMore, bc, computeDensity, computeTemperature, computePressure, \
computeSoundSpeed, computeRadPower, timeStepRestriction, \
hyper, addRadiation]
syncVars_bcBack =[q]
syncVars_hyper = [qnew]
integrationVariablesIn = [q]
integrationVariablesOut = [qnew]
dummyVariables_q = [dummy1, dummy2]
syncAfterSubStep = [qnew]
</Updater>
The timeStepRestriction is included in the multiUpdater (1d, 2d, 3d), but could also be included in a separate UpdateStep depending on how rapidly the power loss is changing.
In this problem divergence cleaning is computed as a separate step using hyperbolicCleanEqn. This approach has benefits because a less diffusive method can be used for the cleaning update and for certain choices of numerical flux, it can reduce overall diffusion of the solution. We use a separate classicMusclUpdater (1d, 2d, 3d) updater as defined below:
<Updater hyperClean>
kind = classicMuscl2d
timeIntegrationScheme = none
variableForm = conservative
cfl = CFL
numericalFlux = hlleFlux
limiter = [muscl]
onGrid = domain
in = [qClean]
out = [qCleanNew]
equations = [clean]
<Equation clean>
kind = hyperbolicCleanEqn
waveSpeed = CORRECTIONSPEED
</Equation>
</Updater>
Which uses the hyperbolicCleanEqn. The hyperbolicCleanEqn takes one parameter waveSpeed which is the speed that the correction wave propagates at. As with the multiUpdater (1d, 2d, 3d) for evolving the rest of the MHD model, we must also include boundary conditions for the cleaning update:
<Updater rkClean>
kind = multiUpdater2d
onGrid = domain
timeIntegrationScheme = rk2
updaters = [cleanCopyBc, hyperClean]
integrationVariablesIn = [qClean]
integrationVariablesOut = [qCleanNew]
dummyVariables_qClean = [dummyClean1,dummyClean2]
syncAfterSubStep = [qCleanNew]
</Updater>
<Updater cleanCopyBc>
kind = copy2d
onGrid = domain
out = [qClean]
entity = ghost
</Updater>
<Updater copyClean>
kind = linearCombiner2d
onGrid = domain
in = [qCleanNew]
out = [qClean]
coeffs = [1.0]
</Updater>
In the above example we used the DataStructAlias which is a pointer to a DataStruct:
<DataStructAlias qClean>
kind = nodalArray
target = q
componentRange = [5,9]
</DataStructAlias>
<DataStructAlias qCleanNew>
kind = nodalArray
target = qNew
componentRange = [5,9]
</DataStructAlias>
In the DataStructAlias we only need to define target which is the DataStruct that the alias is pointing to, the componentRange which defines which components of DataStruct to use. In the case above componentRange=[5,9] means that qClean points to elements [5,6,7,8] of q.
Finally, we put all of this into its own update step:
<UpdateStep cleanStep>
updaters = [rkClean, copyClean]
</UpdateStep>
Which can then be called before the main multiUpdater (1d, 2d, 3d).
The vertexJetUpdater (1d, 2d, 3d) was an updater specifically designed for modeling the Los Alamos Plasma Liner Experiment. It allows one to easily generate a series of directed jets with uniform initial conditions each aligned about their own axis. The vertexJetUpdater (1d, 2d, 3d) works in both 2 and 3 dimensions. The variety of options provided by the updater are described in the reference manual, however we can point out a few key things here. An arbitrary number of jets can be defined by adding additional vertices, but they must be labeled in sequential order vertex0,1,2,3,.... The vertex represents the tip of the jet and that jet will be directed towards the location specified by origin. The width of the jet defines the width about each jet will be evaluated perpendicular to the line from the vertex to the origin. The length specifies the distance from the vertex along the line from the vertex to the origin that the jet will be evaluated. The normalizedDensityFunction is used to define jets with non-uniform density profiles. The density profile is rotated into the coordinate system where X is parallel to the vector from the origin to the vertex and Y and Z are perpendicular. An example block is given below where it is used to update q:
<Updater jetSet>
kind = vertexJetUpdater$DIM$d
origin = [0.0, 0.0, 0.0]
width = JET_INIT_WIDTH
length = JET_INIT_LENGTH
radialVelocity = $-U$
numberDensity = $RHO_JET/MI$
speciesMass = MI
temperature = TKELVIN
<normalizedDensityFunction>
kind = exprFunc
preExprs = JET_DENSITY_PRE_FUNCTION
exprs = JET_DENSITY_FUNCTION
</normalizedDensityFunction>
vertex0 = [RAD, 0.0, 0.0]
vertex1 = [0.0, RAD, 0.0]
vertex2 = [-RAD, 0.0, 0.0]
vertex3 = [0.0, -RAD, 0.0]
vertex4 = [$(1.0/sqrt(2.0))*RAD$, $(1.0/sqrt(2.0))*RAD$, 0.0]
vertex5 = [$-(1.0/sqrt(2.0))*RAD$, $(1.0/sqrt(2.0))*RAD$, 0.0]
vertex6 = [$-(1.0/sqrt(2.0))*RAD$, $-(1.0/sqrt(2.0))*RAD$, 0.0]
vertex7 = [$(1.0/sqrt(2.0))*RAD$, $-(1.0/sqrt(2.0))*RAD$, 0.0]
includeElectronTemperature = false
correctionSpeed = CORRECTIONSPEED
gasGamma = GAMMA
#We use the mhdMunzEqn to initialize energy based on a temperature
model = mhdMunzEqn
mu0 = MU0
onGrid = domain
out = [q]
</Updater>
The input file for the problem Plasma Jet Merging in the USim USimHEDP package demonstrates each of the concepts described above to solve the plasma jet merging problem.