\documentclass{article}
\usepackage{amssymb,amsmath}
\usepackage{verbatim}
\usepackage{graphicx}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage[square, comma, numbers]{natbib}
\usepackage{kpfonts}
\usepackage{tikz}
\usepackage{listings}
\usetikzlibrary{positioning,shapes}

\newcommand*{\vv}[1]{\vec{\mkern0mu#1}}
\newcommand{\myparallel}{{\mkern3mu\vphantom{\perp}\vrule depth 0pt\mkern2mu\vrule depth 0pt\mkern3mu}}

\lstset{language=[90]Fortran,
   basicstyle=\ttfamily,
   keywordstyle=\color{red},
   commentstyle=\color{green},
 %  morecomment=[1]{!\ } % comment only with space after !
}
\begin{document}
\title{Three-Dimensional Anisotropic Thermal Conduction Solver in AstroBEAR 3.0}
\author{Baowei Liu}
\maketitle

\section{Anisotropic Thermal Conduction Equation}
\subsection{Equations}
 Thermal conduction in plasmas with the mean free path much larger than the gyroradius (also know as cyclotron radius, is the radius of the circular motion of a charged particle in the presence of a uniform magnetic field) is anisotropic with respect to the magnetic field line; heat flow primarily along the filed lines with little conduction in the perpendicular direction. In such cases, a divergence of anisotropic heat flux is added to the energy equation. Thermal conduction can modify the characteristic structure of the magnetohydrodynamic equations making it difficult to incorporate into upwind methods. However, thermal conduction can be evolved independently of the MHD equations using operator splitting, as done in [10]. The equation for the evolution of internal energy density due to anisotropic thermal conduction is 
\begin{equation}
\vec{q} = \hat{b}n(\chi_{\myparallel}-\chi_{\perp})(\hat{b}\cdot\nabla)T+n\chi_{\perp}\nabla T
\label{eq:flux}
\end{equation}
where $\hat{b}$ is the unit vector along the field direction, $n$ is the number density of electrons/ions?, $\chi$s are the temperature dependent thermal conductivities. The detailed expressions are
\begin{equation}
\chi_{\myparallel}=\kappa_{\myparallel}T^{5/2}
\label{eq:chi_pr}
\end{equation}
\begin{equation}
\chi_{\perp}=\kappa_{\perp}\frac{n}{B^{2}T^{1/2}}
\label{eq:chi_pd}
\end{equation}
%\begin{equation}
%\begin{cases} \chi_{\myparallel}=\kappa_{\myparallel}T^{5/2}\\ \chi_{\perp}=\kappa_{\perp}\frac{n}{B^2T^{1/2}}
%\end{cases}
%\label{eq:twoeq}
%\end{equation}

%\begin{equation}
%\begin{cases} \chi_{\myparallel}=\kappa_{\myparallel}T^{5/2}\\ \chi_{\perp}=0
%\end{cases}
%\end{equation}
%
%and the heat diffusion equation

%For 2D we should have:
%\begin{equation}
%\begin{cases}
%q_{x}=-b_{x}n(\chi_{\myparallel}-\chi_{\perp})(b_{x}\frac{\partial T}{\partial x}+b_{y}\frac{\partial T}{\partial y})-n\chi_{\perp}\frac{\partial T}{\partial y}\\
%q_{y}=-b_{y}n(\chi_{\myparallel}-\chi_{\perp})(b_{x}\frac{\partial T}{\partial x}+b_{y}\frac{\partial T}{\partial y})-n\chi_{\perp}\frac{\partial T}{\partial x}
%\end{cases}
%\end{equation}

The heat diffusion equation in AstroBEAR
\begin{equation}
\frac{1}{\gamma -1}\rho C_{v}\frac{K_{B}\partial T}{\partial t}=\nabla\cdot\frac{1}{\gamma -1}\hat{q}
\label{eq:diffEq}
\end{equation}

For simplicity we write
\begin{equation}
\frac{\partial T}{\partial t}=\nabla\cdot\hat{q}
\label{eq:HD_case1}
\end{equation}

\subsection{Simpler Case Without Perpendicular Conductivity}
We first consider the case:

\begin{equation}
\begin{cases} \chi_{\myparallel}=\kappa_{\myparallel}T^{5/2}\\ \chi_{\perp}=0
\end{cases}
\end{equation}

Here the parallel conductivity is much larger than the perpendicular conductivity (larger by a factor of $10^9$). These expressions are from the Orlando simulation. In a different scenario (i.e. when considering some dense and high pressure plasma), the expressions may be different. The heat flux is

\begin{equation}
\vec{q} = \hat{b}n\chi_{\myparallel}(\hat{b}\cdot\nabla T)
\label{eq:appoxflux}
\end{equation}
So
\begin{equation}
\begin{aligned}
\nabla\cdot\vec{q}&=\nabla\cdot\left[\hat{b}n\chi_{\myparallel}(\hat{b}\cdot\nabla T) \right]\\
&=n\chi_{\myparallel}(\hat{b}\cdot\nabla)^{2}T+n(\hat{b}\cdot\nabla T)(\hat{b}\cdot\nabla\chi_{\myparallel})
\end{aligned}
\label{eq:fluxdev_case1}
\end{equation}

The first part in \ref{eq:fluxdev_case1} is
\begin{equation}
\begin{aligned}
n\chi_{\myparallel}(\hat{b}\cdot\nabla)^{2}T= &n\chi_{\myparallel}\left[\left(b_{x}^{2}\frac{\partial^{2} T}{\partial x^{2}}+b_{x}b_{y}\frac{\partial^{2} T}{\partial x\partial y}+b_{x}b_{z}\frac{\partial^{2} T}{\partial x\partial z}\right)\right. + \\
&\left(b_{y}b_{x}\frac{\partial^{2} T}{\partial x\partial y}+b^{2}_{y}\frac{\partial^{2} T}{\partial y^{2}}+b_{y}b_{z}\frac{\partial^{2} T}{\partial y\partial z}\right) + \\
&\left.\left(b_{z}b_{x}\frac{\partial^{2} T}{\partial z\partial x}+b_{z}b_{y}\frac{\partial^{2} T}{\partial z\partial y}+b^{2}_{z}\frac{\partial^{2} T}{\partial z^{2}}\right)\right]
\end{aligned}
\label{eq:1stp_fluxdev_case1}
\end{equation}

And since from Eq. \ref{eq:chi_pr}
\begin{equation}
\hat{b}\cdot\nabla\chi_{\myparallel}=2.5\kappa_{\myparallel}T^{3/2}(\hat{b}\cdot\nabla T)
\end{equation}

So the second part in Eq. \ref{eq:fluxdev_case1} is
\begin{equation}
\begin{aligned}
&n(\hat{b}\cdot\nabla T)(\hat{b}\cdot\nabla\chi_{\myparallel})\\
&=2.5n\kappa_{\myparallel}T^{3/2}(\hat{b}\cdot\nabla T)^2\\
&=2.5n\kappa_{\myparallel}T^{3/2}\left[b^{2}_{x}(\frac{\partial T}{\partial x})^{2}+b^{2}_{y}(\frac{\partial T}{\partial y})^{2} + b^{2}_{z}(\frac{\partial T}{\partial z})^{2}\right. \\
&+\left. 2b_{x}b_{y}\frac{\partial T}{\partial x}\frac{\partial T}{\partial y}+2b_{x}b_{z}\frac{\partial T}{\partial x}\frac{\partial T}{\partial z}+2b_{y}bz{y}\frac{\partial T}{\partial y}\frac{\partial T}{\partial z}   \right]
\end{aligned}
\label{eq:2ndp_fluxdev_case1}
\end{equation}

So Eq. \ref{eq:HD_case1} becomes

%\begin{equation}
%\begin{aligned}
%\frac{\partial T}{\partial t}=&2.5\kappa_{\myparallel}T^{3/2}n\left[b^{2}_{x}\left(\frac{\partial T}{\partial x}\right)^{2}+b^{2}_{y}\left(\frac{\partial T}{\partial y}\right)^{2}+b^{2}_{z}\left(\frac{\partial T}{\partial z}\right)^{2}+\right.\\
%&\left.2b_{x}b_{y}\frac{\partial T}{\partial x}\frac{\partial T}{\partial y} + 2b_{x}b_{z}\frac{\partial T}{\partial x}\frac{\partial T}{\partial z} +2b_{y}b_{z}\frac{\partial T}{\partial y}\frac{\partial T}{\partial z}\right]+\\
%&n\kappa_{\myparallel}T^{2.5}\left[\left(b_{x}^{2}\frac{\partial^{2} T}{\partial x^{2}}+b_{x}b_{y}\frac{\partial^{2} T}{\partial x\partial y}+b_{x}b_{z}\frac{\partial^{2} T}{\partial x\partial z}\right)\right. + \\
%&\left(b_{y}b_{x}\frac{\partial^{2} T}{\partial y\partial x}+b^{2}_{y}\frac{\partial^{2} T}{\partial y^{2}}+b_{y}b_{z}\frac{\partial^{2} T}{\partial y\partial z}\right) + \\
%&\left.\left(b_{z}b_{x}\frac{\partial^{2} T}{\partial z\partial x}+b_{z}b_{y}\frac{\partial^{2} T}{\partial z\partial y}+b^{2}_{z}\frac{\partial^{2} T}{\partial z^{2}}\right)\right]
%\end{aligned}
%\label{eq:HD_case1_3}
%\end{equation}
\begin{equation}
\begin{aligned}
\frac{\partial T}{\partial t}=&2.5\kappa_{\myparallel}T^{3/2}n\left[b^{2}_{x}\left(\frac{\partial T}{\partial x}\right)^{2}+b^{2}_{y}\left(\frac{\partial T}{\partial y}\right)^{2}+b^{2}_{z}\left(\frac{\partial T}{\partial z}\right)^{2}+\right.\\
&\left.2b_{x}b_{y}\frac{\partial T}{\partial x}\frac{\partial T}{\partial y} + 2b_{x}b_{z}\frac{\partial T}{\partial x}\frac{\partial T}{\partial z} +2b_{y}b_{z}\frac{\partial T}{\partial y}\frac{\partial T}{\partial z}\right]+\\
&n\kappa_{\myparallel}T^{2.5}\left[\left(b_{x}^{2}\frac{\partial^{2} T}{\partial x^{2}}+2b_{x}b_{y}\frac{\partial^{2} T}{\partial x\partial y}+2b_{x}b_{z}\frac{\partial^{2} T}{\partial x\partial z}\right)\right. + \\
&\left(b^{2}_{y}\frac{\partial^{2} T}{\partial y^{2}}+2b_{y}b_{z}\frac{\partial^{2} T}{\partial y\partial z}\right) + \\
&\left.\left(b^{2}_{z}\frac{\partial^{2} T}{\partial z^{2}}\right)\right]
\end{aligned}
\label{eq:HD_case1_3}
\end{equation}

\subsection{With Perpendicular Conductivity}
Then we consider the case:

\begin{equation}
\chi_{\myparallel}=\kappa_{\myparallel}T^{5/2}
\label{eq:chi_pr2}
\end{equation}
\begin{equation}
\chi_{\perp}=\kappa_{\perp}\frac{n}{B^{2}T^{1/2}}
\label{eq:chi_pd2}
\end{equation}

From Eq. \ref{eq:flux} we have
\begin{equation}
\begin{aligned}
\nabla\cdot \vec{q}&=\nabla\cdot\left[\hat{b}n(\chi_{\myparallel}-\chi_{\perp})\left(\hat{b}\cdot\nabla T\right) \right]+\chi_{\perp}n\nabla^{2}T \\
 &=n(\hat{b}\cdot\nabla T)\left[\hat{b}\cdot\nabla(\chi_{\myparallel}-\chi_{\perp})\right]+n(\chi_{\myparallel}-\chi_{\perp})(\hat{b}\cdot\nabla)^{2}T+n\chi_{\perp}\nabla^{2}T 
\end{aligned}
\label{eq:pd_diff1}
\end{equation}

Since 
\begin{equation}
\hat{b}\cdot\nabla(\chi_{\myparallel}-\chi_{\perp})=\left[2.5\kappa_{\myparallel}T^{3/2}+0.5\kappa_{\perp}\frac{n}{B^{2}T^{3/2}} \right](\hat{b}\cdot\nabla T)
\end{equation}

so similar as \ref{eq:HD_case1_3} the first term on the right of Eq.\ref{eq:pd_diff1} is

\begin{equation}
\begin{aligned}
&n(\hat{b}\cdot\nabla T)\left[\hat{b}\cdot\nabla(\chi_{\myparallel}-\chi_{\perp})\right]\\
&=n\left(2.5\kappa_{\myparallel}T^{3/2}+0.5\kappa_{\perp}\frac{n}{B^{2}T^{3/2}} \right)(\hat{b}\cdot\nabla T)^{2}\\
&=n\left(2.5\kappa_{\myparallel}T^{3/2}+0.5\kappa_{\perp}\frac{n}{B^{2}T^{3/2}} \right)\left[b^{2}_{x}\left(\frac{\partial T}{\partial x}\right)^{2}+b^{2}_{y}\left(\frac{\partial T}{\partial y}\right)^{2}+b^{2}_{z}\left(\frac{\partial T}{\partial z}\right)^{2}+\right.\\
&\left.2b_{x}b_{y}\frac{\partial T}{\partial x}\frac{\partial T}{\partial y} + 2b_{x}b_{z}\frac{\partial T}{\partial x}\frac{\partial T}{\partial z} +2b_{y}b_{z}\frac{\partial T}{\partial y}\frac{\partial T}{\partial z}\right]
\end{aligned}
\label{eq:1stp_pddiff}
\end{equation}

and the second term on the right of Eq.\ref{eq:pd_diff1} is
\begin{equation}
\begin{aligned}
&n(\chi_{\myparallel}-\chi_{\perp})(\hat{b}\cdot\nabla)^{2}T\\
&=n\left(\kappa_{\myparallel}T^{5/2}-\kappa_{\perp}\frac{n}{B^{2}T^{1/2}} \right)(\hat{b}\cdot\nabla)^{2}T\\
&=n\left(\kappa_{\myparallel}T^{5/2}-\kappa_{\perp}\frac{n}{B^{2}T^{1/2}} \right)\left[\left(b_{x}^{2}\frac{\partial^{2} T}{\partial x^{2}}+b_{x}b_{y}\frac{\partial^{2} T}{\partial x\partial y}+b_{x}b_{z}\frac{\partial^{2} T}{\partial x\partial z}\right)\right. + \\
&\left(b_{y}b_{x}\frac{\partial^{2} T}{\partial y\partial x}+b^{2}_{y}\frac{\partial^{2} T}{\partial y^{2}}+b_{y}b_{z}\frac{\partial^{2} T}{\partial y\partial z}\right) + \\
&\left.\left(b_{z}b_{x}\frac{\partial^{2} T}{\partial z\partial x}+b_{z}b_{y}\frac{\partial^{2} T}{\partial z\partial y}+b^{2}_{z}\frac{\partial^{2} T}{\partial z^{2}}\right)\right]
\end{aligned}
\label{eq:2ndp_pddiff}
\end{equation}

and finally the third term on the right of Eq.\ref{eq:pd_diff1} is

\begin{equation}
n\chi_{\perp}\nabla^{2}T = n\kappa_{\perp}\frac{n}{B^{2}T^{1/2}}\left(\frac{\partial^{2} T}{\partial x^{2}}+\frac{\partial^{2} T}{\partial y^{2}}+\frac{\partial^{2} T}{\partial z^{2}}\right)
\label{eq:3rdp_pddiff}
\end{equation}

So Eq. \ref{eq:HD_case1} becomes
\begin{equation}
\begin{aligned}
\frac{\partial T}{\partial t}=&n\left(2.5\kappa_{\myparallel}T^{3/2}+0.5\kappa_{\perp}\frac{n}{B^{2}T^{3/2}} \right)\left[b^{2}_{x}\left(\frac{\partial T}{\partial x}\right)^{2}+b^{2}_{y}\left(\frac{\partial T}{\partial y}\right)^{2}+b^{2}_{z}\left(\frac{\partial T}{\partial z}\right)^{2}+\right.\\
&\left.2b_{x}b_{y}\frac{\partial T}{\partial x}\frac{\partial T}{\partial y} + 2b_{x}b_{z}\frac{\partial T}{\partial x}\frac{\partial T}{\partial z} +2b_{y}b_{z}\frac{\partial T}{\partial y}\frac{\partial T}{\partial z}\right]+\\
&n\left(\kappa_{\myparallel}T^{5/2}-\kappa_{\perp}\frac{n}{B^{2}T^{1/2}} \right)\left[\left(b_{x}^{2}\frac{\partial^{2} T}{\partial x^{2}}+2b_{x}b_{y}\frac{\partial^{2} T}{\partial x\partial y}+2b_{x}b_{z}\frac{\partial^{2} T}{\partial x\partial z}\right)\right. + \\
&+\left(b^{2}_{y}\frac{\partial^{2} T}{\partial y^{2}}+2b_{y}b_{z}\frac{\partial^{2} T}{\partial y\partial z}\right) + \\
&\left.\left(b^{2}_{z}\frac{\partial^{2} T}{\partial z^{2}}\right)\right]+\\
&n\kappa_{\perp}\frac{n}{B^{2}T^{1/2}}\left(\frac{\partial^{2} T}{\partial x^{2}}+\frac{\partial^{2} T}{\partial y^{2}}+\frac{\partial^{2} T}{\partial z^{2}}\right)
\end{aligned}
\label{eq:anisohd}
\end{equation}

%\begin{equation}
%\begin{aligned}
%\frac{\partial T}{\partial t}=&n\left(2.5\kappa_{\myparallel}T^{3/2}+0.5\kappa_{\perp}\frac{n}{B^{2}T^{3/2}} \right)\left[b^{2}_{x}\left(\frac{\partial T}{\partial x}\right)^{2}+b^{2}_{y}\left(\frac{\partial T}{\partial y}\right)^{2}+b^{2}_{z}\left(\frac{\partial T}{\partial z}\right)^{2}+\right.\\
%&\left.2b_{x}b_{y}\frac{\partial T}{\partial x}\frac{\partial T}{\partial y} + 2b_{x}b_{z}\frac{\partial T}{\partial x}\frac{\partial T}{\partial z} +2b_{y}b_{z}\frac{\partial T}{\partial y}\frac{\partial T}{\partial z}\right]+\\
%&n\left(\kappa_{\myparallel}T^{5/2}-\kappa_{\perp}\frac{n}{B^{2}T^{1/2}} \right)\left[\left(b_{x}^{2}\frac{\partial^{2} T}{\partial x^{2}}+b_{x}b_{y}\frac{\partial^{2} T}{\partial x\partial y}+b_{x}b_{z}\frac{\partial^{2} T}{\partial x\partial z}\right)\right. + \\
%&\left(b_{y}b_{x}\frac{\partial^{2} T}{\partial y\partial x}+b^{2}_{y}\frac{\partial^{2} T}{\partial y^{2}}+b_{y}b_{z}\frac{\partial^{2} T}{\partial y\partial z}\right) + \\
%&\left.\left(b_{z}b_{x}\frac{\partial^{2} T}{\partial z\partial x}+b_{z}b_{y}\frac{\partial^{2} T}{\partial z\partial y}+b^{2}_{z}\frac{\partial^{2} T}{\partial z^{2}}\right)\right]+\\
%&n\kappa_{\perp}\frac{n}{B^{2}T^{1/2}}\left(\frac{\partial^{2} T}{\partial x^{2}}+\frac{\partial^{2} T}{\partial y^{2}}+\frac{\partial^{2} T}{\partial z^{2}}\right)
%\end{aligned}
%\label{eq:anisohd}
%\end{equation}

The following Table gives the coefficients for each derivative term
\begin{equation}
\begin{aligned} 
&\frac{\partial^{2} T}{\partial x^{2}} & n\left[b^{2}_{x}\kappa_{\myparallel}T^{5/2}+(1-b^{2}_{x})\kappa_{\perp}\frac{n}{B^{2}T^{1/2}} \right] & = A_{11}\\
&\frac{\partial^{2} T}{\partial y^{2}} & n\left[b^{2}_{y}\kappa_{\myparallel}T^{5/2}+(1-b^{2}_{y})\kappa_{\perp}\frac{n}{B^{2}T^{1/2}} \right] & = A_{22}\\
&\frac{\partial^{2} T}{\partial z^{2}} & n\left[b^{2}_{z}\kappa_{\myparallel}T^{5/2}+(1-b^{2}_{z})\kappa_{\perp}\frac{n}{B^{2}T^{1/2}} \right] & = A_{33}\\
&\frac{\partial^{2} T}{\partial x\partial y} & 2nb_{x}b_{y}\left(\kappa_{\myparallel}T^{5/2}-\kappa_{\perp}\frac{n}{B^{2}T^{1/2}} \right) & = A_{12}\\
%&\frac{\partial^{2} T}{\partial y\partial x} & nb_{y}b_{x}\left(\kappa_{\myparallel}T^{5/2}-\kappa_{\perp}\frac{n}{B^{2}T^{1/2}} \right) & A_{21}\\
&\frac{\partial^{2} T}{\partial x\partial z} & 2nb_{x}b_{z}\left(\kappa_{\myparallel}T^{5/2}-\kappa_{\perp}\frac{n}{B^{2}T^{1/2}} \right) & = A_{13}\\
%&\frac{\partial^{2} T}{\partial z\partial x} & nb_{z}b_{x}\left(\kappa_{\myparallel}T^{5/2}-\kappa_{\perp}\frac{n}{B^{2}T^{1/2}} \right) & A_{31}\\
&\frac{\partial^{2} T}{\partial y\partial z} & 2nb_{y}b_{z}\left(\kappa_{\myparallel}T^{5/2}-\kappa_{\perp}\frac{n}{B^{2}T^{1/2}} \right) & = A_{23}\\
%&\frac{\partial^{2} T}{\partial z\partial y} & nb_{z}b_{y}\left(\kappa_{\myparallel}T^{5/2}-\kappa_{\perp}\frac{n}{B^{2}T^{1/2}} \right) & A_{32}\\
&\left(\frac{\partial T}{\partial x} \right)^2  & nb^{2}_{x}\left(2.5\kappa_{\myparallel}T^{3/2}+0.5\kappa_{\perp}\frac{n}{B^{2}T^{3/2}}  \right) & = B_{11} \\ 
&\left(\frac{\partial T}{\partial y} \right)^2  & nb^{2}_{y}\left(2.5\kappa_{\myparallel}T^{3/2}+0.5\kappa_{\perp}\frac{n}{B^{2}T^{3/2}}  \right) & = B_{22}\\
&\left(\frac{\partial T}{\partial z} \right)^2 & nb^{2}_{z}\left(2.5\kappa_{\myparallel}T^{3/2}+0.5\kappa_{\perp}\frac{n}{B^{2}T^{3/2}}  \right) & = B_{33}\\ 
&\frac{\partial T}{\partial x}\frac{\partial T}{\partial y}   & nb_{x}b_{y}\left(5\kappa_{\myparallel}T^{3/2}+\kappa_{\perp}\frac{n}{B^{2}T^{3/2}} \right)  & =  B_{12} \\ 
&\frac{\partial T}{\partial x}\frac{\partial T}{\partial z}   & nb_{x}b_{z}\left(5\kappa_{\myparallel}T^{3/2}+\kappa_{\perp}\frac{n}{B^{2}T^{3/2}} \right) &  =  B_{13} \\ 
&\frac{\partial T}{\partial y}\frac{\partial T}{\partial z}   & nb_{y}b_{z}\left(5\kappa_{\myparallel}T^{3/2}+\kappa_{\perp}\frac{n}{B^{2}T^{3/2}} \right)  & =  B_{23}\\ 
\end{aligned}
\label{eq:3dco}
\end{equation}

With the coefficients above, Eq. \ref{eq:anisohd} can be re-written as
\begin{equation}
\frac{\partial T}{\partial t}=\sum_{i,j=1,2,3}\left[A_{ij}\frac{\partial^{2} T}{\partial x_{i}\partial x_{j}}+B_{i,j}\frac{\partial T}{\partial x_{i}}\frac{\partial T}{\partial x_{j}}\right]
\end{equation}

\subsection{Discretization and Linearization}

Assume $\Delta x=\Delta y=\Delta z=a$. Discretizing Eq.\ref{eq:anisohd} using the coefficients in Table. \ref{eq:3dco} we have

%%========== equal mixed partial equation =======
%\begin{equation}
%\begin{aligned}
%\frac{T^{l+1}_{i,j,k}-T^{l}_{i,j,k}}{\Delta t}=&[n(\chi_{\myparallel}-\chi_{\perp})b_{x}^{2}+n\chi_{\perp}]\frac{T^{l}_{i+1,j,k}+T^{l}_{i-1,j,k}-2T^{l}_{i,j,k})}{a^{2}}+\\
%  &[n(\chi_{\myparallel}-\chi_{\perp})b_{y}^{2}+n\chi_{\perp}]\frac{T^{l}_{i,j+1,k}+T^{l}_{i,j-1,k}-2T^{l}_{i,j,k})}{a^{2}}+\\
%  &[n(\chi_{\myparallel}-\chi_{\perp})b_{z}^{2}+n\chi_{\perp}]\frac{T^{l}_{i,j,k+1}+T^{l}_{i,j,k-1}-2T^{l}_{i,j,k})}{a^{2}}+\\
% &2n(\chi_{\myparallel}-\chi_{\perp})b_{x}b_{y}\frac{T^{l}_{i+1,j,k}+T^{l}(i,j-1,k)-2T^{l}_{i,j,k}}{a^{2}}+\\
% &2n(\chi_{\myparallel}-\chi_{\perp})b_{x}b_{z}\frac{T^{l}_{i+1,j,k}+T^{l}(i,j,k-1)-2T^{l}_{i,j,k}}{a^{2}}+\\
% &2n(\chi_{\myparallel}-\chi_{\perp})b_{y}b_{z}\frac{T^{l}_{i,j+1,k}+T^{l}(i,j,k-1)-2T^{l}_{i,j,k}}{a^{2}}
%\end{aligned}
%\end{equation}

