In the USimBase tutorials the basic concepts of USim were described. In this HEDP tutorial we show how to use the two-temperature mhd equations with general equation of state on an unstructured grid.
Contents
This problem is performed on an unstructured grid in axisymmetric geometry. The standard grid block is used making sure to set isRadial=true and using the unstructured gmsh grid file called dpf1.msh:
<Grid domain>
kind = unstructured
ghostLayers = 2
isRadial = true
decomposeHalos = false
<Creator ctor>
kind = gmsh
ndim = 2
file = dpf1.msh
</Creator>
</Grid>
The twoTemperatureMhdEosEqn is considerably more complex than any of the ideal MHD models since the electric field is computed externally. The hyperbolic updater is given below using the classicMusclUpdater (1d, 2d, 3d) method:
<Updater hyper>
kind = classicMuscl2d
onGrid = domain
variableForm = conservative
gradientType = leastSquares
timeIntegrationScheme = none
preservePositivity = true
numericalFlux = NUMERICALFLUX
limiterType = component
in = [q, pressure, electronPressure, soundSpeed, J, E, Z]
out = [qnew]
cfl = CFL
limiter = [LIMITER, LIMITER, LIMITER, none, LIMITER, LIMITER, none]
equations = [mhd]
sources = [axisymmetricSource, mhdSrc]
<Equation mhd>
kind = twoTemperatureMhdEosEqn
mu0 = MU0
ionMass = MI
gasGamma = GAMMA
fundamentalCharge = CHARGE
basementDensity = BASEMENTDENSITY
basementPressure = BASEMENTPRESSURE
</Equation>
<Source axisymmetricSource>
kind = mhdSym
symmetryType = cylindrical
model = twoTemperatureMhdEosEqn
inputVariables = [q, pressure, electronPressure, J, E, Z]
fundamentalCharge = CHARGE
mu0 = MU0
ionMass = MI
</Source>
<Source mhdSrc>
kind = mhdSrc
model = twoTemperatureMhdEosEqn
inputVariables = [q, J, E, Z]
fundamentalCharge = CHARGE
mu0 = MU0
ionMass = MI
</Source>
</Updater>
The updater differs from previous MHD type schemes in that we have 7 input variables which are described in the reference manual for twoTemperatureMhdEosEqn. The 7 variables are the conserved variable vector q, the pressure p the electron pressure electronPressure an estimate of the local sound speed soundSpeed the current density (computed from the curl of the magnetic field) J, the electric field E and the local charge state Z which can vary spatially
Limiters must be defined for each input variable, in this case we have
limiter = [LIMITER, LIMITER, LIMITER, none, LIMITER, LIMITER, none]
where each limiter corresponds to the corresponding input variable. Limiters need to be applied to q, p, electronPressure, E and J since derivatives of each of these terms are compute. No limiters are applied to the remaining variable. Limiters could also be applied to soundSpeed and Z however the stability of the scheme does not depend on limiting this terms so we reduce the work load by not limiting these variables.
The Equation block is kind=twoTemperatureMhdEosEqn and is given by
<Equation mhd>
kind = twoTemperatureMhdEosEqn
mu0 = MU0
ionMass = MI
gasGamma = GAMMA
fundamentalCharge = CHARGE
basementDensity = BASEMENTDENSITY
basementPressure = BASEMENTPRESSURE
</Equation>
Where the variables are described in the reference manual. In addition to the Equation block we require two Source blocks. The following line indicates which sources will be called by the updater
sources = [axisymmetricSource, mhdSrc]
The first source adds the axisymmetric terms of the hyperbolic flux so that the derivative is correct in axisymmetric geometry. This block is defined in mhdSym and some specifics outlined here
<Source axisymmetricSource>
kind = mhdSym
symmetryType = cylindrical
model = twoTemperatureMhdEosEqn
inputVariables = [q, pressure, electronPressure, J, E, Z]
fundamentalCharge = CHARGE
mu0 = MU0
ionMass = MI
</Source>
We are interested in cylindrical symmetry so symmetryType=cylindrical is specified. The model being used is model = twoTemperatureMhdEosEqn. The key parameter in the axisymmetricSource block is the inputVariable block. In computing sources, if USim does not include inputVariable block it will assume that the inputVariables will be provided for the updaters in list in the order they are provided. In this case this is not what is desired so we must specify a inputVariable list so that USim uses the correct variables in evaluating the source.
First note that in the classicaMuscl2d updater the in list is given as
in = [q, pressure, electronPressure, soundSpeed, J, E, Z]
The axisymmetric source does not need the variable soundSpeed, so, if no inputVariables list is specified in the axisymmetric source then USim will try and use the first 6 variables q,`pressure`,`electronPressure` soundSpeed,`J` and E in evaluating the source. Fortunately the size of soundSpeed is wrong (size 1 instead of the expected size 3) so USim will throw an exception and won’t run. In order to get around this problem, and to get the correct computed source we specify the inputVariables list. The inputVariables can only contain variables specified in the in list of the classicMusclUpdater2d. In this example the correct inputVariables is
inputVariables = [q, pressure, electronPressure, J, E, Z]
note that there is no soundSpeed variable in this list.
The next Source is the mhdSrc. This source is needed to add the \(E\cdot J_{e}\) term to the electron energy equation since the two-temperature mhd system cannot be written entirely in flux form. The source is desribed in mhdSrc and some details are given here
<Source mhdSrc>
kind = mhdSrc
model = twoTemperatureMhdEosEqn
inputVariables = [q, J, E, Z]
fundamentalCharge = CHARGE
mu0 = MU0
ionMass = MI
</Source>
The kind=mhdSrc is specified and the model is twoTemperatureMhdEosEqn. As with the axisymmetric source the user must specify inputVariables to override the default input variables that would be used based on the in list from the classicaMuscl2d updater. In this case the list is given as
inputVariables = [q, J, E, Z]
The electric field could be computed using a combiner, but USim also has a generalizedOhmsLaw (1d, 2d, 3d) updater which can simplify things. The generalizedOhmsLaw updater is discussed in the reference manual. It allows the user to add electron pressure gradient, resistive hall and ideal terms to the electric field. All derivatives (such as the electron pressure gradient) are computed outside of the generalizedOhmsLaw updater. The generalizedOhmsLaw used in this case is specified below
<Updater computeE>
kind = generalizedOhmsLaw2d
onGrid = domain
in = [q, J, Z]
out = [E]
hallTerm = false
fundamentalCharge = CHARGE
ionMass = MI
electronMass = ME
boltzmannConstant = KB
</Updater>
The updater takes in parameters conserved variable vector q, current density J and the charge state Z. These terms are sufficient to compute the ideal and hall terms. In order to add in resistivity and the electron pressure gradient refer to the manual. In this case hallTerm is set to false as it severly constrains the time step in many explicit simulations. The output is stored in the electric field E which has 3 components regardless of whether a generalizedOhmsLaw (1d, 2d, 3d) 1d, 2d or 3d is used.
Derivatives such as gradient, divergence and curl can be computed using two different types of algorithms in USim. The first approach is a finite difference approach which is fast, but becomes inaccurate on triangular and sufficiently skewed quadrilateral meshes. This approach currently works consistently in 1, 2 and 3 dimension. The following example shows how to compute the curl using a finite difference approach
<Updater computeJ>
kind = curl2d
onGrid = domain
inIndex = [5, 6, 7]
outIndex = [0, 1, 2]
in = [q]
out = [J]
coefficient = $1.0/MU0$
</Updater>
The updater is described in the reference manual. In this case we want the curl of the magnetic field to compute the current density J. The input variable is the conserved variable vector q and the output is the current density J. The magnetic field component stored in the conserved variable vector q occur in components [5,6,7] of the vector q so we must specify inIndex=[5,6,7] so that USim takes the correct components when computing J.
The second approach produces better results on unstructured grids and uses a least squares approach to computing the curl. This method is more accurate than the finite difference approach and is robust in 1 and 2 dimensions, but is slower than the previous method. An example block using the least squares vector (1d, 2d, 3d), with the option curl is given below
<Updater computeJ>
kind = vector2d
onGrid = domain
derivative = curl
numScalars = 1
orderAccuracy = 1
coefficient = $1.0/MU0$
numberOfInterpolationPoints = 5
in = [BPhi]
out = [J]
</Updater>
This updater expects an input vector of size 3 and produces and output vector of size 3. In this case, a second data struct BPhi was created and the magnetic field form q was copied into that variable. In addition, the parameter numberOfInterpolationPoints must be specified. The value chosen here is variable depending on the grid, but should generally be 5 or greater in 2d. Note that the complete list of options available in vector (1d, 2d, 3d) are gradient, curl and divergence.
Discussion of the pressure calculation for this system is included for clarity using a combiner (1d, 2d, 3d). Pressure in this case is total pressure or electron pressure + ion pressure. If we are using the ideal gas law (as is done in this case) the pressure must be broken into ion pressure and electron pressure separately. The combiner below does just that
<Updater computePressure>
kind = combiner2d
onGrid = domain
in = [q, electronPressure]
out = [pressure]
mi = MI
mu0 = MU0
gamma = GAMMA
k=KB
indVars_q = ["rho","mx","my","mz","en","bx","by","bz","ene"]
indVars_electronPressure = ["ene"]
exprs = ["(gamma-1)*(en-ene-(0.5/mu0)*(bx*bx+by*by+bz*bz)-0.5*(mx*mx+my*my+mz*mz)/rho)+(gamma-1.0)*ene"]
</Updater>
The boundary at the inlet is specified using a functionBc (1d, 2d, 3d). The functionBc allows you to specify the boundary condition as a function of time t and spatial variables x,`y`,`z`. This is to be contrasted with the generalBc (1d, 2d, 3d) which also allows boundaries to be specified as functions of other DataStructs. The functionBc used in this example is specified below
<Updater bcSource>
kind = functionBc2d
onGrid = domain
out = [q]
entity = source
<Function func>
kind = exprFunc
#insulator end
iEnd = $0.027-0.81e-2$
iRad = $1.5*0.27e-2$
b0 = B0
riseTime = 2.0e-6
mu0 = MU0
gamma = GAMMA
p=$BETA*P0$
rho0 = $2.0*BASEMENTDENSITY$
preExprs = ["uz=1.0e3","by = (b0*iRad/x)*sin(6.28*t/riseTime)","en=(p/(gamma-1))+0.5*rho0*uz*uz+0.5*by*by/mu0"]
exprs = ["rho0","0.0","0.0","rho0*uz","en","0.0","by","0.0","0.5*p/(gamma-1.0)"]
</Function>
</Updater>
As with all USim boundary conditions an entity must be defined to tell USim where to apply the boundary condition. In addition a Function block must be specified to define the initial conditions. The kind of this block is always exprFunc. Inside this block we specify a preExprs list to provide some work space for performing computations before the values of q are set in the boundary with exprs. The preExprs list is given below
preExprs = ["uz=1.0e3","by = (b0*iRad/x)*sin(6.28*t/riseTime)","en=(p/(gamma-1))+0.5*rho0*uz*uz+0.5*by*by/mu0"]
We’ve specified by as a sin function to mimic current rise and fall. The values defined in the preExprs can then be used in exprs
exprs = ["rho0","0.0","0.0","rho0*uz","en","0.0","by","0.0","0.5*p/(gamma-1.0)"]
In this examples magnetic field is fed into the domain using a small velocity uz and low density rho0. Using this technique we can make sure the magnetic field gets into the domain without defining a “resistive” layer near the wall. Since the model is an MHD model a finite plasma density must exist so that the Alfven wave speed does not become infinite.
The input file for the problem Dense Plasma Focus in the USimHEDP package demonstrates each of the concepts described above to evolve the dense plasma focus problem in 2D axisymmetric geometry using two temperature MHD with a general equation of state.