reflux: the viscous reflux of a non-conservative Laplacian_S tracer is divided by rho_half, although its diffusion equation carries no density #206

Location: Source/NavierStokes.cpp:1786
Severity: High — wrong sync correction for the default tracer in shipped Exec/*/regtest.* hotspot decks (tracer diffusivity 0.001, density 2, two levels, reflux on); the effect is a coarse/fine correction term, so the magnitude is modest
Category: AMR-sync
Based on commit bb697bf5 (line numbers refer to that tree).

Problem

reflux divides the viscous flux-register correction of every NonConservative scalar by the half-time
density before the advective reflux is added:

// Source/NavierStokes.cpp:1786-1792
    for (int istate = AMREX_SPACEDIM; istate < NUM_STATE; istate++)
    {
      if (advectionType[istate] == NonConservative)
      {
          MultiFab::Divide(Ssync,Rh,0,istate-AMREX_SPACEDIM,1,0);
      }
    }

That is right for Temp, whose diffusion form is RhoInverse_Laplacian_S (rho dT/dt = div k grad T,
rho_flag = 1): the register holds dtarea(k grad T), so the increment of T is that amount over rho.
The default tracer is also NonConservative, but NS_setup.cpp:320-321 gives it
diffusionType[Tracer] = Laplacian_S (rho_flag = 0): scalar_diffusion_update solves
dS/dt = div beta grad S with no density anywhere, scalar_advection (line 805) and
MacProj::mac_sync_compute treat Ssync[Tracer] as a rate of S without rho, and the sync solve in
mac_sync (Diffusion::diffuse_scalar, rho_flag = 0) adds Ssync*dt to the RHS unscaled. Only the
viscous reflux part of Ssync[Tracer] is divided by rho_half, so it is wrong by a factor 1/rho_half.

Impact

Any run with amr.max_level > 0, ns.do_reflux = 1 (default), ns.scal_diff_coefs > 0, the default
ns.do_cons_trac = 0, and density different from 1: the coarse/fine correction of the tracer's diffusive
flux is scaled by 1/rho, so the tracer is not conserved across level boundaries and the coarse and fine
answers disagree by O(dt*(1-1/rho)*flux mismatch). Shipped: Exec/run2d/regtest.2d.hotspot,
Exec/run3d/regtest.3d.hotspot, Exec/eb_run2d/regtest.2d.hotspot, Exec/eb_run3d/regtest.3d.hotspot
and Tutorials/Bubble/inputs.2d.bubble* (all prob.density_ic = 2.0, scal_diff_coefs = 0.001,
max_level = 2). Runs with rho == 1 or do_cons_trac = 1 are unaffected. Every build.

Suggested fix

Divide by rho_half only for the diffusion form that has rho in front of the time derivative.

--- a/Source/NavierStokes.cpp
+++ b/Source/NavierStokes.cpp
@@ -1783,9 +1783,14 @@
       }
     }
 
+    //
+    // The viscous register holds the flux of the equation that was solved:
+    // only RhoInverse_Laplacian_S (rho dS/dt = div beta grad S, e.g. Temp)
+    // needs the increment divided by rho. Laplacian_S tracers carry no rho.
+    //
     for (int istate = AMREX_SPACEDIM; istate < NUM_STATE; istate++)
     {
-      if (advectionType[istate] == NonConservative)
+      if (diffusionType[istate] == RhoInverse_Laplacian_S)
       {
           MultiFab::Divide(Ssync,Rh,0,istate-AMREX_SPACEDIM,1,0);
       }