\begin{equation}
\begin{aligned}
\frac{T^{l+1}_{i,j,k}-T^{l}_{i,j,k}}{\Delta t}=&A_{11}\frac{T^{l}_{i+1,j,k}+T^{l}_{i-1,j,k}-2T^{l}_{i,j,k}}{a^{2}}+\\
  &A_{22}\frac{T^{l}_{i,j+1,k}+T^{l}_{i,j-1,k}-2T^{l}_{i,j,k}}{a^{2}}+\\
  &A_{33}\frac{T^{l}_{i,j,k+1}+T^{l}_{i,j,k-1}-2T^{l}_{i,j,k}}{a^{2}}+\\
 &A_{12}\frac{T^{l}_{i+1,j+1,k}-T^{l}_{i-1,j+1,k}-T^{l}_{i+1,j-1,k}+T_{i-1,j-1,k}}{2a^{2}}+\\
 &A_{23}\frac{T^{l}_{i,j+1,k+1}-T^{l}_{i,j+1,k-1}-T^{l}_{i,j-1,k+1}+T_{i,j-1,k-1}}{2a^{2}}+\\
 &A_{13}\frac{T^{l}_{i+1,j,k+1}-T^{l}_{i+1,j,k-1}-T^{l}_{i-1,j,k+1}+T_{i-1,j,k-1}}{2a^{2}}+\\
 &B_{11}\frac{T^{l}_{i+1,j,k}T^{l}_{i+1,j,k}-2T^{l}_{i+1,j,k}T^{l}_{i-1,j,k}+T^{l}_{i-1,j,k}T^{l}_{i-1,j,k}}{4a^{2}} +\\
 &B_{22}\frac{T^{l}_{i,j+1,k}T^{l}_{i,j+1,k}-2T^{l}_{i,j+1,k}T^{l}_{i,j-1,k}+T^{l}_{i,j-1,k}T^{l}_{i,j-1,k}}{4a^{2}} +\\
 &B_{33}\frac{T^{l}_{i,j,k+1}T^{l}_{i,j,k+1}-2T^{l}_{i,j,k+1}T^{l}_{i,j,k-1}+T^{l}_{i,j,k-1}T^{l}_{i,j,k-1}}{4a^{2}} +\\
 &B_{12}\frac{T^{l}_{i+1,j,k}T^{l}_{i,j+1,k}-T^{l}_{i+1,j,k}T^{l}_{i,j-1,k}-T^{l}_{i-1,j,k}T^{l}_{i,j+1,k}+T^{l}_{i-1,j,k}T^{l}_{i,j-1,k}}{4a^{2}} +\\
 &B_{13}\frac{T^{l}_{i+1,j,k}T^{l}_{i,j,k+1}-T^{l}_{i+1,j,k}T^{l}_{i,j,k-1}-T^{l}_{i-1,j,k}T^{l}_{i,j,k+1}+T^{l}_{i-1,j,k}T^{l}_{i,j,k-1}}{4a^{2}} +\\
 &B_{23}\frac{T^{l}_{i,j+1,k}T^{l}_{i,j,k+1}-T^{l}_{i,j+1,k}T^{l}_{i,j,k-1}-T^{l}_{i,j-1,k}T^{l}_{i,j,k+1}+T^{l}_{i,j-1,k}T^{l}_{i,j,k-1}}{4a^{2}} 
