19.1 Diffuse Unit

The Diffuse unit contains four things:

1. Flux-Based Diffusion Solvers: These are explicit diffusion solvers used to compute fluxes for specific processes, such as thermal conduction or magnetic resistivity. These fluxes can then be added to the total fluxes in the hydrodynamic solvers (see Sec:DiffuseFluxBased).
2. General Implicit Diffusion Solvers: The Diffuse unit also contains implicit diffusion solvers that are general purpose. These subroutines can be used to implicitly solve scalar diffusion problems in FLASH (see Sec:DiffuseGeneralSolver).
3. Flux Limiters: The Diffuse unit contains a general purpose subroutine which modifies diffusion coefficients to incorporate flux limiters (see Sec:DiffuseFluxLimiters).
4. Isotropic Thermal Diffusion: The Diffuse unit contains code which uses the general purpose implicit isotropic diffusion solver to solve an equation which models isotropic thermal conduction in a 3T plasma (see Sec:IsoDiffusion)
5. Anisotropic Thermal Diffusion: The Diffuse unit contains code which uses the general purpose implicit anisotropic diffusion solver to solve an equation which models anisotropic thermal conduction in a 3T plasma (see Sec:AnisoDiffusion)

19.1.1 Diffuse Flux-Based implementations

There are some important differences between other physics units and the flux-based Diffuse implementations:

• DiffuseFluxBased does not operate by modifying solution data arrays (like UNK, etc.) directly. DiffuseFluxBased modifies flux arrays instead.
• DiffuseFluxBased therefore depends on another physics unit to
  1. define and initialize the flux arrays (before DiffuseFluxBased calls add to them), and to
  2. apply the fluxes to actually update the solution data in UNK.
• DiffuseFluxBased calls that implement the diffusive effects are not made from Driver_evolveFlash as for other physics units and the DiffuseMain subunit. Rather, those DiffuseFluxBased calls need to be made from cooperating code unit that defines, initializes, and applies the flux arrays. As of FLASH4, among the provided unit implementations only the PPM implementation of the Hydro unit does this (by calls in the hy_ppm_sweep routine).
• The DiffuseFluxBased routines that implement the diffusive effects are called separately for each flux direction for each block of cells.

Other Hydro implementations, in particular the MHD implementations, currently implement some diffusive effects in their own flux-based way that does not use the DiffuseFluxBased unit.

To use DiffuseFluxBased routines of the Diffuse unit, a simulation should

• include the Diffuse unit using an option like -with-unit=physics/Diffuse/DiffuseFluxBased on the setup line, or REQUIRES physics/Diffuse/DiffuseFluxBased in a Config file;
• include a unit that makes Diffuse calls, currently the PPM implementation of Hydro;
• set useXXX runtimes parameters appropriately (see below).

The appropriate calls to DiffuseFluxBased routines will then be made by the following Hydro code (which can be found in hy_ppm_sweep.F90):  

call Diffuse_therm(sweepDir, igeom, blockList(blk), numCells,&
                        blkLimits, blkLimitsGC, primaryLeftCoord, &
                        primaryRghtCoord, tempFlx, tempAreaLeft)


call Diffuse_visc (sweepDir, igeom, blockList(blk), numCells,&
                        blkLimits, blkLimitsGC, primaryLeftCoord,primaryRghtCoord,&
                        tempFlx, tempAreaLeft,secondCoord,thirdCoord)


!!      call Diffuse_species(sweepDir, igeom, blockList(blk), numCells,&
!!                         blkLimits, blkLimitsGC, primaryLeftCoord,primaryRghtCoord,&
!!                         tempFlx, tempFly, tempFlz)

To use the Diffuse unit for heat conduction in a simulation, the runtime parameters useDiffuse and useDiffuseTherm must be set to .true.; and to use the Diffuse unit for viscosity effects in a simulation, the runtime parameters useDiffuse and useDiffuseVisc must be set to .true..

19.1.2 General Implicit Diffusion Solver

