In the USimBase tutorials the basic concepts of USim were described. In this first HEDP tutorial we describe how to solve the fully electromagnetic two-fluid plasma equations using a semi-implicit operator to step over the plasma frequency and cyclotron frequency and electric field diffusion to minimize errors in in the electric field divergence relation.
Contents
To demonstrate how to use USim to solve a problem using the two-fluid plasma system, we will use the well known GEM (geomagnetic environmental modeling) reconnection challenge setup to solve fast reconnection of a current layer. The GEM challenge was originally described in
Birn, J., et al. "Geospace Environmental Modeling (GEM) magnetic reconnection challenge."
Journal of Geophysical Research: Space Physics (1978–2012) 106.A3 (2001): 3715-3719.
This tutorial is based on the GEM Challenge template.
The first thing we need to model two-fluids, is data structures for the electrons, ions and the electromagnetic field. In this case the electrons will be represented by a 5-moment compressible fluid:
<DataStruct electrons>
kind = nodalArray
onGrid = domain
numComponents = 5
</DataStruct>
The ions are also represented by a 5-moment compressible fluid:
<DataStruct ions>
kind = nodalArray
onGrid = domain
numComponents = 5
</DataStruct>
The electromagnetic field variable contains the full Field vector [Ex,Ey,Ez,Bx,By,Bz,Ep,Bp] with the variables Ep and Bp the error correction potentials. As such the electromagnetic field data structure is defined as:
<DataStruct em>
kind = nodalArray
onGrid = domain
numComponents = 8
</DataStruct>
Separate initialization using an initialize (1d, 2d, 3d) updater is performed for each variable, electrons, ions and em:
<Updater initElectrons>
kind = initialize2d
onGrid = domain
out = [electrons]
<Function func>
kind = exprFunc
.
.
.
preExprs = [ \
"rhoe = n0*me*(1.0/(cosh(y/lambda)*cosh(y/lambda))+0.2)", \
"mze = -(me/qe)*b0*(1.0/lambda)*1.0/(cosh(y/lambda)*cosh(y/lambda))", \
"ee = (1.0/12.0)*(1./(gamma-1.0))*b0*b0*(rhoe/me)+0.5*(mze*mze/rhoe)"]
exprs = [ \
"rhoe", "0.0", "0.0", "mze", "ee"]
</Function>
</Updater>
<Updater initIons>
kind = initialize2d
onGrid = domain
out = [ions]
<Function func>
kind = exprFunc
.
.
.
preExprs = [ \
"rhoi = n0*mi*(1.0/(cosh(y/lambda)*cosh(y/lambda))+0.2)", \
"ei = (5.0/12.0)*(1.0/(gamma-1.0))*b0*b0*(rhoi/mi)"]
exprs = ["rhoi", "0.0", "0.0", "0.0", "ei"]
</Function>
</Updater>
<Updater initEm>
kind = initialize2d
onGrid = domain
out = [em]
<Function func>
kind = exprFunc
.
.
.
preExprs = ["bx = b0*tanh(y/lambda)-psi*(pi/Ly)*cos(2.0*pi*x/Lx)*sin(pi*y/Ly)", \
"by = psi*(2.0*pi/Lx)*sin(2.0*pi*x/Lx)*cos(pi*y/Ly)"]
exprs = ["0.0", "0.0", "0.0", "bx", "by", "0.0", "0.0", "0.0"]
</Function>
</Updater>
In the above initialization some variables have been eliminated for conciseness. In addition, every variable, electrons, ions, em, must have their own classicMusclUpdater (1d, 2d, 3d)
<Updater hyperElectrons>
kind = classicMuscl2d
onGrid = domain
timeIntegrationScheme = none
numericalFlux = FLUID_NUMERICAL_FLUX
preservePositivity = true
limiter = [LIMITER, none]
limiterType = characteristic
variableForm = conservative
in = [electrons]
out = [electronsNew]
cfl = CFL
equations = [euler]
<Equation euler>
kind = eulerEqn
gasGamma = GAS_GAMMA
basementDensity = BASEMENT_DENSITY
basementPressure = BASEMENT_PRESSURE
</Equation>
</Updater>
<Updater hyperIons>
kind = classicMuscl2d
onGrid = domain
timeIntegrationScheme = none
numericalFlux = FLUID_NUMERICAL_FLUX
preservePositivity = true
limiter = [LIMITER]
limiterType = characteristic
variableForm = conservative
in = [ions]
out = [ionsNew]
cfl = CFL
equations = [euler]
<Equation euler>
kind = eulerEqn
basementDensity = BASEMENT_DENSITY
basementPressure = BASEMENT_PRESSURE
gasGamma = GAS_GAMMA
</Equation>
</Updater>
<Updater hyperEm>
kind = classicMuscl2d
onGrid = domain
timeIntegrationScheme = none
numericalFlux = fWaveFlux
limiterType = characteristic
limiter = [LIMITER]
variableForm = conservative
in = [em]
out = [emNew]
cfl = CFL
equations = [maxwell]
<Equation maxwell>
kind = maxwellEqn
c0 = SPEED_OF_LIGHT
gamma = BP
chi = 0.0
</Equation>
</Updater>
The coupling between the fields and fluids is provided by Lorentz forces (for the fluid equations) and current sources (for the electromagnetic field). One option is to simply add these to the right hand side of the flux calculation and then integrate, however, this leads to an explicit scheme where the plasma frequency and cyclotron frequency must be resolved. Instead we use a semi-implicit operator as defined in
Harish Kumar and Siddhartha Mishra. “Entropy Stable Numerical Schemes for Two-Fluid Plasma Equations.” Journal of Scientific Computing (2012): 1-25.
The implicit operator twoFluidSrc is a source in USim, it’s applied cell by cell and does not require a global implicit solve. The twoFluidSrc is evaluated using the explicit solution to the electron, ion and em variables and the resulting matrix is multiplied by those same variables to produce the implicit source evaluation with explicit hyperbolic terms. The semi-implicit operator is given below:
<Updater twoFluidLorentz>
kind = equation2d
onGrid = domain
in = [electronsNew, ionsNew, emNew]
out = [electronsNew, ionsNew, emNew]
<Equation twofluidLorentz>
kind = twoFluidSrc
type = 5Moment
electronCharge = ELECTRON_CHARGE
electronMass = ELECTRON_MASS
ionCharge = ION_CHARGE
ionMass = ION_MASS
epsilon0 = EPSILON0
</Equation>
</Updater>
The semi-implicit operator is applied in a special location in the multiUpdater (1d, 2d, 3d). The list of updaters in the multiUpdater (1d, 2d, 3d) define the explicit steps in the multiUpdater (1d, 2d, 3d). A second list of updaters is in the operator list. Updaters in the operator list are applied to integrationVariablesOut after a complete update. In the block below operators = [twoFluidLorentz] applies the operator after the explicit right hand side has been calculated, for example
We want to solve the hyperbolic part of the multi-fluid equations explicitly and the source term implicitly. For a first order scheme the discretization becomes.
and therefore
The term in parentheses is the twoFluidLorentz operator defined in the multiUpdater below. As stated previously, the updaters in operate are applied to the right hand side which is computed in the updaters list
<Updater rkUpdaterMain>
kind = multiUpdater2d
onGrid = domain
in = [em, ions, electrons]
out = [emNew, ionsNew, electronsNew]
<TimeIntegrator rkIntegrator>
kind = rungeKutta2d
ongrid = domain
scheme = RK_SCHEME
</TimeIntegrator>
<UpdateSequence sequence>
loop = [boundaries,hyper,implicit]
</UpdateSequence>
<UpdateStep boundaries>
updaters = [periodicEm, periodicIons, periodicElectrons, electronBcTop, electronBcBottom, \
ionBcTop, ionBcBottom, emBcTop, emBcBottom]
</UpdateStep>
<UpdateStep hyper>
operation = "integrate"
updaters = [hyperIons, hyperElectrons, hyperEm, addSource]
</UpdateStep>
<UpdateStep implicit>
operation = "operate"
updaters = [twoFluidLorentz]
</UpdateStep>
</Updater>
In the multiUpdater (1d, 2d, 3d) above we include 3 sets of in variables, 3 sets of out variables for each of the integration variables.
The results described above are sufficient to create an algorithm that steps over plasma frequency and cyclotron frequency, but this does not show us how to minimze errors in the divergence equations.
The standard approach to divergence preservation in USim is to use hyperbolic divergence cleaning. Hyperbolic divergence cleaning is described for the MHD equation in
Andreas Dedner, et al. “Hyperbolic divergence cleaning for the MHD equations.” Journal of Computational Physics 175.2 (2002): 645-673.
And for Maxwell’s equation in
Munz, C-D., et al. “Divergence correction techniques for Maxwell solvers based on a hyperbolic model.” Journal of Computational Physics 161.2 (2000): 484-511.
Unfortunately, for the two-fluid system hyperbolic cleaning of the electric field is often inadequate, or results in larger errors than we started with. Instead we use electric field diffusion. In this section we describe how we perform the divergence cleaning in the GEM challenge problem.
First off, the magnetic field can be cleaned using the hyperbolic approach. In the hyperbolic updated of the electromagnetic field the Hyperbolic Equations block defines the speed of light c0 as well as the wave speeds of the correction potentials. gamma is the magnetic field correction potential factor and chi is the electric field correction factor. The speed of the magnetic field correction potential is gamma*c0 and of the electric field correction potential chi*c0. In the case below we give gamma a finite value (typically 1.0) and we set chi to 0 so that we can use a different correction method for the electric field:
<Updater hyperEm>
kind = classicMuscl2d
onGrid = domain
timeIntegrationScheme = none
numericalFlux = fWaveFlux
limiterType = characteristic
limiter = [LIMITER]
variableForm = conservative
in = [em]
out = [emNew]
cfl = CFL
equations = [maxwell]
<Equation maxwell>
kind = maxwellEqn
c0 = SPEED_OF_LIGHT
gamma = BP
chi = 0.0
</Equation>
</Updater>
Electric field diffusion is quite simple. We add a diffusion term to the electric field given as
where \(\gamma\) is the electric field diffusion coefficient. You can see that the diffusion term never kicks in unless there is a numerical error gauss’ law. How do we go about implementing this in USim? First of all we define a set of data structures just for electric field divergence cleaning
The first data structure stores the characteristic cell length for the diffusion coefficient:
<DataStruct cellDx>
kind = nodalArray
onGrid = domain
numComponents = 1
writeOut = 0
</DataStruct>
We define a data structure for storing the error computed from \(\nabla\cdot E-\frac{\rho_{c}}{\epsilon_{0}}\):
<DataStruct residual>
kind = nodalArray
onGrid = domain
numComponents = 1
writeOut = 0
</DataStruct>
We define a data structure for storing the divergence of E:
<DataStruct divE>
kind = nodalArray
onGrid = domain
numComponents = 1
writeOut = 0
</DataStruct>
We then store the actual diffusion term in the last data structure:
<DataStruct gradDiv>
kind = nodalArray
onGrid = domain
numComponents = 3
writeOut = 0
</DataStruct>
Along with the data structures, we have a series of updaters that are used to fill up the data structures. The first updater characteristicCellLength (1d, 2d, 3d) simply computes a characteristic length for each cell. This updater only needs to be called at startup since the cell length does not change.
<Updater computeCellDx>
kind = characteristicCellLength2d
onGrid = domain
out = [cellDx]
coefficient = 1.0
</Updater>
The next updater computes the \(\nabla\cdot E\) from the electric field. It takes as in the array em. The vectorDivergence operator assumes the 3 vector of interest occurs in the 3 components of em which happens to be correct in this case as those components correspond to Ex,`Ey`,`Ez`. In addition, a parameter coeffs is provided which multiplies the resulting divergence by the factor 1.0. This is a simple way to reverse the sign of the divergence or multiply by some other factor:
<Updater computeDivE>
kind = vectorDivergence2d
onGrid = domain
in = [em]
out = [divE]
coeffs = [1.0]
</Updater>
The next updater uses an Hyperbolic Equations to compute the residual \(\nabla\cdot E-\rho_{c}/\epsilon_{0}\). The source computeChargeError expects as input \(\nabla\cdot E\) and then expectes the remaining variables to be fluid mass densities, any number of species can be added. Along with the species mass densities we provide a list indicating the species mass and species charge for each of these species densities. We also specify the permittivity epsilon0 so that the residual can be calculated.
<Updater computeResidual>
kind = equation2d
onGrid = domain
in = [divE, electrons, ions]
out = [residual]
<Equation>
kind = computeChargeError
speciesCharge = [ELECTRON_CHARGE, ION_CHARGE]
speciesMass = [ELECTRON_MASS, ION_MASS]
epsilon0 = EPSILON0
</Equation>
</Updater>
The final diffusion term including gamma factor is computed using a scalar gradient calculator. The scalar gradient takes 2 inputs. It takes the gradient of the first input (residual in this case) and multiplies that gradient by the second input cellDx. The result is the multiplied by coefficient which is constant for all space. This particular form of diffusion is such that the diffusion is stable to explicit time stepping. If super time stepping is used or subcycling, the diffusion coefficient can be increased to do a better job of error cleaning.
<Updater gradient>
kind = scalarGradient2d
gradientType = leastSquares
onGrid = domain
in = [residual, cellDx]
out = [gradDiv]
coefficient = 0.5
</Updater>
Once the diffusion term is computed it needs to be added to the right hand side of the update equation. The term is added after the hyperbolic explicit terms are computed. We use a combiner2d to add the diffusion term to the equation. In this example the updater takes two input data structures, emNew and gradDiv and dumps the output into emNew. Every input variable requires an indVars_inputName block which provides a way to access each component of the input variable. These names can then be used in the output expression exprs. In the multiUpdater (1d, 2d, 3d), the addSource block is called after the hyperbolic terms so that it is not overwritten by updaters that are called earlier.
<Updater addSource>
kind = combiner2d
onGrid = domain
in = [emNew, gradDiv]
out = [emNew]
indVars_emNew = ["ex","ey","ez","bx","by","bz","phiE","phiB"]
indVars_gradDiv = ["dx","dy","dz"]
c = SPEED_OF_LIGHT
exprs = ["ex+c*dx","ey+c*dy","ez+c*dz","bx","by","bz","phiE","phiB"]
</Updater>
Finally, boundary conditions must be provided to the residual. In this problem the residual on the boundary should be 0 so we use a functionBc (1d, 2d, 3d) to explicitly set the residual on the boundary to zero. The variable out specifies the data that the boundary condition will be applied to, while entity tells the boundary condition which boundary it should be applied to. In this case entity=ghost means that the boundary is applied to all boundary cells.
<Updater residualBc>
kind = functionBc2d
onGrid = domain
<Function func>
kind = exprFunc
exprs = ["0.0"]
</Function>
out = [residual]
entity = ghost
</Updater>
This simulation has periodic boundaries in the x direction so we also apply a periodic boundary condition to the residual. This boundary condition is called after residualBc so it overwrites the boundry conditions in the X direction set by residualBc
<Updater periodicResidual>
kind = periodicCartBc2d
onGrid = domain
in = [residual]
out = [residual]
</Updater>
Now that we’ve added a bunch of new updaters we need to modify the multiUpdater to include the changes that have been made. The new updater added to the updater list include computeDivE, computeResidual, residualBc,`periodicResidual`,`gradient` in the UpdateStep‘s compute and clean. These updaters are evaluated after the boundary conditions for electrons,`ions` and em are applied, but before the hyperbolic updaters are called hyperIons, hyperElectrons and hyperEm. Once again, for parallel simulations it’s very important to get the synchronization correct. We’ve added one synchronizations. The new synchronization occurs directly after periodicResidual. periodicResidual is the last updater that is applied before a derivative gradient is computed on residual so we need to synchronize residual at this point.
<Updater rkUpdaterMain>
kind = multiUpdater2d
onGrid = domain
in = [em, ions, electrons]
out = [emNew, ionsNew, electronsNew]
<TimeIntegrator rkIntegrator>
kind = rungeKutta2d
ongrid = domain
scheme = RK_SCHEME
</TimeIntegrator>
<UpdateSequence sequence>
loop = [boundaries,compute,clean,hyper,implicit]
</UpdateSequence>
<UpdateStep boundaries>
updaters = [periodicEm, periodicIons, periodicElectrons, electronBcTop, electronBcBottom, \
ionBcTop, ionBcBottom, emBcTop, emBcBottom]
</UpdateStep>
<UpdateStep compute>
updaters = [ computeDivE, computeResidual, residualBc, periodicResidual]
</UpdateStep>
<UpdateStep clean>
updaters = [ gradient]
</UpdateStep>
<UpdateStep hyper>
operation = "integrate"
updaters = [hyperIons, hyperElectrons, hyperEm, addSource]
</UpdateStep>
<UpdateStep implicit>
operation = "operate"
updaters = [twoFluidLorentz]
</UpdateStep>
</Updater>
Finally, recall that the characteristic cell lengths only need to be calculated once. As a result we calculate them during the initialization step.
<UpdateStep initStep>
updaters = [initElectrons, initIons, initEm, computeCellDx]
syncVars = [electrons, ions, em]
</UpdateStep>
In the GEM Challenge simulation we’ve added passive diagnostic to compute the line integral across the domain. In order to compute the line integral we add two data structures. The first data structure is called a bin. The bin is a uniform grid that overlays the exiting grid. The extents of the bin match the extents of the grid, but a bin is rectangular regardless of the shape of the grid. The bin has two important parameters, the first is onGrid which specifies which grid the bin will use to define itself, the second is the scale. scale tells roughly how many bins there are per grid cell in domain. The bin is used for fast lookup so that USim can quickly tell which cell a particular point in space is in. The bin is given as
<DataStruct cellBin>
kind = bin
onGrid = domain
scale = 2.0
</DataStruct>
In addition we need to fill the bin with data. In this particular case we want to fill each bin with a a list of indexes corresponding to the cells that fall inside each cell of the bin. The binCells (1d, 2d, 3d) updater does exactly that. This updater only needs to be called at startup since the grid does not change with time.
<Updater fillBin>
kind = binCells2d
onGrid = domain
out = [cellBin]
</Updater>
The second variable is a dynVector. A dynVector is a dataStructure which is a vector and has the same value on all domains in parallel simulations. dynVectors store a time series of data, so the value of the dynVector is dumped at every time step so it can be used to store integrated quantities. The dynVector had numComponents to specify how long the vector is. In this case it stores only 1 value. writeOut=true is set so the dynVector is written.
<DataStruct integratedFlux>
kind = dynVector
numComponents = 1
writeOut = 1
</DataStruct>
These two data structures are then used to compute the line integral with the lineIntegral (1d, 2d, 3d) updater. The lineIntegral (1d, 2d, 3d) behaves like a combiner (1d, 2d, 3d) except that it takes in nodalArray‘s and fills up a dynVector.
<Updater computeLineIntegral>
kind = lineIntegral2d
onGrid = domain
startPosition = [XMIN, 0.0]
endPosition = [XMAX, 0.0]
numberOfSamples = 1000
layout = [cellBin]
in = [em]
indVars_em = ["ex","ey","ez","bx","by","bz","phiE","phiB"]
exprs = ["0.5*abs(by)"]
out = [integratedFlux]
</Updater>
This updater should not be called within an rkUpdater otherwise it will store multiple values per time step. Instead we call this updater in its own UpdateStep
<UpdateStep diagnosticStep>
updaters = [computeLineIntegral]
</UpdateStep>
The input file for the problem GEM Challenge in the USim USimHEDP package demonstrates each of the concepts described above to solve fast magnetic reconnection.