\end{aligned}
\label{eq:dischd1}
\end{equation}

Or in a different form

\begin{equation}
\begin{aligned}
&\frac{T^{l+1}_{i,j,k}-T^{l}_{i,j,k}}{\Delta t}\\
&=\sum_{\mu=1,2,3}A_{\mu\mu}\frac{T^{l}_{i+\delta_{\mu,1},j+\delta_{\mu,2},k+\delta_{\mu,3}}+T^{l}_{i-\delta_{\mu,1},j-\delta_{\mu,2},k-\delta_{\mu,3}}-2T^{l}_{i,j,k}}{a^{2}}\\
    & +\sum_{\mu,\nu=1,2,3}^{\mu <\nu}\frac{T^{l}_{i+\delta_{\mu,1},j+\delta_{\mu,2}+\delta_{\nu,2},k+\delta_{\nu,3}}-T^{l}_{i+\delta_{\mu,1},j+\delta_{\mu,2}-\delta_{\nu,2},k-\delta_{\nu,3}}-T^{l}_{i-\delta_{\mu,1},j-\delta_{\mu,2}+\delta_{\nu,2},k+\delta_{\nu,3}}+T^{l}_{i-\delta_{\mu,1},j-\delta_{\mu,2}-\delta_{\nu,2},k-\delta_{\nu,3}}}{a^{2}}\\
    &+