FLASH contains a general purpose implicit diffusion solver that can be used to include physics such as electron thermal conduction and multigroup diffusion. The actual code for these solvers is located in GridSolvers. Two implementations of these solvers exist in the Diffuse unit. An older and less useful version lies under the Split directory. A newer solver which uses the HYPRE library is contained in the Unsplit directory. As the name suggests, the split solver is directionally split (default: Strang) to interact with the pencil grid infrastructure of the GridSolvers/Pfft implementation of the Grid unit. The runtime parameter useDiffuse must be set to .true. to activate the general purpose implicit solver.

The typical diffusion equation that can be solved by this unit takes the form,

$$A \frac{\partial f}{\partial t} + Cf = \nabla \cdot B \nabla f + D,$$ (19.1)

where $f$ is the variable to be diffused in time, and $A$, $B$, $C$, $D$ are parameters with spatial variation. The infrastructure is provided for a general diffusion equation and can be customized to solve specific equations. The API function Diffuse_solveScalar provides a comprehensive interface to solve different types of equations.

In the unsplit solver the equation is implicitly advanced as follows,

$$A^{n} \frac{f^{n+1}-f^{n}}{\delta t} + C^{n}f^{n+1} = \theta(\nab... ...cdot B^{n} \nabla f^{n+1})+(1-\theta)(\nabla \cdot B^{n} \nabla f^{n}) + D,$$ (19.2)

In the split solver, the solution is advanced along each direction separately. i.e,

$$A^{n} \frac{f^{n+\frac{1}{3}}-f^{n}}{\delta t} + \frac{1}{3}C^{n}... ...{\partial }{\partial x}(B^{n}\frac{\partial }{\partial x}f^{n}) + \frac{1}{3} D$$ (19.3)

$$A^{n} \frac{f^{n+\frac{2}{3}}-f^{n+\frac{1}{3}}}{\delta t} + \fra... ...\partial y}(B^{n}\frac{\partial }{\partial y}f^{n+\frac{1}{3}}) + \frac{1}{3} D$$ (19.4)

$$A^{n} \frac{f^{n+1}-f^{n+\frac{2}{3}}}{\delta t} + \frac{1}{3}C^{... ...\partial z}(B^{n}\frac{\partial }{\partial z}f^{n+\frac{2}{3}}) + \frac{1}{3} D$$ (19.5)

$\theta$	Scheme
0.0	Explicit
0.5	Crank Nicholson
1.0	Backward Euler

Implicitness: It should be noted that the parameters $A$, $B$, $C$ are time lagged and evaluated at $t^{n}$. Whereas keeping with the spirit of the implicit scheme these should have been evaluated at $t^{n+1}$. Evaluating them at $t^{n+1}$ makes the problem non-linear for non-constant values of $A$, $B$ and $C$. These equations are much more complicated to solve. Although time lagged coefficient essentially linearized the non-linear problem, the trade off is the scheme is no longer completely implicit. The scheme is fully implicit only in the constant case. This could possibly impact the overall stability characteristics of the scheme. In general, implicit schemes are numerically stable for any time step, so the runtime parameter dt_diff_factor could be set to a large value to ignore the diffusion time. However, there are cases in which HYPRE will have trouble converging on a solution, and decreasing dt_diff_factor less than 1.0 may resolve the issue.

The Laplacian term is discretized using a central difference scheme as follows,

$$\nabla \cdot B \nabla f = \frac{B_{i+\frac{1}{2},j,k}(f_{i+1,j,k}-f_{i,j,k})-B_{i-\frac{1}{2},j,k}(f_{i,j,k}-f_{i-1,j,k})}{\delta x^{2}}$$

$$+\frac{B_{i,j+\frac{1}{2},k}(f_{i,j+1,k}-f_{i,j,k})-B_{i,j-\frac{1}{2},k}(f_{i,j,k}-f_{i,j-1,k})}{\delta y^{2}}$$

$$+\frac{B_{i,j,k+\frac{1}{2}}(f_{i,j,k+1}-f_{i,j,k})-B_{i,j,k-\frac{1}{2}}(f_{i,j,k}-f_{i,j,k-1})}{\delta z^{2}}$$ (19.6)

