Diffusion Solvers in AstroBEAR

Astrobear supports both implicit and explicit thermal diffusion that is both isotropic and anistropic. We first derive the equations for the implicit-anisotropic case, and then discuss simplifications that arise for isotropic and explicit versions.

Basic Equations

The anisotropic conduction equation is

\rho c_v\frac{\partial T}{\partial t} = \nabla \cdot \left [ n \hat{b} \left ( \chi_\parallel - \chi_\perp \right ) \left ( \hat{b} \cdot \nabla T \right ) + n \chi_\perp \nabla T \right ]

however, the conduction coefficients have a Temperature dependence.

\chi_\parallel = \kappa_\parallel T^{5/2} n^{-1}

\chi_\perp = \kappa_\perp \frac{n}{B^2 T^{1/2}}

Just FYI, some references will define variables a bit differently, and write these equations like this:

\rho c_v\frac{\partial T}{\partial t} = \nabla \cdot \left [\hat{b} \left ( \kappa_\parallel - \kappa_\perp \right ) \left ( \hat{b} \cdot \nabla T \right ) + \kappa_\perp \nabla T \right ]

\kappa_\parallel = c_\parallel T^{5/2}

\kappa_\perp = c_\perp \frac{n^2}{B^2 T^{1/2}}

where the c's are coefficients usually given as some number in references.

We will use \kappa_{\parallel} = c_{\parallel} T^{\lambda_{\parallel}} and the same for \kappa_{\perp} where \lambda_{\parallel}=5/2 and \lambda_{\perp}=-1/2

Collecting power's of T

So we can rewrite the equations as

\rho c_v \frac{\partial T}{\partial t} = \nabla \cdot \left [ \hat{b} \left ( \kappa_\parallel - \frac{n^2 \kappa_\perp}{B^2} \right ) \left ( \hat{b} \cdot \nabla T \right ) + \frac{n^2 \kappa_\perp}{B^2} \nabla T \right ]

Einstein simplification

or in Einstein notation

\rho c_v \partial_t T = \partial_i \left [ b_i \left (\kappa_\parallel \left ( b_j \partial_j T\right ) - \frac{n^2 \kappa_\perp}{B^2 } \left ( b_j \partial_j T \right ) \right ) + \frac{n^2 \kappa_\perp}{B^2 } \partial_i T \right ]

or

\rho c_v \partial_t T = \partial_i \kappa_\parallel b_i b_j \partial_j T + \partial_i \kappa_\perp \left ( \delta_{ij} - b_i b_j \right ) \partial_j T

Implicitization

Now we can rewrite the equation

\partial_t T = \displaystyle \sum_{\parallel, \perp} A \partial_i B_{ij} \partial_j T

where the \parallel or \perp subscript on A, B_{ij}, and \lambda is implied.

A_\parallel = \frac{c_\parallel}{\rho c_v } and A_\perp = \frac{c_\perp}{\rho c_v }

and B_{\parallel, ij} = T^{\lambda_{\parallel}}b_i b_j and B_{\perp, ij} = n^2 B^{-2} T^{\lambda_\perp} \left ( \delta_{ij}-b_i b_j \right )

Now to solve this implicitly, we need to replace T with T_* where