\end{aligned}
\end{equation}
To linearize the non-linear terms in Eq.\ref{eq:dischd1} we use Taylor expansion of the term at time $l+1$. For example, for the terms of $B_{11}$ we have
\begin{equation}
\begin{aligned}
\left(T^{l+1}_{i\pm 1,j,k}\right)^{2}&=(T^{l}_{i\pm 1, j,k}+2T^{l}_{i,j,k}(T^{l+1}_{i\pm 1,j,k}-T^{l}_{i\pm 1,j,k}) \\
&=2T^{l}_{i\pm 1, j,k}T^{l+1}_{i\pm 1, j,k}-(T^{l}_{i\pm 1,j,k})^{2}
\end{aligned}
\label{eq:B111}
\end{equation}

\begin{equation}
\begin{aligned}
2T^{l+1}_{i+1,j,k}T^{l+1}_{i-1,j,k}&=2\left[T^{l}_{i+1, j,k}+(T^{l+1}_{i+1,j,k}-T^{l}_{i+1,j,k})\right]\left[T^{l}_{i-1, j,k}+(T^{l+1}_{i-1,j,k}-T^{l}_{i-1,j,k})\right] \\
&\approx 2\left[T^{l}_{i+1, j,k}T^{l}_{i-1,j,k}+T^{l}_{i+1, j,k}(T^{l+1}_{i-1,j,k}-T^{l}_{i-1,j,k})+(T^{l+1}_{i+1,j,k}-T^{l}_{i+1,j,k})T^{l}_{i-1,j,k}\right] \\
&=2\left[T^{l}_{i+1,j,k}T^{l+1}_{i-1,j,k}+T^{l}_{i-1,j,k}T^{l+1}_{i+1,j,k}-T^{l}_{i+1,j,k}T^{l}_{i-1,j,k}\right]
\end{aligned}
\label{eq:B112}
\end{equation}
Where we dropped the higher order terms..

We do the same for the $B_{22}$ and $B_{33}$ terms:
\begin{equation}
\left(T^{l+1}_{i,j \pm 1,k}\right)^{2} =2T^{l}_{i,j\pm 1,k}T^{l+1}_{i,j\pm 1,k}-(T^{l}_{i,j\pm 1,k})^{2}
\label{eq:B221}
\end{equation}

\begin{equation}
2T^{l+1}_{i,j+1,k}T^{l+1}_{i,j-1,k}=2\left[T^{l}_{i,j+1,k}T^{l+1}_{i,j-1,k}+T^{l}_{i,j-1,k}T^{l+1}_{i,j+1,k}-T^{l}_{i,j+1,k}T^{l}_{i,j-1,k}\right]
\label{eq:B222}
\end{equation}

\begin{equation}
\left(T^{l+1}_{i,j,k\pm 1}\right)^{2} =2T^{l}_{i, j,k\pm 1}T^{l+1}_{i, j,k\pm 1}-(T^{l}_{i,j,k\pm 1})^{2}
\label{eq:B331}
\end{equation}

\begin{equation}
2T^{l+1}_{i,j,k+1}T^{l+1}_{i,j,k-1}=2\left[T^{l}_{i,j,k+1}T^{l+1}_{i,j,k-1}+T^{l}_{i,j,k-1}T^{l+1}_{i,j,k+1}-T^{l}_{i,j,k+1}T^{l}_{i,j,k-1}\right]
\label{eq:B332}
\end{equation}
\section{Anisotropic Solver in AstroBEAR}

\subsection{Two Dimensional}
In 2D, if it's isotropic conductivity, the neighbor index are shown on the left below and those for the anisotropic conductivity case are shown on the right.\\

\begin{tikzpicture}[center_node/.style={rectangle,fill=red!50,draw,minimum size=1em,inner sep=3pt}]
\tikzstyle{neighbor}=[circle,draw=blue!40,fill=yellow!30,minimum size=1em]
\tikzstyle{neighbor2}=[circle,draw=red!30,fill=white,minimum size=1em]
\tikzstyle{neighbor3}=[circle,fill=blue!20,minimum size=1em]
%20\% blue

\begin{scope}[xshift=1cm] %2D isotropic
\node[center_node] (1) at (0,0) {0}; %no text
\node[neighbor] (2) at (-1,0) {1}; 
\node[neighbor] (3) at (1,0) {2};
\node[neighbor] (4) at (0,-1) {3};
\node[neighbor] (5) at (0,1) {4};
\node[neighbor3] (6) at (1,1) {};
\node[neighbor3] (7) at (-1,1) {};
\node[neighbor3] (8) at (-1,-1) {};
\node[neighbor3] (9) at (1,-1) {};

%\draw (1) -- (2);
%\draw (1) -- (3);
%\draw (1) -- (4);
%\draw (1) -- (5);
\draw (1) -- (2)--(7)--(5)--(1);
\draw (1) -- (3)--(6)--(5);
\draw (1) -- (4)--(8)--(2);
\draw (3) -- (9)--(4);
\end{scope}

\begin{scope}[xshift=8cm]  %2D anisotropic

\node[center_node] (1) at (0,0) {0}; %no text
\node[neighbor] (2) at (-1,0) {1}; 
\node[neighbor] (3) at (1,0) {2};
\node[neighbor] (4) at (0,-1) {3};
\node[neighbor] (5) at (0,1) {4};
\node[neighbor2] (6) at (1,1) {5};
\node[neighbor2] (7) at (-1,1) {6};
\node[neighbor2] (8) at (-1,-1) {7};
\node[neighbor2] (9) at (1,-1) {8};

\draw (1) -- (2)--(7)--(5)--(1);
\draw (1) -- (3)--(6)--(5);
\draw (1) -- (4)--(8)--(2);
\draw (3) -- (9)--(4);
\end{scope}
\end{tikzpicture}


So corresponding the offsets matrix in the code for 2D case is

\begin{lstlisting}
offsets(:,1)=(/-1,0/)
offsets(:,2)=(/1,0/)
offsets(:,3)=(/0,-1/)
offsets(:,4)=(/0,1/)
offsets(:,5)=(/1,1/)
offsets(:,6)=(/-1,1/)
offsets(:,7)=(/-1,-1/)
offsets(:,8)=(/1,-1/)
\end{lstlisting}


\subsection{Three Dimensional}

In 3D, if it's isotropic conductivity, the neighbor index are shown on the left below and those for the anisotropic conductivity case are shown on the right.\\

\begin{tikzpicture}
\tikzstyle{center_node}=[rectangle,fill=red!50,draw,minimum size=1em,inner sep=3pt]
\tikzstyle{neighbor}=[circle,draw=blue!40,fill=yellow!30,minimum size=1em]
\tikzstyle{neighbor2}=[circle,draw=red!30,fill=white,minimum size=1em]

\def \dx{1.5};
\def \dy{2};
\def \dz{1.5};
\def \nbx{3}; % this is y-axis in normal 3D coordinates
\def \nby{3}; % this is y-axis in normal 3D coordinates
\def \nbz{3}; % this is x-axis in normal 3D coordinates

\begin{scope}

% draw all the nodes in light color on a 3D lattice
\foreach \x in {1,...,\nbx} {
    \foreach \y in {1,...,\nby} {
        \foreach \z in {1,...,\nbz} {
            \node at (\x*\dx,\y*\dy,\z*\dz) [circle, fill=blue!20] {};
        }
    }
}


% z lines
\foreach \x in {1,...,\nbx} {
    \foreach \z in {1,...,\nbz}{
        \draw[dashed] (\x*\dx,\dy,\z*\dz) -- ( \x*\dx,\nby*\dy,\z*\dz);
    }
}

% x lines 
\foreach \y in {1,...,\nbx} {
    \foreach \z in {1,...,\nbz}{
        \draw[dashed] (\dx,\y*\dy,\z*\dz) -- ( \nbx*\dx,\y*\dy,\z*\dz);
    }
}

% y lines
\foreach \x in {1,...,\nbx} {
    \foreach \y in {1,...,\nbz}{
        \draw[dashed] (\x*\dx,\y*\dy,\dz) -- ( \x*\dx,\y*\dy,\nbz*\dz);
    }
}

% ====== draw solid lines on the top ======
\foreach \x in {1,...,\nbx} {
  \draw (\x*\dx,\nby*\dy,\dz) -- ( \x*\dx,\nby*\dy,2*\dz);
  \draw (\x*\dx,\nby*\dy,2*\dz) -- ( \x*\dx,\nby*\dy,3*\dz);
}

\foreach \z in {1,...,\nbz} {
  \draw (\dx,\nby*\dy,\z*\dz) -- ( 2*\dx,\nby*\dy,\z*\dz);
  \draw (2*\dx,\nby*\dy,\z*\dz) -- ( 3*\dx,\nby*\dy,\z*\dz);
}
% ========================

% ====== draw solid lines on the front =====
\foreach \x in {1,...,\nbx} {
  \draw (\x*\dx,\dy,\nbz*\dz) -- ( \x*\dx,2*\dy,\nbz*\dz);
  \draw (\x*\dx,2*\dy,\nbz*\dz) -- ( \x*\dx,3*\dy,\nbz*\dz);
}

\foreach \y in {1,...,\nby} {
  \draw (\dx,\y*\dy,\nbz*\dz) -- ( 2*\dx,\y*\dy,\nbz*\dz);
  \draw (2*\dx,\y*\dy,\nbz*\dz) -- ( 3*\dx,\y*\dy,\nbz*\dz);
}
%===============================

% ====== draw solid lines on the right =====
\foreach \z in {1,...,\nbz} {
  \draw (\nbx*\dx,\dy,\z*\dz) -- ( \nbx*\dx,2*\dy,\z*\dz);
  \draw (\nbx*\dx,2*\dy,\z*\dz) -- ( \nbx*\dx,3*\dy,\z*\dz);
}

\foreach \y in {1,...,\nby} {
  \draw (\nbx*\dx,\y*\dy,\dz) -- ( \nbx*\dx,\y*\dy,2*\dz);
  \draw (\nbx*\dx,\y*\dy,2*\dz) -- ( \nbx*\dx,\y*\dy,3*\dz);
}
%===============================
%re-draw center node and nearest neighbor nodes
\node[center_node] (0) at (2*\dx,2*\dy,2*\dz) {0}; % center point
\node[neighbor] (1) at (2*\dx,2*\dy,\dz) {1}; % normal coord,  x-1
\node[neighbor] (2) at (2*\dx,2*\dy,3*\dz) {2};  % normal coord, x+1
\node[neighbor] (3) at (\dx,2*\dy,2*\dz) {3}; % normal coord, y-1
\node[neighbor] (4) at (3*\dx,2*\dy,2*\dz) {4}; % normal coord, y+1
\node[neighbor] (5) at (2*\dx,\dy,2*\dz) {5}; % normal coord, z-1
\node[neighbor] (6) at (2*\dx,3*\dy,2*\dz) {6}; % normal coord, z+1

\end{scope}

\begin{scope} [xshift=7cm]

% draw all the nodes in light color on a 3D lattice
\foreach \x in {1,...,\nbx} {
    \foreach \y in {1,...,\nby} {
        \foreach \z in {1,...,\nbz} {
            \node at (\x*\dx,\y*\dy,\z*\dz) [circle, fill=blue!20] {};
        }
    }
}


% z lines
\foreach \x in {1,...,\nbx} {
    \foreach \z in {1,...,\nbz}{
        \draw[dashed] (\x*\dx,\dy,\z*\dz) -- ( \x*\dx,\nby*\dy,\z*\dz);
    }
}

% x lines 
\foreach \y in {1,...,\nbx} {
    \foreach \z in {1,...,\nbz}{
        \draw[dashed] (\dx,\y*\dy,\z*\dz) -- ( \nbx*\dx,\y*\dy,\z*\dz);
    }
}

% y lines
\foreach \x in {1,...,\nbx} {
    \foreach \y in {1,...,\nbz}{
        \draw[dashed] (\x*\dx,\y*\dy,\dz) -- ( \x*\dx,\y*\dy,\nbz*\dz);
    }
}

% ====== draw solid lines on the top ======
\foreach \x in {1,...,\nbx} {
  \draw (\x*\dx,\nby*\dy,\dz) -- ( \x*\dx,\nby*\dy,2*\dz);
  \draw (\x*\dx,\nby*\dy,2*\dz) -- ( \x*\dx,\nby*\dy,3*\dz);
}

\foreach \z in {1,...,\nbz} {
  \draw (\dx,\nby*\dy,\z*\dz) -- ( 2*\dx,\nby*\dy,\z*\dz);
  \draw (2*\dx,\nby*\dy,\z*\dz) -- ( 3*\dx,\nby*\dy,\z*\dz);
}
% ========================

% ====== draw solid lines on the front =====
\foreach \x in {1,...,\nbx} {
  \draw (\x*\dx,\dy,\nbz*\dz) -- ( \x*\dx,2*\dy,\nbz*\dz);
  \draw (\x*\dx,2*\dy,\nbz*\dz) -- ( \x*\dx,3*\dy,\nbz*\dz);
}

\foreach \y in {1,...,\nby} {
  \draw (\dx,\y*\dy,\nbz*\dz) -- ( 2*\dx,\y*\dy,\nbz*\dz);
  \draw (2*\dx,\y*\dy,\nbz*\dz) -- ( 3*\dx,\y*\dy,\nbz*\dz);
}
%===============================

% ====== draw solid lines on the right =====
\foreach \z in {1,...,\nbz} {
  \draw (\nbx*\dx,\dy,\z*\dz) -- ( \nbx*\dx,2*\dy,\z*\dz);
  \draw (\nbx*\dx,2*\dy,\z*\dz) -- ( \nbx*\dx,3*\dy,\z*\dz);
}

\foreach \y in {1,...,\nby} {
  \draw (\nbx*\dx,\y*\dy,\dz) -- ( \nbx*\dx,\y*\dy,2*\dz);
  \draw (\nbx*\dx,\y*\dy,2*\dz) -- ( \nbx*\dx,\y*\dy,3*\dz);
}
%===============================
%re-draw center node and nearest neighbor nodes
\node[center_node] (0) at (2*\dx,2*\dy,2*\dz) {0}; % center point

\node[neighbor] (1) at (2*\dx,2*\dy,\dz) {1}; % normal coord,  x-1
\node[neighbor] (2) at (2*\dx,2*\dy,3*\dz) {2};  % normal coord, x+1
\node[neighbor] (3) at (\dx,2*\dy,2*\dz) {3}; % normal coord, y-1
\node[neighbor] (4) at (3*\dx,2*\dy,2*\dz) {4}; % normal coord, y+1
\node[neighbor] (5) at (2*\dx,\dy,2*\dz) {5}; % normal coord, z-1
\node[neighbor] (6) at (2*\dx,3*\dy,2*\dz) {6}; % normal coord, z+1

\node[neighbor2] (7) at (2*\dx,1*\dy,1*\dz) {7}; % normal coord,  x-1
\node[neighbor2] (8) at (1*\dx,1*\dy,2*\dz) {8};  % normal coord, x+1
\node[neighbor2] (9) at (2*\dx,1*\dy,3*\dz) {9}; % normal coord, y-1
\node[neighbor2] (10) at (3*\dx,1*\dy,2*\dz) {10}; % normal coord, y+1

\node[neighbor2] (11) at (1*\dx,2*\dy,1*\dz) {11}; % normal coord, z-1
\node[neighbor2] (12) at (1*\dx,2*\dy,3*\dz) {12}; % normal coord, z+1
\node[neighbor2] (13) at (3*\dx,2*\dy,3*\dz) {13}; % normal coord,  x-1
\node[neighbor2] (14) at (3*\dx,2*\dy,1*\dz) {14};  % normal coord, x+1
\node[neighbor2] (15) at (2*\dx,3*\dy,1*\dz) {15}; % normal coord, y-1
\node[neighbor2] (16) at (1*\dx,3*\dy,2*\dz) {16}; % normal coord, y+1
\node[neighbor2] (17) at (2*\dx,3*\dy,3*\dz) {17}; % normal coord, z-1
\node[neighbor2] (18) at (3*\dx,3*\dy,2*\dz) {18}; % normal coord, z+1

\end{scope}
\end{tikzpicture}

So corresponding the offsets matrix in the code for 3D case is

\begin{lstlisting}
offsets(:,1)=(/-1,0,0/)
offsets(:,2)=(/1,0,0/)
offsets(:,3)=(/0,-1,0/)
offsets(:,4)=(/0,1,0/)
offsets(:,5)=(/0,0,-1/)
offsets(:,6)=(/0,0,1/)
offsets(:,7)=(/-1,0,-1/)
offsets(:,8)=(/0,-1,-1/)
offsets(:,9)=(/1,0,-1/)
offsets(:,10)=(/0,1,-1/)
offsets(:,11)=(/-1,-1,0/)
offsets(:,12)=(/1,-1,0/)
offsets(:,13)=(/1,1,0/)
offsets(:,14)=(/-1,1,0/)
offsets(:,15)=(/-1,0,1/)
offsets(:,16)=(/0,-1,1/)
offsets(:,17)=(/1,0,1/)
offsets(:,18)=(/0,1,1/)
\end{lstlisting}

In the code we define array neighbor(dir, edge, subdir, subedge) according the the diagram above. The subdir is defined using function $mod(l+q-1, nDim)+1$. So in 2D
\begin{equation}
subdir=\begin{cases} 2\mbox{ or y, for x(l=1,q=1) } \\
                     1\mbox{ or x, for y(l=2,q=2) }
     \end{cases}
\end{equation}
and in 3D, since every direction has two sub-directions (for example, x direction has y and z as sub-directions)
\begin{equation}
subdir=\begin{cases} 2\mbox{ or y, for x (l=1, q=1) }\\
                     3\mbox{ or z, for x (l=1, q=2) }\\
                     3\mbox{ or z, for y (l=2, q=1) }\\
                     1\mbox{ or x, for y (l=2, q=2) }\\
                     1\mbox{ or x, for z (l=3, q=1) }\\
                     2\mbox{ or y, for z (l=3, q=2) }
     \end{cases}
\end{equation}
\section{Multiphysics Integration}


\end{document}