Where, $f$ and $B$ are cell centered variables, and

$$B_{i+\frac{1}{2},j,k} = \frac{B_{i+1,j,k}+B_{i,j,k}}{2},$$ (19.7)

$$B_{i+\frac{1}{2},j,k} = \frac{B_{i-1,j,k}+B_{i,j,k}}{2},$$ (19.8)

$$B_{i,j+\frac{1}{2},k} = \frac{B_{i,j+1,k}+B_{i,j,k}}{2},$$ (19.9)

$$B_{i,j-\frac{1}{2},k} = \frac{B_{i,j-1,k}+B_{i,j,k}}{2},$$ (19.10)

$$B_{i,j,k+\frac{1}{2}} = \frac{B_{i,j,k+1}+B_{i,j,k}}{2},$$ (19.11)

$$B_{i,j,k-\frac{1}{2}} = \frac{B_{i,j,k-1}+B_{i,j,k}}{2}$$ (19.12)

Paramesh: In coming up with a numerical scheme in AMR, it is important to ensure that fluxes are conserved at fine-coarse boundary. This is typically done in FLASH using the Grid_conserveFluxes routine. However this is an option which is not available to the implicit solver. Since it uses HYPRE to solve the linear system, the flux conversation should be embedded in the system being solved by HYPRE.

Looking at standard diffusion equation mentioned earlier, for simplicity we set $C = D = 0$.

So,

$$A \frac{\partial f}{\partial t} = \nabla \cdot B \nabla f,$$ (19.13)

We introduce a new variable $\vec F$ defined as,

$$\vec F = B \nabla f,$$ (19.14)

$$A \frac{\partial f}{\partial t} = \nabla \cdot \vec F,$$ (19.15)

We discretize the equation as follows,

$$A^{n} \frac{f^{n+1}-f^{n}}{\delta t} \partial v = \theta(\vec F^{n+1}\cdot \partial \vec s)+(1-\theta)(\vec F^{n} \cdot \partial \vec s)$$ (19.16)

The idea is then to compute $\vec F\cdot \partial \vec s$ on the fine-coarse boundary so as to conserve flux. In a 1D problem the $\vec F\cdot \partial \vec s$ on a fine-coarse boundary is computed as

$$(F^{L})_{i+\frac{1}{2}} = (F^{R})_{i+\frac{1}{2}} = B_{i+\frac{1}{2}} \frac {f_{i+1} - f_{i}}{\frac{\delta x_{i+1} + \delta x_{i}}{2}},$$ (19.17)

$$(F^{L})_{i-\frac{1}{2}} = (F^{R})_{i-\frac{1}{2}} = B_{i-\frac{1}{2}} \frac {f_{i} - f_{i-1}}{\frac{\delta x_{i-1} + \delta x_{i}}{2}}.$$ (19.18)

The idea is similarly extended to 2D where coarse cell flux is set to sum of computed fine cell fluxes of its neighbors.

$$F_{c} ds_{c} = F_{1} ds_{1} + F_{2} ds_{2}$$ (19.19)

Where,

$$F_{1} = B_{face} \frac {f_{1} - f_{c}}{\frac{\delta x_{1} + \delta x_{c}}{2}}$$ (19.20)

$$F_{2} = B_{face} \frac {f_{2} - f_{c}}{\frac{\delta x_{2} + \delta x_{c}}{2}}$$ (19.21)

The value of $B_{face}$ can then be evaluated from surrounding cells using paramesh interpolation.

Usage: To use the “Split” implementation of the Diffuse unit, a simulation should include the Diffuse unit using an option like -with-unit=physics/Diffuse/DiffuseMain/Split on the setup line, or add REQUIRES physics/Diffuse/DiffuseMain/Split in the Config file. To get the “Unsplit” implementation, replace Split with Unsplit. Split and Unsplit solver are mutually exclusive implementations and cannot exist in the same simulation.

Table 19.1: comparison of capabilities of Split and Unsplit diffusion.