T_* = T + \phi ( T' - T ) = (1-\phi) T + \phi T'

Note for Backward Euler, \phi = 1 and for Crank Nicholson, \phi = 1/2

\partial_t T = \displaystyle \sum_{\parallel, \perp} A \partial_i B_{ij} \partial_j T_*

\partial_t T = \displaystyle \sum_{\parallel, \perp} A \partial_i B_{ij} \left [ C \partial_j T + D \partial_j T' \right]

where C= \left ( 1 - \phi \right ) and D = \phi

Now we can expand the derivatives and get

Discretization

We then need expressions for

\partial_t T = \frac{ T'_0 - T_0}{\Delta t}

\partial_i B_{ij} \partial_j T = \left [ \frac{B^{ij}_{\hat{i}} T_{\hat{i}+\hat{j}} - B^{ij}_{\hat{i}}T_{\hat{i}-\hat{j}}-B^{ij}_{-\hat{i}} T_{-\hat{i}+\hat{j}} + B^{ij}_{-\hat{i}}T_{-\hat{i}-\hat{j}}}{4 \Delta x^2} \right ] \left ( 1-\delta_{ij} \right ) + \left [ \frac{B^{ij}_{\hat{i/2}} T_{\hat{i}} - B^{ij}_{\hat{i/2}}T_{0}-B^{ij}_{-\hat{i/2}} T_{0} + B^{ij}_{-\hat{i/2}}T_{\hat{-i}}}{\Delta x^2} \right ] \delta_{ij}

Using the above definitions, we can write the discretized equation as

T_0 + \Delta t \displaystyle \sum_{\parallel, \perp} AC \left [ \alpha_0 T_0 + \displaystyle \sum_{\pm j} \alpha_{\pm j} T_{\pm \hat{j}} + \sum_{\pm i, \pm j,i \ne j} \alpha_{\pm i, \pm j} T_{\pm \hat{i} \pm \hat{j}} \right ] = T'_0 - \Delta t \displaystyle \sum_{\parallel, \perp} AD \left [ \alpha_0 T'_0 + \displaystyle \sum_{\pm j} \alpha_{\pm j} T'_{\pm \hat{j}} + \sum_{\pm i, \pm j,i \ne j} \alpha_{\pm i, \pm j} T'_{\pm \hat{i} \pm \hat{j}} \right ]

where

\alpha_0 = - \displaystyle \sum_{i} \frac{B^{ii}_{\hat{i}} + B^{ii}_{-\hat{i}}}{\Delta x^2}

\alpha_{\pm i} = \frac{B^{ii}_{\pm \hat{i}}}{\Delta x^2}

\alpha_{\pm i, \pm j} = \pm \pm \frac{B^{ij}_{\pm \hat{i}}}{4 \Delta x^2} where the \pm in B_{\pm \hat{i}} corresponds to the \pm i term.

Also note, that the indexing of the temperatures is commutative, and that all of the diagonal temperature terms will have two contributions. One from \alpha_{\pm i \pm j} and one from \alpha_{\pm j \pm i}. We can rewrite

\displaystyle \sum_{\pm i, \pm j,i \ne j} \alpha_{\pm i, \pm j} T^{\lambda+1}_{\pm \hat{i} \pm \hat{j}} = \sum_{\pm i, \pm j,i < j} \left ( \alpha_{\pm i, \pm j} + \alpha_{\pm j, \pm i} \right ) T^{\lambda+1}_{\pm \hat{i} \pm \hat{j}} = \sum_{\pm i, \pm j,i < j} \alpha*_{\pm i, \pm j}T^{\lambda+1}_{\pm \hat{i} \pm \hat{j}}

where

\alpha*_{\pm i \pm j} = \pm \pm \frac{B^{ij}_{\pm \hat{i}}+B^{ji}_{\pm \hat{j}}}{4 \Delta x^2}

Summary

	\parallel	\perp
A	\frac{c_\parallel}{\rho c_v }	\frac{c_\perp}{\rho c_v }
B_{ij}	b_i b_j T^{\lambda_{\parallel}}	n^2 B^{-2} T^{\lambda_{\perp}} \left (\delta_{ij} - b_i b_j \right )
C	\left ( 1 - \phi \right )	\left ( 1 - \phi \right )
D	\phi	\phi

\alpha_0	- \displaystyle \sum_{i} \frac{B^{ii}_{\hat{i}} + B^{ii}_{-\hat{i}}}{\Delta x^2}
\alpha_{\pm i}	\frac{B^{ii}_{\pm \hat{i}}}{\Delta x^2}
\alpha*_{\pm i \pm j}	\pm \pm \frac{B^{ij}_{\pm \hat{i}}+B^{ji}_{\pm \hat{j}}}{4 \Delta x^2}

Effective diffusion coeffiient ==

The traditional diffusion equation looks like

\frac{\partial T}{\partial t} = D \nabla^2 T

where D is the diffusion coefficient. For uniform density, magnetic field, and for \lambda = 0, the equations revert to the traditional diffusion equation with

D = \frac{n \chi}{\rho c_v}

Isotropic Case

For the Isotropic Case, the diffusion equation looks like

\rho c_v \frac{\partial T}{\partial t} = \nabla \cdot \kappa \nabla T

or in Einstein notation

\rho c_v \partial_t T = \partial_i \kappa \partial_i T

which we can rewrite as

\partial_t T = A \partial_i B \partial i T

\partial_t T = A \partial_i B \delta_{ij} \partial j T

This is the exact same form as before, but now

B_{ij}=\kappa \delta{ij}

is diagonal so the cross terms vanish

and \alpha_{\pm i, \pm j} = 0

and it is independent of direction, so it is just a scalar.

So in summary we have

A	\frac{1}{\rho c_v }
B_{ij}	\kappa \delta{ij}
C	\left ( 1 - \phi \right )
D	\phi

and

\alpha_0	- \displaystyle \sum_{i} \frac{B_{\hat{i}} + B_{-\hat{i}}}{\Delta x^2}
\alpha_{\pm i}	\frac{\kappa_{\pm \hat{i}}}{\Delta x^2}
\alpha*_{\pm i \pm j}	0

And our equation becomes

T_0 + \Delta t \displaystyle \sum_{\parallel, \perp} AC \left [ \alpha_0 T_0 + \displaystyle \sum_{\pm j} \alpha_{\pm j} T_{\pm \hat{j}} \right ] = T'_0 - \Delta t \displaystyle \sum_{\parallel, \perp} AD \left [ \alpha_0 T'_0 + \displaystyle \sum_{\pm j} \alpha_{\pm j} T'_{\pm \hat{j}} \right ]

Now in the code, \kappa is assumed constant

\mbox{A} = \kappa A

\mbox{C} = \frac{\Delta t A C}{\Delta x^2}

\mbox{D} = \frac{\Delta t A D}{\Delta x^2}

\mbox{Tlambda} = T_{\pm \hat{j}}^\lambda

\mbox{T} = T_{\pm \hat{j}}

\mbox{T0} = T

\mbox{T0lambda} = T^{\lambda}

\mbox{alpha} = 1

And stencil(0) has the coefficient of T_0', and stencil(1:2*ndim) has the coefficients of T_{\pm \hat{j}}' and source has everything else. Note that the coefficients \alpha_0 has a term for +\hat{i} and -\hat{i} but the way the code is written, it loops over each neighbor and adds the contribution from -\hat{i} and +\hat{i} on separate iterations.

Time stepping

The implicit method requires approximating

T*^{\lambda+1} = T^{\lambda + 1} + \left ( \lambda + 1 \right ) T^{\lambda} \left [ \left ( T_{*} - T \right ) + \frac{\lambda}{2} \frac{T_{*} - T}{T} \right ]

but solving this accurately with a linear method requires

\frac{T_{*}-T}{T} << 1

but this restricts our time step. Using the instantaneous derivatives, we can calculate the maximum time step as

T_0'-T_0 \approx \displaystyle \sum_{\parallel, \perp} \left [ \alpha_0 T^{\lambda+1}_0 + \displaystyle \sum_{\pm j} \alpha_{\pm j} T^{\lambda+1}_{\pm \hat{j}} + \sum_{\pm i, \pm j,i \ne j} \alpha_{\pm i, \pm j} T^{\lambda+1}_{\pm \hat{i} \pm \hat{j}} \right ]

so

\Delta t < \epsilon \frac{T_0}{\displaystyle \sum_{\parallel, \perp} \left [ \alpha_0 T^{\lambda+1}_0 + \displaystyle \sum_{\pm j} \alpha_{\pm j} T^{\lambda+1}_{\pm \hat{j}} + \sum_{\pm i, \pm j,i \ne j} \alpha_{\pm i, \pm j} T^{\lambda+1}_{\pm \hat{i} \pm \hat{j}} \right ]}

Explicit Case

For the Explicit case, we just set \phi = 0 which sets D = 0 and C=1 and we have

T_0 + \Delta t \displaystyle \sum_{\parallel, \perp} A \left [ \alpha_0 T^{\lambda+1}_0 + \displaystyle \sum_{\pm j} \alpha_{\pm j} T^{\lambda+1}_{\pm \hat{j}} + \sum_{\pm i, \pm j,i \ne j} \alpha_{\pm i, \pm j} T^{\lambda+1}_{\pm \hat{i} \pm \hat{j}} \right ] = T'_0

Explicit time stepping

For the explicit scheme to be stable, we must have

\Delta t < \frac{1}{2\displaystyle \max_{\parallel, \perp, i} A \left [ | \alpha_i T_i^{\lambda} | \right ]}

Flux limiting

Limiting the fluxes in the explicit case is not too difficult, however in the implicit case, we cannot limit the flux directly, but must instead, limit the conductivities.

n \chi_\parallel \hat{b} \cdot \nabla T < 5 \phi \rho c_s^3 \rightarrow \chi_\parallel < \frac{5 \phi \rho c_s^3}{n \hat{b} \cdot \nabla T} \equiv \chi_{\max}

We can limit \chi using \chi = \min{\left ( \chi, \chi_{\max} \right )} or do something like

\chi = \chi_{\max} \tanh \frac{\chi}{\chi_{\max}}

Boundary Terms

Explicit

Here, the boundary terms are quite straightforward. We can just use the regular boundary conditions (or user specified boundary conditions). The only difficulty arises when trying to set a heat flux since this involves solving a non-linear equation. One solution is to set the heat flux explicitly - instead of indirectly through the temprature. Then use extrapolating (pr reflecting) BC's for temperature so there is no heat flux initially calculated along the boundary - and then add in the flux in a separate step.

Implicit

For reflecting - or extrapolating boundary terms, we can just transfer the coefficients for the boundary terms and apply them to their mirror zone. The only thing to be careful about - is that for the anisotropic case, there are corner zones who's coefficients may need to be transferred twice. This just requires transferring corner zone coefficients, to edge zones, before transferring edge zones, to the central cell in the stencil.