===================================

Viscosity and Resistivity in AstroBEAR

===================================

Viscosity

The viscosity term in the Navier-Stokes Equation can be separated as:

dv/dt = \nabla \cdot T
For the applications we are considering, the fluid is isotropic. Under this assumption, the tensor can be expressed as the following combination of two scalar dynamic viscosities \mu and \mu':

T = 2 \mu \nabla v + \mu' \nabla \cdot v I

where I is the unit tensor.
It can be proven that the stress tensor should have zero trace, which leads to the following relation between \mu and \mu':

\mu' = -2/3\mu

Combining these relations, we can write down the final form of the separated viscous equation:

dv/dt=\nabla \cdot (2 \mu \nabla v) - 2/3 \nabla (\mu \nabla \cdot v)

Notice that this equation already counts the effect of all the elements in a viscosity tensor, assuming the fluid itself is isotropic.

Resistivity

The resistivity equation can be written as:
\frac{\partial \textbf{B}}{\partial t}=\nabla \times (\eta \nabla \times \textbf{B})

Here, the magnetic field are located at face centers. Using the face centered magnetic field, we can calculate the edge centered currents, which then give a source term.
The curl of the edge centered source term gives a face centered \delta B, which can then be applied to the face centered magnetic field.

===================================

Implementation

===================================

Viscosity

The viscous source term is calculated rigorously by first finding the tensor elements at the surrounding face centers of a cell center considered: The \mu sources 2\mu \nabla v are naturally located at the cell centers, the \mu' sources however needs an averaging to obtain the cell centered velocity divergence.
Therefore the scheme requires a 20 point stencil in 3-D configuration, which makes the implicit scheme a bit complicated. The following graph shows such a stencil in 2-D setting:



Using this stencil setting, one can get face centered velocity divergence for the cell considered, and then calculate the gradient components located at the cell center.
In the viscous plasma context, the physical viscosity can be calculated rigorously. The plasma viscosity is proportional to the product of ion mean free path and the ion thermal energy:

\mu \sim \lambda_{ifp}v_{thermal}

In our code, we offer the option to calculate the realistic viscosity using the following equation:

\mu = \frac{m_{i}}{3 \sqrt{2} \pi d_{i}^2} (\frac{8kT}{\pi m_{i}})^{\frac{1}{2}}

The ion mass and local temperature should be readily available, but a good estimate of ion diameter d_{i} usually requires look up a ionic radius table.
As an estimate, one can use that of hydrogen, about 0.1~0.2 nm.
One can see that the viscosity itself is anisotropic and depends on the local temperature.
That is why we have both the \mu and \mu' terms in our viscosity equation, and have to calculate face centered source terms before applying the final divergence.

The viscosity interface is as follows:
The viscosity source term can be turned on in astrobear by putting the viscosity flag to true:
lViscous = .true.
The second parameter VisType determines which type of viscosity one requires: 1 for constant viscosity; 2 for real physical viscosity.
The third parameter viscosity is the user defined constant viscosity.
The viscosity module is tested to be working compatibly with implicit thermal conduction module and uniform gravity module.

Resistivity

Implicit Implementation of Resistivity

The resistive equation can be written as:

\frac{\partial \textbf{B}}{\partial t}=\nabla \times (\eta \nabla \times \textbf{B})

The magnetic diffusivity eta is a function of temperature as we all know.

Now for the first equation, what if we want to include the temperature dependence into the code? From the next section, you can see that when a field configuration in equilibrium is subject to strong diffusion, usually heating would occur and surppress the local resistivity and thus the diffusivity. By expanding the first equation, we have the form:

\frac{\partial \textbf{B}}{\partial t}=\eta \nabla^2 \textbf{B} + \nabla \eta \times (\nabla \times \textbf{B})

Now the second term depends on all of the three components of B. So we end up with equations in which the time variance of Bx, By and Bz depend on each other.