Capability	Split	Unsplit
Paramesh	no	yes
Uniform grid	yes	yes
Cartesian 1D, 2D	yes	yes
Cartesian 3D	yes	yes
Cylindrical 1D	yes	yes
Cylindrical 2D (R-Z)	yes	yes
Spherical 1D	yes	yes

19.1.2.1 Boundary Conditions

Unlike much of the rest of FLASH, the general purpose implicit diffusion solver does not implement boundary conditions explicitly through guard cell fills. Instead the diffusion matrix is modified to directly include boundary conditions. The bcTypes and bcValues arguments to Diffuse_solveScalar are used to indicate the boundary condition to use. Each of these arguments is an array of length six corresponding to the six faces of the domain: imin, imax, jmin, imax, kmin, kmax in that order. The bcTypes argument is used to set the type of boundary condition to use. Several integer values are acceptable that are defined in the “constants.h” header file. The allowed integer values and the mathematical effect are shown in Table 19.2. When a Dirichlet (constant value) boundary conditions is required, bcTypes stores the value to use.

Table 19.2: General Implicit Diffusion Solver Boundary Conditions

Name	Integer	Effect
Neumann	GRID_PDE_BND_NEUMANN	$\nabla f \cdot \boldsymbol n = 0$
Dirichlet	GRID_PDE_BND_DIRICHLET	$f = c$
Vacuum	VACUUM	something like $\frac{1}{2} f + \frac{1}{4} \nabla f \cdot \boldsymbol n = 0$
Outstream	OUTSTREAM	free streaming radiation at outer boundary (1D spherical only)

19.1.3 Flux Limiters

The Diffuse unit contains the subroutine Diffuse_fluxLimiter which modifies a mesh variable representing a diffusion coefficient over the entire mesh to include a flux-limit. The subroutine modifies the diffusion coefficient to smoothly transition to some maximum flux, $q_{\max}$. The subroutine takes the following arguments:  

subroutine Diffuse_fluxLimiter(idcoef, ifunc, ifl, mode, blkcnt, blklst)

where:

• idcoef [integer,input] Mesh variable representing the diffusion coefficient (e.g. COND_VAR for electron thermal conduction)
• ifunc [integer,input] Mesh variable representing the function being diffused (e.g. TELE_VAR for thermal conduction)
• ifl [integer,input] Mesh variable representing the flux limit (e.g. FLLM_VAR for electron thermal conduction)
• mode [integer,input] The type of flux-limiter to apply
• blkcnt [integer,input] Number of blocks
• blklst [integer array, input] List of blocks

The subroutine modifies the mesh variable whose index is given by idcoef by applying a flux-limit. Many methods of incorporating flux-limits into the diffusion coefficient exist. Diffuse_fluxLimiter currently supports three of these: “harmonic”, “min/max”, “Larsen”, and “none”. The mode argument selects the specific method to use. The modes are integers defined in the “constants.h” header file. String representations also exist and can be converted to the integer values using the RuntimeParameters_mapStrToInt subroutine. The string representations can be used to create runtime parameters which represent a flux-limiter mode. Table 19.3 shows the integer constants used as arguments for Diffuse_fluxLimiter as well as the string representations of these constants.

Table 19.3: String and Integer Representations of Flux-Limiter Modes

Mode	Integer Constant	String Representation
None	FL_NONE	“fl_none”
Harmonic	FL_HARMONIC	“fl_harmonic”
Min/Max	FL_MINMAX	“fl_minmax”
Larsen	FL_LARSEN	“fl_larsen”
Levermore-Pomraning	FL_LEVPOM	“fl_levermorepomraning1981”

Each flux limiter mode modifies the existing diffusion coefficient to incorporate the flux limiter in a different way. If $D$ is the original diffusion coefficient and $D_{fl}$ is the flux-limited diffusion coefficient then:

• the “none” flux limiter sets:

  $$D_{fl} = D$$ (19.22)

  
  

• the “harmonic” flux limiter sets:

  $$D_{fl} = \frac{1}{ \frac{1}{D} + \frac{\vert\nabla f\vert}{q_{\max}} }$$ (19.23)

  
  

