Changes between Version 5 and Version 6 of ThermalConduction

Timestamp:

10/13/15 13:20:05 (9 years ago)

Author:

Jonathan

Comment:

—

ThermalConduction

@@ -1 @@
-== Thermal Conduction Solver ==
-
-=== !!Update 5.22.2015 ===
-In the newer version of the code, the unit of temperature in the code or Info%q(iTemp) was changed to Kelvin from C.U.. However the following equations and source code keep the same except the term "gamma7*rho" has to be changed to function getdEdT() which gives C.U energy/Kelvin.. and no gamma7 appears...
-
-=== Before temperature unit was changed to Kelvin ===
-
-The equation we are trying to solve in AstroBEAR is
-
-{{{
-#!latex
-\begin{equation}
-\rho c_v \frac{\partial T}{\partial t} = \nabla \left ( \kappa T^{n} \nabla T \right )
-\end{equation}
-}}}
-
-which can be "Laplacetasized" as 
-
-{{{ 
-#!latex
-\setcounter{equation} 1
-\begin{equation} 
-\rho c_v \frac{\partial T}{\partial t} =\nabla^2 \left ( \frac{1}{n+1} \kappa T^{n+1} \right )
-\end{equation}
-}}}
-
-For backward euler we evaluate the RHS at time [[latex($j+1$)]] although we still have to linearize 
-
-[[latex($T_i^{{j+1}^{n+1}}$)]]
-
-So we can Taylor expand about [[latex($T_i^j$)]]
-
-
-[[latex($T_i^{{j+1}^{n+1}} = T_i^{j^{n+1}} + \left (n + 1 \right ) T_i^{j^n} \left ( T_i^{j+1} - T_i^j \right ) $)]]
-
-now if we substitute that back into the discretized form of equation 2 we get
-
-{{{
-#!latex
-\setcounter{equation} 2
-\begin{eqnarray}
-\rho_i c_v \left ( T_i^{j+1}-T_i^j \right ) & =  & \frac{\kappa \Delta t}{\Delta x^2}\left (T_{i+1}^{j^n}  T_{i+1}^{j+1} - 2 T_{i}^{j^n}  T_{i}^{j+1} + T_{i-1}^{j^n}  T_{i-1}^{j+1} \right ) \\
- & - & \frac{n \kappa \Delta t}{ \left ( n+ 1 \right ) \Delta x^2} \left (T_{i+1}^{j^{n+1}} - 2 T_{i}^{j^{n+1}} + T_{i-1}^{j^{n+1}} \right )  
-\end{eqnarray}
-}}}
-
-== Normal CK scheme ==
-We can identify the left and right coefficients as
-
-[[latex($A_\pm = -\frac{\kappa \Delta t}{\Delta x^2}T_{i\pm1}^{j^n} $)]]
-
-and the center coefficient as 
-
-[[latex($A_0 =  2 \frac{\kappa \Delta t}{ \Delta x^2}T_{i}^{j^n}  + \rho c_v$)]]
-
-and the RHS left and right contributions
-
-[[latex($B_\pm = -\frac{n \kappa \Delta t}{ \left ( n+ 1 \right ) \Delta x^2} T_{i\pm 1}^{j^{n+1}}$)]]
-
-and the RHS center term contribution as 
-
-[[latex($B_0 = 2 \frac{n \kappa \Delta t}{ \left ( n+ 1 \right ) \Delta x^2} T_{i}^{j^{n+1}} + \rho c_v T_i^j$)]]
-
-and finally in the code the values for [[latex($A$)]] and [[latex($B$)]] are negated:
-
-{{{
-    kx = kappa1*dt_diff/(dx**2)
-    CALL getTemp(Info,i,j,k,T)
-    stencil_fixed(0) = -2.0*REAL(nDim,8)*kx*T(0)**ndiff
-    DO l = 1, 2*nDim
-       stencil_fixed(l) = kx*T(l)**ndiff
-    END DO
-    source = ((ndiff)/(ndiff+1.0))*dot_product(T(0:2*nDim),stencil_fixed(0:2*nDim)) - T(0)*gamma7*Info%q(i,j,k,1)
-    stencil_fixed(0) = stencil_fixed(0) - gamma7*Info%q(i,j,k,1)
-}}}
-
-also note that the equation is multipled by the mean molecular weight... which is why [[latex($c_v \rightarrow c_v \times  \chi = \gamma_7$)]] and [[latex($\kappa \rightarrow \kappa \times \chi = \kappa_1$)]]
-
-
-== a note about flux boundary conditions ==
-If we have a specified flux at a boundary - then the equation is just 
-
-[[latex($\rho c_v \frac{\partial T}{\partial t} = -\nabla Q $)]]
-
-which can be discretized as 
-
-[[latex($\rho_i c_v \left (T_i^{j+1}-T_i^j \right ) = \frac{-\Delta t}{\Delta x} \left ( Q_{i+1/2} - Q_{i-1/2} \right ) $)]]
-
-which leads to a center coefficient of 
-
-[[latex($A_0 = \rho c_v$)]]
-
-and a left right source term contribution
-
-[[latex($B_\pm = \mp \frac{\Delta t}{\Delta x} Q_{i\pm 1/2}$)]]
-
-=== Now if the left boundary is a flux boundary, but the right side is normal, then we have to adjust the coefficients ===
-
-[[latex($A_+ = -\frac{\kappa \Delta t}{\Delta x^2}T_{i+1}^{j^n} $)]]
-
-[[latex($A_- = 0 $)]]
-
-and the center coefficient as 
-
-[[latex($A_0 =  \frac{2\kappa \Delta t}{ \Delta x^2}T_{i}^{j^n}  + \rho c_v$)]]
-
-and the RHS left and right contributions
-
-[[latex($B_+ = -\frac{n \kappa \Delta t}{ 2  \left ( n+ 1 \right ) \Delta x^2} T_{i+1}^{j^{n+1}}$)]]
-
-[[latex($B_-  =  \frac{\Delta t}{\Delta x} Q_{i-1/2}$)]]
-
-and the RHS center term contribution as 
-
-[[latex($B_0 = \frac{2n \kappa \Delta t}{ \left ( n+ 1 \right ) \Delta x^2} T_{i}^{j^{n+1}} + \rho c_v T_i^j$)]]
-
-
-Of course - remembering that in AstroBEAR everything is negated...
-
-This is accomplished by
-{{{
-     source = source + (ndiff/(ndiff+1.0))*(T(0)*kx*T(0)**ndiff)
-     source = source - (ndiff/(ndiff+1.0))*(T(p)*kx*T(p)**ndiff)
-     stencil_fixed(0)=stencil_fixed(0)+kx*T(0)**ndiff
-     stencil_fixed(p)=0.0
-     source = source - flb*dt_diff/dx
-}}}
-
-== If the right boundary is a zero-slope boundary, we have ==
-[[latex($T_{i+1}^{j+1}=T_{i}^{j+1}$)]]
-and 
-[[latex($T_{i+1}^{j}=T_{i}^{j}$)]]
-
-So now on the right boundary, the equation becomes
-{{{
-#!latex
-\setcounter{equation} 4
-\begin{eqnarray}
-\rho_i c_v \left ( T_i^{j+1}-T_i^j \right ) & =  & \frac{\kappa \Delta t}{\Delta x^2}\left (- T_{i}^{j^n}  T_{i}^{j+1} + T_{i-1}^{j^n}  T_{i-1}^{j+1} \right ) \\
- & - & \frac{n\kappa \Delta t}{ \left ( n+ 1 \right ) \Delta x^2} \left (-  T_{i}^{j^{n+1}} + T_{i-1}^{j^{n+1}} \right )  
-\end{eqnarray}
-}}}
-And the coeeficients becomes
-
-[[latex($A_+ = 0 $)]]
-
-[[latex($A_-=-\frac{\kappa \Delta t}{\Delta x^2}T_{i-1}^{j^n} $)]]
-
-[[latex($A_0=\frac{\kappa \Delta t}{\Delta x^2}T_{i}^{j^n}+ \rho c_v$)]]
-
-and the RHS left and right contributions
-
-[[latex($B_+=0 $)]]
-
-[[latex($B_-=-\frac{n \kappa \Delta t}{ \left ( n+ 1 \right ) \Delta x^2} T_{i\pm 1}^{j^{n+1}}$)]]
-
-[[latex($B_0=\frac{n \kappa \Delta t}{ \left ( n+ 1 \right ) \Delta x^2} T_{i}^{j^n}+\rho C_vT_{i}^{j}$)]]
-
-Again in AstroBEAR everything is negated
-
-This is accomplished by
-
-{{{
-stencil_fixed(0) =stencil_fixed(0)+kx*info%q(iq(1),iq(2),iq(3),itemp)**ndiff ! store the stencil value to temp
-ource = source -(ndiff)/(ndiff+1.0)*(stencil_fixed(p)*info%q(ip(1),ip(2),ip(3),itemp)+kx*info%q(iq(1),iq(2),iq(3),itemp)**(ndiff+1))
-stencil_fixed(p) = 0.0 ! zero out the stencil
-}}}
-
-== 2nd order in Time CK scheme ==
-Now to make the scheme 2nd order in time we can set [[latex($T_i^{j+1}=\frac{T_i^{j}+T_i^{j+1}}{2}$)]] and we get
-
-{{{
-#!latex
-\setcounter{equation} 4
-\begin{eqnarray}
-\rho_i c_v \left ( T_i^{j+1}-T_i^j \right ) & =  & \frac{\kappa \Delta t}{2 \Delta x^2}\left (T_{i+1}^{j^n}  T_{i+1}^{j+1} - 2 T_{i}^{j^n}  T_{i}^{j+1} + T_{i-1}^{j^n}  T_{i-1}^{j+1} \right ) \\
- & - & \frac{\left (n - 1 \right ) \kappa \Delta t}{ 2 \left ( n+ 1 \right ) \Delta x^2} \left (T_{i+1}^{j^{n+1}} - 2 T_{i}^{j^{n+1}} + T_{i-1}^{j^{n+1}} \right )  
-\end{eqnarray}
-}}}
-
-Where we can identify the left and right coefficients as
-
-[[latex($A_\pm = -\frac{\kappa \Delta t}{2 \Delta x^2}T_{i\pm1}^{j^n} $)]]
-
-and the center coefficient as 
-
-[[latex($A_0 =  2 \frac{\kappa \Delta t}{ 2\Delta x^2}T_{i}^{j^n}  + \rho c_v$)]]
-
-and the RHS left and right contributions
-
-[[latex($B_\pm = -\frac{\left (n - 1 \right ) \kappa \Delta t}{ 2  \left ( n+ 1 \right ) \Delta x^2} T_{i\pm 1}^{j^{n+1}}$)]]
-
-and the RHS center term contribution as 
-
-[[latex($B_0 = 2 \frac{\left (n - 1 \right ) \kappa \Delta t}{ 2  \left ( n+ 1 \right ) \Delta x^2} T_{i}^{j^{n+1}} + \rho c_v T_i^j$)]]
-
-and finally in the code the values for [[latex($A$)]] and [[latex($B$)]] are negated:
-
-{{{
-    kx = (0.5*kappa1*dt_diff/(dx**2))
-    CALL getTemp(Info,i,j,k,T)
-    stencil_fixed(0) = -2.0*REAL(nDim,8)*kx*T(0)**ndiff
-    DO l = 1, 2*nDim
-       stencil_fixed(l) = kx*T(l)**ndiff
-    END DO  
-    source = ((ndiff-1.0)/(ndiff+1.0))*dot_product(T(0:2*nDim),stencil_fixed(0:2*nDim)) - T(0)*gamma7*Info%q(i,j,k,1)
-    stencil_fixed(0) = stencil_fixed(0) - gamma7*Info%q(i,j,k,1)
-}}}
-
-also note that the equation is multipled by the mean molecular weight... which is why [[latex($c_v \rightarrow c_v \times  \chi = \gamma_7$)]] and [[latex($\kappa \rightarrow \kappa \times \chi = \kappa_1$)]]
-
-
-== a note about flux boundary conditions ==
-If we have a specified flux at a boundary - then the equation is just 
-
-[[latex($\rho c_v \frac{\partial T}{\partial t} = -\nabla Q $)]]
-
-which can be discretized as 
-
-[[latex($\rho_i c_v \left (T_i^{j+1}-T_i^j \right ) = \frac{-\Delta t}{\Delta x} \left ( Q_{i+1/2} - Q_{i-1/2} \right ) $)]]
-
-which leads to a center coefficient of 
-
-[[latex($A_0 = \rho c_v$)]]
-
-and a left right source term contribution
-
-[[latex($B_\pm = \mp \frac{\Delta t}{\Delta x} Q_{i\pm 1/2}$)]]
-
-=== Now if the left boundary is a flux boundary, but the right side is normal, then we have to adjust the coefficients ===
-
-[[latex($A_+ = -\frac{\kappa \Delta t}{2 \Delta x^2}T_{i+1}^{j^n} $)]]
-
-[[latex($A_- = 0 $)]]
-
-and the center coefficient as 
-
-[[latex($A_0 =  \frac{\kappa \Delta t}{ 2\Delta x^2}T_{i}^{j^n}  + \rho c_v$)]]
-
-and the RHS left and right contributions
-
-[[latex($B_+ = -\frac{\left (n - 1 \right ) \kappa \Delta t}{ 2  \left ( n+ 1 \right ) \Delta x^2} T_{i+1}^{j^{n+1}}$)]]
-
-[[latex($B_-  =  \frac{\Delta t}{\Delta x} Q_{i-1/2}$)]]
-
-and the RHS center term contribution as 
-
-[[latex($B_0 = \frac{\left (n - 1 \right ) \kappa \Delta t}{ 2  \left ( n+ 1 \right ) \Delta x^2} T_{i}^{j^{n+1}} + \rho c_v T_i^j$)]]
-
-
-Of course - remembering that in AstroBEAR everything is negated...
-
-This is accomplished by
-{{{
-     source = source + ((ndiff-1.0)/(ndiff+1.0))*(T(0)*kx*T(0)**ndiff)
-     source = source - ((ndiff-1.0)/(ndiff+1.0))*(T(p)*kx*T(p)**ndiff)
-     stencil_fixed(0)=stencil_fixed(0)+kx*T(0)**ndiff
-     stencil_fixed(p)=0.0
-     source = source - flb*dt_diff/dx
-}}}
-
-== If the right boundary is a zero-slope boundary, we have ==
-[[latex($T_{i+1}^{j+1}=T_{i}^{j+1}$)]]
-and 
-[[latex($T_{i+1}^{j}=T_{i}^{j}$)]]
-
-So now on the right boundary, the equation becomes
-{{{
-#!latex
-\setcounter{equation} 4
-\begin{eqnarray}
-\rho_i c_v \left ( T_i^{j+1}-T_i^j \right ) & =  & \frac{\kappa \Delta t}{2 \Delta x^2}\left (- T_{i}^{j^n}  T_{i}^{j+1} + T_{i-1}^{j^n}  T_{i-1}^{j+1} \right ) \\
- & - & \frac{\left (n - 1 \right ) \kappa \Delta t}{ 2 \left ( n+ 1 \right ) \Delta x^2} \left (-  T_{i}^{j^{n+1}} + T_{i-1}^{j^{n+1}} \right )  
-\end{eqnarray}
-}}}
-And the coeeficients becomes
-
-[[latex($A_+ = 0 $)]]
-
-[[latex($A_-=-\frac{\kappa \Delta t}{2 \Delta x^2}T_{i-1}^{j^n} $)]]
-
-[[latex($A_0=\frac{\kappa \Delta t}{2 \Delta x^2}T_{i}^{j^n}+ \rho c_v$)]]
-
-and the RHS left and right contributions
-
-[[latex($B_+=0 $)]]
-
-[[latex($B_-=-\frac{\left (n - 1 \right ) \kappa \Delta t}{ 2  \left ( n+ 1 \right ) \Delta x^2} T_{i\pm 1}^{j^{n+1}}$)]]
-
-[[latex($B_0=\frac{\left (n - 1 \right ) \kappa \Delta t}{ 2  \left ( n+ 1 \right ) \Delta x^2} T_{i}^{j^n}+\rho C_vT_{i}^{j}$)]]
-
-Again in AstroBEAR everything is negated
-
-This is accomplished by
-
-{{{
-stencil_fixed(0) =stencil_fixed(0)+kx*info%q(iq(1),iq(2),iq(3),itemp)**ndiff ! store the stencil value to temp
-!source = source -((ndiff-1.0)/(ndiff+1.0))*(stencil_fixed(p)*info%q(ip(1),ip(2),ip(3),itemp)-kx*info%q(iq(1),iq(2),iq(3),itemp)**(ndiff+1))
-source = source -((ndiff-1.0)/(ndiff+1.0))*(stencil_fixed(p)*info%q(ip(1),ip(2),ip(3),itemp)+kx*info%q(iq(1),iq(2),iq(3),itemp)**(ndiff+1))
-stencil_fixed(p) = 0.0 ! zero out the stencil
-}}}
-
-where the commented out line is the old-version which got the sign before the "kx" wrong.
-
-== Limiting case ==
-
-In the limit of [[latex($T \rightarrow 0$)]] with a flux from the left boundary:
-
-our equation becomes [[latex($\rho_i c_v T_i^{j+1} = \rho_i c_v T_i^j + \frac{\Delta t}{\Delta x} Q \delta_i^1$)]]
-
-and the total energy [[latex($E^{j+1}=\rho c_v \sum T_i^{j+1} \Delta x = E^{j} + \Delta t Q$)]]
-
-or [[latex($\dot{E} = Q$)]]
+= 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}}$
+
+== 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 T^{\lambda_\parallel} - \frac{n^2 \kappa_\perp}{B^2} T^{\lambda_\perp} \right ) \left ( \hat{b} \cdot \nabla T \right ) + \frac{n^2 \kappa_\perp}{B^2} T^{\lambda_\perp} \nabla T \right ]$
+
+or 
+
+$\rho c_v \frac{\partial T}{\partial t} = \nabla \cdot \left [ \hat{b} \left ( \frac{\kappa_\parallel}{\lambda_\parallel+1} \left ( \hat{b} \cdot \nabla T^{\lambda_\parallel+1} \right )  - \frac{n^2 \kappa_\perp}{B^2 \left ( \lambda_\perp+1 \right )} \left ( \hat{b} \cdot \nabla T^{\lambda_\perp + 1} \right )  \right ) + \frac{n^2 \kappa_\perp}{B^2 \left ( \lambda_\perp +1 \right )} \nabla T^{\lambda_\perp+1} \right ]$
+
+== Einstein simplification ==
+
+or in Einstein notation
+
+$\rho c_v \partial_t T = \partial_i \left [ b_i \left ( \frac{\kappa_\parallel}{\lambda_\parallel+1} \left ( b_j \partial_j T^{\lambda_\parallel+1} \right )  - \frac{n^2 \kappa_\perp}{B^2 \left ( \lambda_\perp+1 \right )} \left ( b_j \partial_j  T^{\lambda_\perp + 1} \right )  \right ) + \frac{n^2 \kappa_\perp}{B^2 \left ( \lambda_\perp +1 \right )} \partial_i T^{\lambda_\perp+1} \right ]$
+
+or 
+
+$\rho c_v  \partial_t T = \partial_i \left [ b_i \left ( \frac{\kappa_\parallel}{\lambda_\parallel+1} \left ( b_j \partial_j T^{\lambda_\parallel+1} \right )  - \frac{n^2 \kappa_\perp}{B^2 \left ( \lambda_\perp+1 \right )} \left ( b_j \partial_j  T^{\lambda_\perp + 1} \right )  \right ) + \frac{n^2 \kappa_\perp}{B^2 \left ( \lambda_\perp +1 \right )} \delta_{ij} \partial_j T^{\lambda_\perp+1} \right ]$
+
+or 
+
+$\rho c_v \partial_t T = \partial_i b_i \frac{\kappa_\parallel}{\lambda_\parallel+1} b_j \partial_j T^{\lambda_\parallel+1} +   \partial_i \left ( \delta_{ij} - b_i b_j \right ) \frac{n^2 \kappa_\perp}{B^2 \left ( \lambda_\perp + 1\right )} \partial_j T^{\lambda_\perp+1} $
+
+or 
+
+$\rho c_v  \partial_t T = \frac{\kappa_\parallel}{\lambda_\parallel+1}  \partial_i b_i  b_j \partial_j T^{\lambda_\parallel+1} +   \frac{\kappa_\perp}{\left ( \lambda_\perp + 1\right )} \partial_i n^2 B^{-2} \left ( \delta_{ij} - b_i b_j \right )  \partial_j T^{\lambda_\perp +1} $
+
+== Implicitization ==
+
+Now we can rewrite the equation 
+
+$\partial_t T = \displaystyle \sum_{\parallel, \perp} A \partial_i  B_{ij} \partial_j T^{\lambda+1}$
+
+where the $\parallel$ or $\perp$ subscript on $A$, $B_{ij}$, and $\lambda$ is implied.  
+
+$A_\parallel = \frac{\kappa_\parallel}{\rho c_v  \left ( \lambda_\parallel + 1\right ) }$ and $A_\perp = \frac{\kappa_\perp}{\rho c_v  \left (\lambda_\perp + 1 \right )}$
+
+and $B_{\parallel, ij} = b_i b_j$ and $B_{\perp, ij} = n^2 B^{-2} \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_*^{\lambda+1}$
+
+It is simpler to linearize the equation in terms of $T_*$ and then subsitute then vice-versa - though both give the same answer
+
+$T_*^{\lambda + 1} \approx T^{\lambda + 1} +  \left ( \lambda + 1 \right ) T^{\lambda} \left ( T_{*} - T \right ) =  T^{\lambda + 1} +  \left ( \lambda + 1 \right ) T^{\lambda} \phi \left ( T' - T \right )$
+$= \left(1 - \phi \left ( \lambda + 1\right ) \right ) T^{\lambda + 1}+  \phi \left ( \lambda + 1 \right ) T^{\lambda} T' $
+
+
+So we have
+
+$\partial_t T = \displaystyle \sum_{\parallel, \perp}{A \partial_i  B_{ij} \partial_j \left [ \left ( 1 - \phi \left ( \lambda + 1 \right ) \right )T^{\lambda + 1} + \phi \left ( \lambda + 1 \right ) T^{\lambda} T' \right ]}$
+
+$\partial_t T = \displaystyle \sum_{\parallel, \perp}{A \partial_i  B_{ij} \left [ \left ( 1 - \phi \left ( \lambda + 1 \right ) \right ) \partial_j T^{\lambda + 1} + \phi \left ( \lambda + 1 \right ) \partial_j T^{\lambda} T' \right ]}$
+
+$\partial_t T = \displaystyle \sum_{\parallel, \perp}{A \partial_i  B_{ij} \left [ C \partial_j T^{\lambda + 1} + D \partial_j T^{\lambda} T' \right ]}$
+
+where $ C= \left ( 1 - \phi \left ( \lambda + 1 \right ) \right )$ and $ D = \phi \left ( \lambda + 1 \right )$
+
+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^{\lambda+1} = \left [ \frac{B^{ij}_{\hat{i}} T^{\lambda+1}_{\hat{i}+\hat{j}} - B^{ij}_{\hat{i}}T^{\lambda+1}_{\hat{i}-\hat{j}}-B^{ij}_{-\hat{i}} T^{\lambda+1}_{-\hat{i}+\hat{j}} + B^{ij}_{-\hat{i}}T^{\lambda+1}_{-\hat{i}-\hat{j}}}{4 \Delta x^2} \right ] \left ( 1-\delta_{ij} \right )$
+$ + \left [ \frac{B^{ij}_{\hat{i/2}} T^{\lambda+1}_{\hat{i}} - B^{ij}_{\hat{i/2}}T^{\lambda+1}_{0}-B^{ij}_{-\hat{i/2}} T^{\lambda+1}_{0} + B^{ij}_{-\hat{i/2}}T^{\lambda+1}_{\hat{-i}}}{\Delta x^2} \right ] \delta_{ij} $
+
+$ \partial_i B_{ij} \partial_j T^{\lambda}T' =  \left [ \frac{B^{ij}_{\hat{i}} T^{\lambda}_{\hat{i}+\hat{j}}T'_{\hat{i}+\hat{j}} - B^{ij}_{\hat{i}}T^{\lambda}_{\hat{i}-\hat{j}}T'_{\hat{i}-\hat{j}}-B^{ij}_{-\hat{i}} T^{\lambda}_{-\hat{i}+\hat{j}}T'_{-\hat{i}+\hat{j}} + B^{ij}_{-\hat{i}}T^{\lambda}_{-\hat{i}-\hat{j}}T'_{-\hat{i}-\hat{j}}}{4 \Delta x^2} \right ]  \left ( 1-\delta_{ij} \right )  $
+$+ \left [ \frac{B^{ij}_{\hat{i/2}} T^{\lambda}_{\hat{i}}T'_{\hat{i}} - B^{ij}_{\hat{i/2}}T^{\lambda}_{0}T'_{0}-B^{ij}_{-\hat{i/2}} T^{\lambda}_{0}T'_{0} + B^{ij}_{-\hat{i/2}}T^{\lambda}_{-\hat{i}}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^{\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 - \Delta t \displaystyle \sum_{\parallel, \perp} AD \left [  \alpha_0 T^{\lambda}_0T'_0 + \displaystyle \sum_{\pm j} \alpha_{\pm j} T^{\lambda}_{\pm \hat{j}} T'_{\pm \hat{j}} + \sum_{\pm i, \pm j,i \ne j} \alpha_{\pm i, \pm j} T^{\lambda}_{\pm \hat{i} \pm \hat{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{\kappa_\parallel}{\rho c_v  \left ( \lambda_\parallel + 1 \right ) }$ || $\frac{\kappa_\perp}{\rho c_v  \left (\lambda_\perp + 1\right )}$ ||
+|| $B_{ij} $ || $b_i b_j$ || $n^2 B^{-2} \left (\delta_{ij} - b_i b_j \right )$ ||
+|| $C$ ||  $\left ( 1 - \phi \left ( \lambda_\parallel + 1 \right ) \right )$ ||  $ \left ( 1 - \phi \left ( \lambda_\perp + 1 \right ) \right )$ ||
+|| $D$ ||  $\phi \left ( \lambda_\parallel + 1 \right )$ ||  $ \phi \left ( \lambda_\perp + 1 \right ) $||
+
+||$\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 T^\lambda \nabla T$
+
+or
+
+$\rho c_v \frac{\partial T}{\partial t} = \nabla \cdot \frac{\kappa}{\lambda+1} \nabla T^{\lambda+1}$
+
+or in Einstein notation
+
+$\rho c_v \partial_t T = \partial_i \frac{\kappa}{\lambda+1} \partial_i T^{\lambda+1}$
+
+which we can rewrite as
+
+$\partial_t T = A \partial_i B \partial i T^{\lambda+1}$
+
+$\partial_t T = A \partial_i B \delta_{ij} \partial j T^{\lambda+1}$
+
+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.
+
+
+||$\alpha_0$ || $- \displaystyle \sum_{i} \frac{B_{\hat{i}} + B_{-\hat{i}}}{\Delta x^2}$ ||
+||$\alpha_{\pm i}$ || $ \frac{B_{\pm \hat{i}}}{\Delta x^2}$ ||
+|| $\alpha*_{\pm i \pm j}$ || $ 0$ ||
+
+Also if $\kappa$ is a constant, then we have
+
+||$\alpha_0$ || $- \displaystyle \sum_{i} \frac{2B}{\Delta x^2}$ ||
+||$\alpha_{\pm i}$ || $ \frac{B}{\Delta x^2}$ ||
+|| $\alpha*_{\pm i \pm j}$ || $ 0$ ||
+
+
+
+=== 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.