So the coefficient array is thus no longer a tri-diagonal matrix. For cases where resistive speed is slow (which is usually the case), we can use explicit solver to treat the problem instead.

One may wonder if it is possible to throw away the second term on the right hand side of the diffusion equation to just let the diffusivity to vary with position but ignore its own spacial variation. This would give us a form of resistive MHD similar to that of the thermal diffusion case but with temperature dependence built in.

\frac{\partial \textbf{B}}{\partial t}=\eta (T) \nabla^2 \textbf{B}

Unfortunately, this does not work because the diffusion equation itself has to be divergence free. Our roughest approximation, treating the resistivity as a constant satisfies the requirement as long as the divergence and the laplacian are commutable. If we try to do the same thing to the above equation, we end up getting:

\frac{\partial \textbf{B}}{\partial t}=\eta \nabla^2 (\nabla \cdot \textbf{B}) + (\nabla \eta \cdot \nabla)\textbf{B}

The first term on the right hand side is zero but the second term is not, especially at sharp temperature fronts where grad(T) is large.
So this approximation only works under slowly varying temperature situation.

Upon a closer look at the simulations conducted below, we can observe that the energy inside the domain does not conserve. There is a small increase of total energy during the evolution, especially for the force-free field case. This phenomenon is explained below.
In the case of resistive MHD, the energy can be dissipated in the form of Joule heat, comparing to the infinite conductivity case, where the voltage inside the fluid is everywhere zero, and no heat is generated by the current.
If we dot the resistive induction equation with the magnetic field B, we obtain the time evolution equation for magnetic energy:

\partial (B^2)/\partial t + \nabla \cdot \textbf{S} = - \textbf{J}\cdot\textbf{E} = -j^2/\eta

where \textbf{S}=\textbf{J} \times \textbf{B} is the magnetic energy flux caused by resistive diffusion and j=|\textbf{J}| is the magnitude of the diffusive current.

In this equation, the \textbf{S} term accounts for the redistribution of magnetic energy (and thus the redistribution of total energy), and the j^2/\eta term accounts for the loss of magnetic energy due to reconnection.
The total energy change for the resistive step is therefore:

\partial \epsilon/\partial t + \nabla \cdot \textbf{S} = 0

Here the j^2/\eta dissipation term is absent because the dissipation of magnetic energy does not change the total energy: the loss of magnetic energy is converted into thermal energy (j^2/\eta is indeed the heat generated by current j in the plasma).
In the code, the \textbf{S} term is calculated explicitly by calculating \textbf{J} \times \textbf{B} at each face centers. This additional flux is added to the total energy flux, while the j^2/\eta dissipation term is automatically accounted.

Here I overview two approaches to solve the resistive MHD equations.

(1)Solving the decomposed format with divergence cleaning (Implicit)
This is the simplest method in terms of solving scheme. This method is based on the simplification that in some cases eta is a constant in the region of interest. This will simplify the above equation into the following form:

\frac{\partial \textbf{B}}{\partial t}=\eta \nabla^2 \textbf{B}

This equation holds true for each separate field component, and the form resembles the linear isotropic thermal conduction equation. We can solve the elliptic equation for three times or we can write a expanded linear system so that the solver is only called once:

In the constant diffusivity case, the matrix M are identical for the Bx, By and Bz components. Define C=[-6k,k,k,k,k,k,k], where k = \eta dt/dx^2. One can immediately recognize that k is the diffusion CFL number. Then the expanded matrix is still "tridiagonal" in a 3-D sense and can be blocked into three identical parts: M_i=C-[1,0,0,0,0,0,0] on each row. The RHS vector are source terms that can be easily obtained by the following the CR scheme:

S_{xi}=(-C-[1,0,0,0,0,0,0]) \cdot B_{xi}

