switch CTU to conservative diffusion

[zingale]

Currently we treat diffusion as an external source and compute the old-time diffusion from S_old, use this in the predictor, and then compute the new-time diffusion after the conservative update, from S_new. The final source is then 1/2(diffusion_old + diffusion_new). However, this is not conservative, and no correction is done at C-F interfaces.

We should instead compute the diffusive flux, -k_th grad T and add this to the energy flux. To be second order, we would still need the old-time diffusion in the predictor to the interface states, and we would also need to evaluate this flux from the Godunov state (including the transport coefficients). This would then automatically be conservative and participate in the flux correction.

This approach will also address #369

[zingale]

These are the steps that need to be done:

• [ ] we need to remove the diffusion source term from construct_old_source and construct_new_source (because we don't want this term included in the apply_source_to_state

• [ ] we still need to get the diffusion source and add it to old_source before the FillPatchAdd to sources_for_hydro, so it will be seen by the predictor for the time-centered fluxes.

• [ ] to include the diffusive flux in the energy flux, we would modify compute_flux_q. There we have the Riemann state on the interface to compute the thermal conductivity, but how to we compute grad T?

[zingale]

We can still use the mlmg operator to get the diffusive flux, see:

    /**
    * \brief For ``(alpha * a - beta * (del dot b grad)) phi = rhs``, flux means ``-b grad phi``
    *
    * \param a_flux
    * \param a_loc
    */
    void getFluxes (const Vector<Array<MultiFab*,AMREX_SPACEDIM> >& a_flux,
                    Location a_loc = Location::FaceCenter);
    void getFluxes (const Vector<Array<MultiFab*,AMREX_SPACEDIM> >& a_flux,
                    const Vector<MultiFab*> & a_sol,
                    Location a_loc = Location::FaceCenter);
    void getFluxes (const Vector<MultiFab*> & a_flux, Location a_loc = Location::CellCenter);
    void getFluxes (const Vector<MultiFab*> & a_flux, const Vector<MultiFab*> & a_sol, Location a_loc = Location::CellCenter);```

[zingale]

I think that this is the flowchart that works for CTU, MOL, and SDC:

1. at the start of a step (or substep) compute the "old source", as we do now. For CTU, store the diffusion term in the sources MF that will be used in the prediction. For all integration types, also get the flux and store this, either in a set of MFs for the diffusive flux (one for each coordinate direction) or by already adding it to the main flux array
2. (CTU only), the diffusion source term will be used in the prediction of the interface states as it is now
3. for MOL/SDC we will add the diffusive flux to the energy fluxes before the conservative update -- this will then correctly incorporate it as a diffusion source, and will also get store in the boundary registers for refluxing
4. for CTU, we will need to do a predictor corrector on the diffusion, by computing the fluxes after the energy update and them time-centering them for the final flux. This will likely be done in the new source term stuff, although we will not actually compute the new source term and add it