• the “Larsen” flux limiter sets:

  $$D_{fl} = \frac{1}{ \left[ \left(\frac{1}{D}\right)^2 + \left(\frac{\vert\nabla f\vert}{q_{\max}}\right)^2 \right]^{1/2}}$$ (19.24)

  
  

• the “min/max” flux limiter sets:

  $$D_{fl} = \min \left( D, \frac{q_{\max}}{\vert\nabla f\vert} \right)$$ (19.25)

  
  

where $f$ is the value of the function being diffused. The gradient of $f$ is computed using a second order central difference.

19.1.4 Isotropic Thermal Diffusion

The Diffuse unit also contains a subroutine which solves the diffusion equations for the cases of electron and ion thermal conduction. This code is in the diff_advanceTherm subroutine that is called by Diffuse. The isotropic diffusion equation for electrons has the form:

$$\rho c_{v,\mathrm{ele}} \frac{\partial T_\mathrm{ele}}{\partial t} = \nabla \cdot K_\mathrm{ele} \nabla T_\mathrm{ele}$$ (19.26)

where $\rho$ is the mass density, $c_{v,\mathrm{ele}}$ is the electron specific heat, and $K_\mathrm{ele}$ is the electron thermal conductivity. Likewise, the isotropic diffusion equation for ions is:

$$\rho c_{v,\mathrm{ion}} \frac{\partial T_\mathrm{ion}}{\partial t} = \nabla \cdot K_\mathrm{ion} \nabla T_\mathrm{ion}$$ (19.27)

where each variable is as previously described, but now for ions.

The conductivity is computed by the Conductivity unit (see Sec:conductivity for more information). This equation is solved implicitly over a time step to avoid overly restricted time step constraints. In many cases a flux-limiter is needed in situations where overly large temperature gradients are present. The flux-limit used for electron thermal conduction, $q_{\max,\mathrm{ele}}$, is defined as:

$$q_{\max,\mathrm{ele}} = \alpha_\mathrm{ele} n_\mathrm{ele} k_B T_\mathrm{ele} \sqrt{\frac{k_B T_\mathrm{ele}}{m_\mathrm{ele}}}$$ (19.28)

where:

• $\alpha_\mathrm{ele}$ is the electron conductivity flux-limiter coefficient
• $n_\mathrm{ele}$ is the electron number density
• $k_B$ is the Boltzmann constant
• $m_\mathrm{ele}$ is the electron mass

There is an analogous flux-limit for ion thermal conduction as well.

The Diffuse_solveScalar subroutine is used for the implicit solve. The use of the isotropic thermal diffusion solver is demonstrated in several simulations. See Sec:SimShafShock and Sec:LaserSlab for examples. For thermal diffusion to function, the useConductivity and useDiffuse runtime parameters must be set to .true.. The following runtime parameters control the behavior of the electron and ion thermal conduction package:

• diff_useEleCond: Set to .true. to activate electron thermal conduction (default behavior is isotropic)
• diff_useIonCond: Set to .true. to activate ion thermal conduction (default behavior is isotropic)
• diff_eleFlMode: String indicating which flux-limiter mode to use for electron conduction. Valid values are described in Sec:DiffuseFluxLimiters.
• diff_ionFlMode: String indicating which flux-limiter mode to use for ion conduction. Valid values are described in Sec:DiffuseFluxLimiters.
• diff_eleFlCoef: Sets the value of $\alpha_\mathrm{ele}$
• diff_ionFlCoef: Sets the value of $\alpha_\mathrm{ion}$
• diff_thetaImplct: A number between 0 and 1. When set to 0, an explicit treatment is used for the electron thermal diffusion solve. When set to 1, a fully implicit solve is used.
• diff_ionThetaImplct: A number between 0 and 1. When set to 0, an explicit treatment is used for the ion thermal diffusion solve. When set to 1, a fully implicit solve is used.

The runtime parameters:  

diff_eleXlBoundaryType
diff_eleXrBoundaryType
diff_eleYlBoundaryType
diff_eleYrBoundaryType
diff_eleZlBoundaryType
diff_eleZrBoundaryType