The matrix equation can be solved using linear solvers at hand in parallel.
The field obtained above resides at cell centers. We can do interpolations to retrieve the new magnetic field at each cell faces. But since we are treating the x, y and z fields independently and applying sources that is not necessarily divergence free to them, we are running the risk of introducing divergence to the face centered magnetic field components which should be divergence free. Experiments show that this decomposing method can introduce large divergence to the face centered field. Thus a divergence cleaning scheme is required once we finish the diffusion update. Here we introduce the divergence cleaning scheme implemented in AstroBEAR. Considering a cell with face centered field components:

We can find the divergence source term: d=\nabla \cdot B. It is obvious that if B is face centered, then its divergence should be cell centered.
Now we solve the following equation system: \nabla^2 \phi = d by using d as a source. This equation system is again well defined, and can be constructed using the same C vector described before with k = 1. Once we obtained phi, we can go back to clean the face centered field:

B'=B-\nabla \phi

We can easily see that the new field is divergence free:

\nabla \cdot B' = \nabla \cdot B - \nabla \cdot \nabla \phi = \nabla \cdot B - \nabla^2 \phi = 0

Since phi is cell centered, its gradient is face centered which matches the face centered field components. So this method is exact.

The advantages and disadvantages are listed below:

(A1) Simple scheme, works fine with slowly varying temperature situations.

(A2) The divergence cleaning scheme is simple and exact.

(D1) The scheme requires constructing cell centered field components and then reconstruct the face centered field components from the updated cell centered field. The cell centered field components are not divergence free to start with, which is diffused using a decomposed scheme. The process may introduce divergence that is large enough that the cleaning mechanism could not work well enough.

(D2) The divergence cleaning scheme has a flexible boundary condition. Normally in the code, we take the Neumann boundary condition. This condition impose the requirement that the face centered field components at the domain boundaries do not need to be adjusted. The condition in itself may not be optimistic, and is likely to be depending on the specific problem.

(D3) Cannot treat fast varying temperature situation where the temperature length scale is comparable or smaller comparing to the field length scale.

(2)Solving the coupled format with varying temperature (Explicit)

Here we consider the full equation:

\frac{\partial \textbf{B}}{\partial t}=\nabla \times (\eta \nabla \times \textbf{B})

Here we treat the face centered field components directly by looking at the following diagram:

]]

We use Bz at the top face of the cell as an example for demonstration. Here, the changing of Bz does not only depend on the Bz alone, but also the field components of Bx and By. The second order accurate scheme requires 12 closest points surrounding the face center. In the diagram, F denotes face centers, X denotes edge centers. The face centers F1 to F8 have distance of sqrt{1/2}, and the face centers F9 to F12 have distance of 1. A larger lattice cell can be considered using further face centers when a higher accuracy is required. The correlation vector C in this case has 13 elements which relates the face center we are considering F0 to the 12 adjacent face centers F1 to F12. Its elements can be written as:

C_0 = -(k_{X0} + ... + k_{X4})
C_1 = -k_{X0}
C_2 = k_{X0}
C_3 = k_{X2}
C_4 = -k_{X2}
C_5 = k_{X3}
C_6 = -k_{X3}
C_7 = -k_{X4}
C_8 = k_{X4}
C_9 = k_{X1}
C_{10} = k_{X2}
C_{11} = k_{X3}
C_{12} = k_{X4}

where k are the CFL number at the edge centers. Since the CFL are defined at the cell centers initially, we need to interpolate the edge centered CFL using the four nearest cell center temperatures. The relation between temperature and the diffusivity is described in the next section. Remember that because the x, y, z field components now reside on different face centers, the value of C vector for a certain cell will be different for different field components.

The matrix equation has the form of:

]]

