\rho c_v \frac{\partial T}{\partial t} = \nabla \left ( \kappa T^{n} \nabla T \right )

\rho c_v \frac{\partial T}{\partial t} =\nabla^2 \left ( \frac{1}{n+1} \kappa T^{n+1} \right )

For backward euler we evaluate the RHS at time j+1 although we still have to linearize

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 )

or equivalently

T_i^{{j+1}^{n+1}} = T_i^{j^n} \left (T_i^{j+1} + n \left (T_i^{j+1} - T_i^{j} \right ) \right )

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}

Now to make the scheme 2nd order in time we can set T_i^{j+1}=\frac{T_i^{j}+T_i^{j+1}}{2} and we get

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

A_\pm = -\frac{\kappa \Delta t}{2 \Delta x^2}T_{i\pm1}^{j^n}

and the center coefficient as

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

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

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 A and 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 c_v \rightarrow c_v \times \chi = \gamma_7 and \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

\rho c_v \frac{\partial T}{\partial t} = -\nabla Q

which can be discretized as

\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

A_0 = \rho c_v

and a left right source term contribution

B_\pm = \mp \frac{\Delta t}{\Delta x} Q_{i\pm 1/2}

and if the left boundary is a flux boundary, but the right side is normal, then we have to adjust the coefficients

A_+ = -\frac{\kappa \Delta t}{2 \Delta x^2}T_{i+1}^{j^n}

and the center coefficient as

A_0 = \frac{\kappa \Delta t}{ 2\Delta x^2}T_{i}^{j^n} + \rho c_v

and the RHS left and right contributions

B_+ = -\frac{\left (n - 1 \right ) \kappa \Delta t}{ 2 \left ( n+ 1 \right ) \Delta x^2} T_{i+1}^{j^{n+1}}

B_- = \frac{\Delta t}{\Delta x} Q_{i-1/2}

and the RHS center term contribution as

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