diff_ionXlBoundaryType
diff_ionXrBoundaryType
diff_ionYlBoundaryType
diff_ionYrBoundaryType
diff_ionZlBoundaryType
diff_ionZrBoundaryType

are strings which are used to set the boundary conditions for electron and ion conduction on each of the six faces of the domain. Acceptable values are:

• “DIRICHLET”, “Dirichlet”, or “dirichlet” to indicate a Dirichlet boundary condition where the appropriate $T$ is set to zero on the boundary
• “neumann”, “zero-gradient”, to indicate a Neumann boundary condition where the gradient of the appropriate $T$ is set to zero on the boundary

19.1.5 Anisotropic Thermal Diffusion

The Diffuse unit also has the capability to solve the anisotropic thermal diffusion equations for both electrons and ions. The anisotropic equations are quite different from the isotropic ones of (19.28) and (19.29) above. The general anisotropic thermal diffusion equation has the form:

$$\rho c_{v} \frac{\partial T}{\partial t} = {\bf\nabla} \cdot \big... ...{\bf\nabla} T \times {\bf b}) + K_{\wedge} ({\bf b} \times {\bf\nabla} T) \big)$$ (19.29)

where ${\bf b}$ is a unit vector in the direction of the magnetic field, $K_{\parallel}$ is the thermal conductivity in the direction of the field, $K_{\perp}$ is the thermal conductivity perpendicular to the field, and $K_{\wedge}$ is thermal conductivity perpendicular to both the field and the temperature gradient. There are different values of these conductivities for both electrons and ions, which govern the diffusion of the electron temperature and ion temperature, respectively (in a 3T plasma). The particular values used for the conductivities are calculated by the Conductivity unit (see Sec:conductivity for more information). The specifc heat $c_v$ will also have a different value for electrons and ions. The only boundary condition currently supported for anisotropic diffusion is the homogeneous Neumann condition (see 19.2).

The anisotropic equation (19.31) is solved within the Unsplit scheme of the Diffuse unit. Since it requires magnetic fields, it is also necessary for the code to be using the Unsplit Staggered Mesh (USM) MHD solver (see Sec:usm_algorithm for more details). A setup shortcut has been provided (+supportAnisoCond), which should be used to enable all functionality available for anisotropic thermal diffusion. Additional runtime parameters control the behavior of the anisotropic thermal conduction package:

• diff_anisoCondForEle: Set to .true. to activate anisotropic electron thermal conduction in 3T case (or for general anisotropic thermal diffusion in 1T case)
• diff_anisoCondForIon: Set to .true. to activate anisotropic ion thermal conduction (3T only)
• diff_useCrossCond: Set to .true. to include the $K_{\wedge}$ terms in the anisotropic diffusion equation
• gr_hypreSlopeLimType: Controls the type of slope limiter used for transverse temperature gradients. The default value is HYPRESL_MC and other acceptable values are: HYPRESL_NONE, HYPRESL_MINMOD

It is possible to have both isotropic and anisotropic thermal diffusion. For example, if diff_anisoCondForIon is set to .false., but diff_useIonCond is .true. and diff_anisoCondForEle is .true., then the ion temperature will diffuse isotropically and the electron temperature will diffuse anisotropically. Note: the runtime parameter diff_anisoCondForEle also controls anisotropic thermal diffusion in 1T applications.

Slope-limiting for transverse temperature gradients was implemented, because without it, the ${\bf b} \cdot {\bf\nabla} T$ term in (19.31) can cause heat to flow from a cold cell to a hotter one (a violation of the 2nd law of thermodynamics). In extreme cases, this can lead to negative temperatures, but gr_hypreusefloor should be set to .true. with an appropriate value for gr_hyprefloor (default values are already set to automatocally avoid negative temperatures during diffusion). If using AMR, the slope-limiter is not applied at fine/coarse boundaries, so inaccurate temperature decreases are still technically possible. The default slope-limiter is HYPRESL_MC. HYPRESL_MINMOD is more diffusive, so it can be used if inaccurate temperature decreases are problematic. To turn off this slope-limiter feature, use HYPRESL_NONE. In general, a solver such as HYPRE_GMRES should be used when using a slope-limiter, because slope-limiting can lead to an asymmetric matrix during the implicit diffusion solve done by HYPRE. The solver type is set with the runtime parameter gr_hypresolvertype, and the default value is HYPRE_GMRES when +supportAnisoCond has been used during setup. If the HYPRE_PCG solver is used, slope-limiting will need to be turned off.