where M_x, M_y and M_z are the diagonal contributions from the face centers F9 to F12: M = [C_0, C_9, C_10, C_11, C_12, ….]-[1,0,0,0…]; and the blocked M is different for x, y and z components as stated before. There are also nondiagonal blocks which connect Bx block to By and etc. These are marked by D, which comes from the C_1 to C_8 contributions. The source vector is obtained by construct the matrix with -C_0-1; C1; C2; … C12 and dot it with the field vector: [Bx1 … Bxn; By1 … Byn; Bz1 … Bzn]. The scheme involves looking at each face centers and find out it C vector. Once the C vector is obtained, it can be used to construct the matrix and the source vector.
This matrix equation is still solvable using linear solver package in parallel, but the construction of the linear system is much more complicated comparing to approach (1).

(A1) Overall more sophisticated than (1), can treat rapid changing temperature.
(A2) Method can be extended to treating viscosity.
(A3) Using a curl field as the source thus preventing unwanted field divergence from generating.
(D1) Computationally much more intense while still requires subcycling. The CFL is instead given by \eta dt/l^2

Microphysical Resistivity as Explicit Source Term

For each cell, the resistive source term is calculated by taking the curl of face centered magnetic field. This results in a 12 dimension array that stores the current source terms running over the 12 edges of a cell in 3-D. This array is then used to calculate the new face centered magnetic field. The stencil for this explicit solver is a 3 by 3 cube surrounding the cell we want to update. See the following diagram:

In cgs units, the real magnetic diffusivity can be written as:

\eta = c^2/4\pi \sigma

where \sigma is the material conductivity which depends on the electron-ion collision rate \nu_{ei} which in term depends on the electron temperature T^{-3/2}.
The cross field diffusivity can thus be expressed as function of electron temperature if we assume the plasma is fully ionized:

\eta = 8.2243 \times 10^5 Z_{eff} F(Z_{eff}) \ln \Lambda T_{e}^{-3/2}

where T is the electron temperature in eV, Zeff is the effective ion charge. The Coulomb logarithm is given by:

THe Coulomb logarithm for interested number density is plotted below:

Here, T is the same as before, n is the electron number density in the units of per cubic centimeter.

F(Z_{eff}) is a function of effective ion charge:

F(Z) = \frac{1+1.198Z+0.222Z^2}{1+2.996Z+0.753Z^2}

Since the electron's maxwellian will be distorted when there is a strong outer magnetic field applied, the realistic resistivity can depend on the local field orientation.
The anisotropicity is around 2 for the parallel and cross field resistivity. This requires the electron mean free path to be much longer than the gyroradius under such a field:
R_B << \lambda_{mfp}

or, the gyro frequency to be much greater than the mean electron ion collision frequency.
This effect is not considered currently.

Resistivity Interface

Currently in the code, the microphysical resistivity gets calculated at each cell edges by averaging the temperatures and density information surrounding that edge.

The resistive source term can be turned on in astrobear by putting the resistive flag to true:
lresistive = .true.
The second parameter ResType determines which type of resistivity one requires: 1 for constant computational resistivity; 2 for microphysical resistivity.
The third parameter resistivity is the user defined constant resistivity. When ResType = 2, the resistivity is determined by the following expression:
\eta = resistivity \times T^{-3/2}
Here, all variables are in computational units.
If ResType = 3, the resistivity is microphysical and is automatically calculated using the scaling parameters. the "resistivity" parameter will be ignored
The F(Z_{eff}) function gets calculated automatically via and is a global const.
If ResType = 3, The microphysical resistivity automatically scales with the input scaling parameter.

Thermal Conduction

The anisotropic thermal conduction is determined by the following equations:

\kappa_{align} = 5.6 \times 10^{-7} T^{5/2}
\kappa_{perp} = 3.3 \times 10^{-16} \frac{n^2}{T^{1/2}B^2}

The ratio between the cross and parallel diffusion can be written as:

r = 1.7 \times 10^{-3} \frac{n \beta}{T^4}

where T is in per 100 Kelvin, n is in per cubic centimeter.

Another user controlled isotropic thermal diffusion is added to make the code stubborn. There are . The microphysical thermal conduction is also flux limited, heat flow at each cell face
centers cannot surpass a fraction of electron mean thermal speed (Cowie, Mckee 1977):

