In this tutorial we show how to use USim to solve a problem with both an induced and imposed magnetic field using the gas dynamic form of the mhd equations.
Contents
This tutorial follows many of the same concepts as the previous tutorial Using USim to Solve MHD with General Equation of State. In this case we are solving the MHD system written in gas dynamic form. This form is significant since it is not conservative, but simplifies some considerations when there are both imposed and induces magnetic fields.
First of all we want to split the field into an imposed field generated by external coils and and induced field generated by the plasma motion. To do this we have 3 data structures, the first, backgroundB stores the field generated by a magnetic field coil
<DataStruct backgroundB>
kind = nodalArray
onGrid = domain
numComponents = 3
</DataStruct>
The second variable stores the conserved variables along with the total magnetic field, i.e., the induced field + the background field
<DataStruct qModified>
kind = nodalArray
onGrid = domain
numComponents = NUMCOMP
</DataStruct>
The 3rd variable stores the conserved variables and only the induced magnetic field
<DataStruct q>
kind = nodalArray
onGrid = domain
numComponents = NUMCOMP
</DataStruct>
The imposed magnetic field is generated at startup using a wire source wireFieldEqn and a Hyperbolic Equations updater. The Hyperbolic Equations updater works with all USim Algebraic Equations blocks. Since the problem is not axisymmetric we use a wire, but could use a coilFieldEqn to model a current carrying loop. The wireFieldEqn initialization of the background magnetic field is given below
<Updater initMagField>
kind = equation2d
onGrid = domain
in = []
out = [backgroundB]
<Equation coil>
kind = wireFieldEqn
outRange = [0, 1, 2]
point = [$COILRAD$, -0.001, 0.0]
mu0 = MU0
normal = [0.0, 0.0, 1.0]
current = $COILCURRENT$
</Equation>
</Updater>
Now the fields are split into induces q and imposed fields backgroundB. The imposed field backgroundB does not evolve in time, but does affect the flow of fluid through the \(J\times B\) force which is treated as a source in this example. Taking this into account we apply the hyperbolic update to q (which ignores the imposed field) not qModified which includes the imposed field. The classicMusclUpdater (1d, 2d, 3d) updater is used with the input file seen below
<Updater hyper>
kind = classicMuscl2d
timeIntegrationScheme = none
gradientType = leastSquares
numericalFlux = hlleFlux
variableForm = conservative
preservePositivity = true
limiter = [muscl, none, muscl, none]
onGrid = domain
in = [q, J, E, Z]
out = [qnew]
cfl = CFL
equations = [mhd]
<Equation mhd>
kind = gasDynamicMhdEqn
gasGamma = ION_GAMMA
mu0 = MU0
correctionSpeed = 0.0
ionMass = ION_MASS
chargeState = ZRATIO
fundamentalCharge = FUNDAMENTAL_CHARGE
</Equation>
</Updater>
Since this equation system does not use the general equation of state the inputs are somewhat simpler than in the case of the twoTemperatureMhdEosEqn. The inputs are the conserved variables q, the total current J, the electric field E and the charge state Z. Also, note that we have not included the mhdSrc which was included in the case of twoTemperatureMhdEosEqn. Instead we add the mhdSrc in a separate step. First we compute qModified which includes the conserved variables and both the induced and imposed field (we use a combiner (1d, 2d, 3d))
<Updater computeQMod>
kind = combiner2d
onGrid = domain
in = [q, backgroundB]
out = [qModified]
indVars_q = ["rho", "mx", "my", "mz", "en", "bx", "by", "bz", "phi", "ee"]
indVars_backgroundB = ["B0x","B0y","B0z"]
exprs = ["rho", "mx", "my", "mz", "en", "bx+B0x", "by+B0y", "bz+B0z", "phi", "ee"]
</Updater>
The current density is computed from the perturbed magnetic field only \(J=\frac{1}{\mu_{0}}\nabla \times B\) since it is already known that that the curl of the imposed field is zero everywhere inside the domain
<Updater computeJ>
kind = vector2d
onGrid = domain
derivative = curl
numScalars = 1
orderAccuracy = 2
coefficient = 1.0
numberOfInterpolationPoints = 8
in = [B]
out = [J]
</Updater>
Now that we’ve computed qModified we use this to compute the mhdSrc through the use of a Hyperbolic Equations so that both the induced and imposed fields provide a force to the fluids. The source term is stored in src.
<Updater computeSource>
kind = equation2d
in = [qModified, J, E, Z]
out = [src]
onGrid = domain
<Equation gdSrc>
kind = mhdSrc
model = gasDynamicMhdEqn
gasGamma = ION_GAMMA
electronGasGamma = ELECTRON_GAMMA
ionMass = ION_MASS
chargeState = ZRATIO
fundamentalCharge = FUNDAMENTAL_CHARGE
correctionSpeed = 0.0
mu0 = MU0
</Equation>
</Updater>
The source term is then added to the update q, qnew through the use of a uniformCombiner (1d, 2d, 3d). The uniformCombiner (1d, 2d, 3d) acts just like a combiner (1d, 2d, 3d) except that it assumes all input variables and output variables are of the same size and that the same operation is being applied to all components. In the block below we have added src generated by the computeSource Updater above
<Updater updateQ>
kind = uniformCombiner2d
onGrid = domain
in = [qnew, src]
out = [qnew]
indVars_qnew = ["qn"]
indVars_src = ["s"]
exprs = ["qn+s"]
</Updater>
As in the previous tutorial the electric field is computed using a generalizedOhmsLaw (1d, 2d, 3d) updater. Notice that the conserved variables including the total magnetic field qModified is used in computing the electric field. In addition, we’ve included a resistivity term by defining resistivity=eta. eta is a nodalArray defined using a combiner (1d, 2d, 3d)
<Updater computeE>
kind = generalizedOhmsLaw2d
onGrid = domain
in = [qModified, J, Z]
out = [E]
resistivity = eta
hallTerm = false
fundamentalCharge = FUNDAMENTAL_CHARGE
ionMass = ION_MASS
electronMass = ME
boltzmannConstant = KB
mu0 = MU0
</Updater>
USim has a Time Step Restriction for explicit diffusion type terms. The restriction can be applied to viscous, thermal diffusion and resistive terms. In order to compute the restriction properly the diffusion coefficient needs to be computed properly. In the case of resistivity we want to compute the diffusion coefficient \(\gamma\) in the proper form
It turns out gamma=eta/mu_{0} in the case of magnetic field diffusion. In this simulation we first compute the resistivity
<Updater initEta>
kind = initialize2d
onGrid = domain
out = [eta]
<Function func>
kind = exprFunc
eta0 = RESISTIVITY
exprs = ["eta0"]
</Function>
</Updater>
and then compute etaBymu0 = eta0/mu0 so that we can compute the explicit time step constraint
<Updater initEtaBymu0>
kind = combiner2d
onGrid = domain
in = [eta]
out = [etaBymu0]
indVars_eta = ["eta0"]
mu0 = $MU0$
exprs = ["eta0/mu0"]
</Updater>
etaByMu0 can then be used in the quadratic (1d, 2d, 3d) time step restriction updater
<Updater timeStepRestriction>
kind = timeStepRestrictionUpdater2d
in = [etaBymu0]
onGrid = domain
restrictions = [quadratic]
<TimeStepRestriction quadratic>
kind = quadratic
cfl = 0.5
</TimeStepRestriction>
</Updater>
Note that diffusion terms have explicit time step restriction given by \(\Delta t \approx \Delta x^{2}\) so that the time step reduces quadatically as the grid is refined. Currently we can get around this problem in USim by using super time stepping (see Advanced Time-Stepping Methods in USim).
Divergence cleaning with imposed fields is accomplished by cleaning the induced field only. It is already known that the imposed field is divergence free so that part can be ignored. Hyperbolic divergence cleaning in this case, but only applied to the perturbed field. A DataStructAlias is used to point just to the magnetic field and correction potential components in the conserved variable vector
<DataStructAlias qClean>
kind = nodalArray
target = q
componentRange = [5,9]
writeOut = 0
</DataStructAlias>
This DataStructAlias is used in the hyperbolic cleaning step which consists of a classicMusclUpdater (1d, 2d, 3d) updater, a multiUpdater (1d, 2d, 3d) and associated boundary conditions.
<Updater hyperClean>
kind = classicMuscl2d
timeIntegrationScheme = none
gradientType = leastSquares
variableForm = conservative
cfl = CFL
numericalFlux = hlleFlux
limiter = [muscl]
onGrid = domain
in = [qClean]
out = [qCleanNew]
equations = [clean]
syncAfterSubStep = [qCleanNew]
<Equation clean>
kind = hyperbolicCleanEqn
waveSpeed = CORRECTION_SPEED
</Equation>
</Updater>
<Updater hyperClean2>
kind = multiUpdater2d
onGrid = domain
timeIntegrationScheme = rk2
updaters = [cleanCopyBc, cleanInflow, hyperClean]
integrationVariablesIn = [qClean]
integrationVariablesOut = [qCleanNew]
dummyVariables_qClean = [dummyClean1,dummyClean2]
syncAfterSubStep = [qCleanNew]
</Updater>
In this case the boundary conditions consist of copying out the magnetic field and reversing the sign of the correction potential. In can be the case that field build up can occur at the boundary and in these cases simply setting the perturbed magnetic field to zero resolves the issue. Boundary conditions used for the hyperbolic cleaning step with the perturbed magnetic field are shown below using a generalBc (1d, 2d, 3d)
<Updater cleanCopyBc>
kind = generalBc2d
onGrid = domain
in = [qClean]
dynVectors = []
indVars_qClean = ["Bx","By","Bz","Phi"]
exprs = ["Bx","By","Bz","-Phi"]
out = [qClean]
entity = ghost
</Updater>
<Updater cleanInflow>
kind = generalBc2d
onGrid = domain
in = [qClean]
dynVectors = []
indVars_qClean = ["Bx","By","Bz","Phi"]
exprs = ["Bx","By","Bz","-Phi"]
out = [qClean]
entity = ghost
</Updater>
The input file for the problem Magnetic Nozzle in the USimHEDP package demonstrates each of the concepts described above.