19.1.6 Magnetic Diffusion

FLASH is able to solve the resistive MHD equations with contributions to the induction equation of the form

$$\frac{\partial \mathbf{B}}{\partial t} = \nabla \times \mathbf{E}$$ (19.30)

$$= \nabla \times \left (\frac{\eta_\parallel}{\epsilon_0 \omega^2_... ...a^2_{pe} \tau_{ei}} (\mathbf{b} \times \mathbf{J}) \times \mathbf{b} \right ) .$$ (19.31)

FLASH provides two methods to add resistive MHD to a simulation, either as a contribution to the unsplit Hydro fluxes, or implicitly as a an operator split update using the Unified HYPRE solver. This is controlled with the runtime parameter resistivitySolver. When set to "implicit" the implicit solver is used.

The exact form of the generalized Ohm's Law used in Eq. 19.32 is determined by the parameter resistivityForm and described in Table 19.4. Currently the "implicit" solver only uses the "parallel" form for isotropic resistivity.

Table 19.4: Options for the runtime parameter resistivityForm and the resultant contribution to the generalized Ohm's Law.

resistivityForm	$\mathbf{E}$
"parallel"	$\frac{\eta_\parallel}{\epsilon_0 \omega^2_{pe} \tau_{ei}}\mathbf{J}$
"perpendicular"	$\frac{\eta_\perp}{\epsilon_0 \omega^2_{pe} \tau_{ei}} (\mathbf{b} \times \mathbf{J}) \times \mathbf{b}$
"fullAniso"	$\frac{\eta_\parallel}{\epsilon_0 \omega^2_{pe} \tau_{ei}}\mathbf{J} + \frac{(... ...lon_0 \omega^2_{pe} \tau_{ei}} (\mathbf{b} \times \mathbf{J}) \times \mathbf{b}$

The "implicit" Magnetic Diffusion solver uses the #t#>o advance Eq. 19.32 implicitly. All the runtime parameters in Table 8.7 that are specific to this solver are prepended with "gr_uhypreMag_". Eq. 19.32 can be expressed in the form of Eq. 8.93 with the tensors

$$A_i = 1$$ (19.32)

$$B_{ijkl} = \eta \left( \delta_{il}\delta_{kj} - \delta_{ij}\delta_{kl}\right)$$ (19.33)

$$C_{ij} =0$$ (19.34)

$$D_i =0$$ (19.35)

$$E_{ijk} = 0.$$ (19.36)

In order to supress the generation of spurious $\nabla \cdot \mathbf{B}$ errors Eq. 19.32 may be augmented with a parabolic divergence cleaning term (Dedner et. al. 2002) that damps the divergence of the magnetic field,

$$\frac{\partial \mathbf{B}}{\partial t} = c^2_p \nabla \left ( \nabla \cdot \mathbf{B} \right ).$$ (19.37)

Here the constant $c_p$ is taken as the maximum value over the whole domain of $\max(\eta, c_f)$, where $c_f$ is the fast magnetosonic speed. This can be expressed as a $B_{ijkl}$ tensor

$$B_{ijkl} = c^2_p\delta_{ik}\delta_{jl}.$$ (19.38)

In cylindrical geometries the evolution equations for the radial and axial magnetic field components may still be expressed in conservative form and so we only need the $B_{ijkl}$ tensor in Eqs 19.34 for $i=r,z$. The azimuthal field equation though will introduce some geometric source terms. These terms can be expressed in the general form of Eq. 8.98, but in particular for problems with spatially varying resistivity, $\eta$, have been observed to introduce numerical artifacts and instability, especially at fine-coarse boundaries. The conservation form may be restored by evolving $\frac{B_\phi}{r}$, which will introduce an additional $E$-factor term,

$$\partial_t \frac{B_\phi}{r} \sim \frac{1}{r}\partial_r \left [ r\eta \left ( \partial_r \frac{B_\phi}{r} + \frac{2}{r}\frac{B_\phi}{r} \right ) \right ]$$ (19.39)

$$E_{\phi\phi r} = 2\eta.$$ (19.40)

This form of the equations will be evolved when the grid geometry is "cylindrical" and the runtime parameter diff_conserveAngfield = .TRUE.. This should not be confused with the Hydro runtime parameter conserveAngField, which serves a similar purpose in the MHD solver.

While the above described conservative formulation works well for evoling the magnetic field in cylindrical coordinates, for problems with orders of magnitude difference in resistivity across material interfaces, as is typical in Z-pinch simulations the resultant fields can cause some small scale discrepencies in the calculated ohmic heating at fine-coarse boundaries that can be significantly detrimental to the overall simulation.

To overcome this the finite-volume discretization of Sec. 8.10.4.2 needs to be modified so that rather than the flux of $\frac{B_\phi}{r}$ being conserved at fine-coarse boundaries the current density $\mathbf{J}=\nabla \times \mathbf{B}$ is conserved.

$$\partial_t rB_\phi \sim r \left [ \frac{\eta}{r} \partial_r (rB_\phi) \right ]$$ (19.41)

$$\int_{V_\alpha}\partial_t rB_\phi dV \approx \frac{V_\alpha r_\alpha}{\Delta r}\left [ \frac{\eta_{\al... ...r_{\alpha-1/2}} \left . \partial_r(rB_\phi) \right \vert _{\alpha-1/2} \right ]$$ (19.42)

The discretization in Eq. 19.44 is the default for cylindrical geometry.

The boundary conditions imposed on the HYPRE matrix are determined by the the ab_boundary_type runtime parameters Table.

Table 19.5: Unified HYPRE solver boundary conditions used for magnetic diffusion

ab_boundary_type	Description
periodic	Periodic (`wrap-around')
reflect,reflecting	the normal vector components change sign
outflow	Zero-gradient boundary condition in normal direction
circuit	Dirichlet boundary condition taken from the Circuit unit
user	Dirichlet boundary condition

19.1.6.1 Ohmic Heating

The Ohmic heating from magnetic diffusion when using the "implicit" solver is computed from the updated magnetic fields using the implicit $\theta$ parameter,

$$\partial_t \rho_i \epsilon_i = (1-\theta)Q_i^n + \theta Q_i^{n+1} .$$ (19.43)

In the above $Q_i^n$ is the rate of heat added to the cell $i$ at time $t^n$. The solver supports two methods for computing $Q$ depending on the runtime parameter ${\tt diff_useQOhmEdgeMethod}$. When set to .FALSE. $Q$ is calculated as a direct centered difference of

$$Q^n = \nabla \times \left (\eta \nabla \times \mathbf{B}\right ).$$ (19.44)

In Z-pinch simulations where an artificially high vacuum resistivity is used this discretization this can bias heating into the vacuum at a conductor-vacuum interface since the model doesn't take into account the relative resistence of neighboring cells. When diff_useQohmEdgeMethod = .TRUE. the solver will use a parallel resistor model at cell edges to distribute the current density at an edge to the neighboring cells. A diagram of this model is shown in Fig. 19.2. The equivalent resistence is determined from the resistivity and the geometry of the cell quadrants that meet at the edge. The current density is calculated from finite-differences at the edge center and the total current through the cross-sectional area made up of the joining quadrants. The parallel resistor is solved for the current through each quadrant where the dissipated energy is taken as the ohmic heating in the quadrant. The heating in the entire cell is summed over the contributions from all edges on the cell.

Figure: 19.2(a) cell quadrants that share the edge at $i+\frac {1}{2},j+\frac {1}{2}$ with current density $J_{i+1/2,j+1/2}$. 19.2(b) The parallel resistor model used to calculate the Ohmic heating in the quadrants of 19.2(a)

[a] [b]