q_{sat} = 5\phi \rho c_s^3

However, the user controller isotropic thermal conduction is not included in this limiting.
It is supposed to be several orders of magnitudes smaller than the microphysical conductivity, and should maintain isotropic.

The thermal conduction is implemented in a similar manner as the resistivity.
First, the heat flux at each cell centers are calculated: on x interfaces, qx are directly obtained, qy is obtained by averaging the adjacent 4 qy, Same for y interfaces.
We do projections of heat flux twice: first onto the field direction, then to the face normal.

We also calculate the saturation flux and project to the face normal. Remember, the saturation flux is always aligned with the field.
Both the projected heat flux and saturation flux are feed into the saturation function to calculate the final flux. The saturation starts at a "saturation point", which is controllable.
Below the saturation point, there is no saturation so there is no additional modification to the original heat flux. Above the saturation point, the flux gets saturated and gradually reach
the saturation flux.

Thermal Conduction Interface
The anisotropic thermal conduction can be turned on by putting:
lresistive = .true.
The CondType controls the type of conduction:
CondType = 1: computational conductivity, linear, only alone the field lines. CondType = 2: computational conductivity, nonlinear, only alone the field lines.
CondType = 3: microphysical conductivity, automatically calculated and scaled.
conductivity: effective only when the conductivity is computational.
lisodiffusion: a bool that controls whether there is artificial thermal diffusion. This diffusion is isotropic.
isodiffusionratio: controls the ratio between artificial diffusion and actual diffusion. For instance, setting this ratio to 0.01 will add an artificial isotropic thermal diffusion
that is 100 times weaker than the minimal actual diffusion on the grid.

===================================

Multiphysics Integration

===================================

Taking viscosity as an example. The goal for the multiphysics modules is to calculate a change of velocisty inside a grid, using the known hydro vars.
This integration can be done before any hyperbolic updates. We denote the hyperbolic update, integrated source term calculation and explicit source term calculation as H, S and E.
The operators applied are as follows:

E S^{1/2} H S^(1/2)

But this operator does not commute with its own transpose, meaning that

E S^{1/2} H S^(1/2) \neq S^{1/2} H S^(1/2) E

This is obvious with for example, viscosity, in which the E operator raises temperature and lowers the velocity profile. The hyperbolic solver will have a different response to the data modified by E comparing to that of unmodified.
A symmetric operator splitting can be achieved if we do

E^{1/2} S^{1/2} H S^(1/2)E^{1/2}

Since the commutator of this operator with its transpose will be zero.
For any given grid, the external explicit solver requires two steps of updating. The tricky part of the solver is to figure out what to use at the grid boundaries. The three different boundary situations are summarized below:
(1)For any physical boundary, use the physical boundary condition.(outflow: impose a fixed velocity; and so on)
(2)For any internal boundary to direct neighbors, get the needed data directly by communication.
(3)For any internal boundary to parent level, we need the following ghosting scheme:



When advancing dt, we have 10 ghost zones, we also need to explicitly store \dot{v}_a dt/2, the time derivative of velocity advection to be used for adjacent child grids.
The steps are:
1) apply E(dt/2) to internal + ghost 1~8 using ghost 9 and 10.
2) apply H(dt/2) to internal + ghost 1~4 using ghost 5~8.
3) apply another H(dt/2) to internal using ghost 1~4.
4) manually update ghost 1~2 using \dot{v}_a dt/2.
5) apply E(dt/2) to internal again using the updated ghost 1~2.

From here, we can see two problems of putting the explicit operator into the source term integration:
1) each complete update requires two more ghost zones on parent-child boundaries which does not fit into the source term subcyclilng. (N subcycling requires 2N ghost zones)
2) for internal boundary to direct neighbors, the two grids need to communicate in order to . This again does not fit into the source term subcycling. (there is no communications happening during subcycling in the current code)