ComputeDt explicit-diffusion limit uses the constant m_mu: ignores the non-Newtonian viscosity and the tracer/temperature diffusivities, so diffusion_type=0 gets an unstable dt

[WeiqunZhang]

Severity: medium · Category: physics-numerics
Locations: src/incflo_compute_dt.cpp:73-90, src/incflo_compute_dt.cpp:132-146, src/incflo_compute_dt.cpp:19, src/rheology/incflo_rheology.cpp:30-40, src/rheology/incflo_rheology.cpp:140-154
Based on commit 46de3367 (line numbers refer to that tree).

The defect

The header of ComputeDt states the diffusive term as V = 2 * max(eta/rho) * (1/dx^2 + ...)
(compute_dt.cpp:19). The implementation, in both the EB branch (73-90) and the regular branch (132-146),
computes max(1/rho) and then does diff_lev *= m_mu; -- the constant Newtonian coefficient read from
incflo.mu. It never evaluates the viscosity actually used by the explicit viscous term:

• compute_viscosity_at_level (rheology.cpp:68-138) returns mu + tau_0*expterm(sr/papa_reg)/papa_reg for
  Bingham, mu*sr^(n-1) for power-law, etc. For Bingham at low strain rate this is mu + tau_0/papa_reg, i.e.
  1001 for the shipped poiseuille_plane_bingham parameters (mu=1, tau_0=1, papa_reg=1e-3) instead of 1.
• The tracer diffusivities m_mu_s[n] (rheology.cpp:140-147) and the thermal diffusivity m_mu_T/(rho*cp)
  (rheology.cpp:149-154) are diffused explicitly by exactly the same m_diff_type switch
  (tracer_explicit_update, update_temperature) but do not enter the dt at all.

The state needed to evaluate the real coefficient is available: PR #209 moved the fillpatch of velocity and
density to the top of Advance() before ComputeDt, so velocity has its ghost cells filled when
ComputeDt runs, which is what the strain-rate evaluation needs.

Why it matters

With incflo.diffusion_type = 0 the code is forward-Euler in the diffusive terms and the dt limit is the only
thing keeping it stable. Forward Euler requires dt * nu * sum_d(1/dx_d^2) <= 1/2. The code's normalisation
(diff_cfl = 2*nu*sum(1/dx^2), dt = 2*cfl/comb_cfl) yields dt * nu * sum(1/dx^2) = cfl/2 when convection
and forcing are negligible, i.e. 0.245 for the shipped cfl = 0.49: a safety factor of about two, relative to
m_mu only. Any diffusivity larger than 2*m_mu therefore makes the run blow up. For a yield-stress fluid the
under-estimate is three orders of magnitude (eta(sr->0) = 1001 vs mu = 1); for the shipped Boussinesq bubble
inputs mu_s = 3e-5 is 3x mu = 1e-5, a tracer diffusion number of 3*0.245 = 0.735 > 0.5.

How to reach it

Any build (2D/3D, EB or not, CPU or GPU).

• test_2d/benchmark.poiseuille_plane_bingham with incflo.diffusion_type = 0 added and incflo.fixed_dt
  removed: dx = 2/32, m_mu = 1 gives dt ~ 20.5/(212256) ~ 1e-3, while the Bingham
  eta(sr -> 0) = 1001 requires dt ~ 1e-6. The explicit viscous update diverges within a few steps.
• test_no_eb_3d/benchmark.bouss_bubble_god with incflo.diffusion_type = 0 and incflo.fixed_dt removed:
  the tracer is advanced explicitly with mu_s = 3e-5 but dt is chosen from mu = 1e-5, a diffusive number of
  3*0.245 = 0.735 > 0.5 for the tracer.
  (Runs that set incflo.fixed_dt bypass the computed dt and only print the CFL warning, which is itself computed
  from the same wrong bound.)

Suggested fix

Compute the effective explicit diffusivity per cell -- the strain-rate dependent eta from
compute_viscosity_at_level, together with max_n mu_s[n] and mu_T/cp -- divide by rho (a non-conservative
tracer diffuses with mu_s itself, hence the max(1, 1/rho) factor), zero it in covered cells, and take its
max in place of m_mu * max(1/rho).

--- a/src/incflo_compute_dt.cpp
+++ b/src/incflo_compute_dt.cpp
@@ -51,6 +51,50 @@
 
        compute_vel_forces_on_level (lev, vel_forces, vel, rho, tra_o, tra);
 
+       // Explicit-diffusion bound: the largest diffusivity that is applied
+       // explicitly, divided by rho.  Use the strain-rate dependent viscosity
+       // (not the constant m_mu) and cover the tracer and temperature
+       // diffusivities too.  Covered cells are set to zero.
+       MultiFab nu;
+       if (explicit_diffusion) {
+           nu.define(grids[lev], dmap[lev], 1, 0, MFInfo(), Factory(lev));
+           compute_viscosity_at_level(lev, &nu, &m_leveldata[lev]->density,
+                                      &m_leveldata[lev]->velocity, geom[lev],
+                                      m_cur_time, 0);
+           Real mu_s_max = Real(0.0);
+           if (m_advect_tracer) {
+               for (int n = 0; n < m_ntrac; ++n) {
+                   mu_s_max = amrex::max(mu_s_max, m_mu_s[n]);
+               }
+           }
+           Real mu_T_eff = m_use_temperature ? m_mu_T / m_cp : Real(0.0);
+#ifdef AMREX_USE_EB
+           auto const& flags = EBFactory(lev).getMultiEBCellFlagFab();
+#endif
+#ifdef _OPENMP
+#pragma omp parallel if (Gpu::notInLaunchRegion())
+#endif
+           for (MFIter mfi(nu,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
+               Box const& bx = mfi.tilebox();
+               Array4<Real> const& nu_a = nu.array(mfi);
+               Array4<Real const> const& r = rho.const_array(mfi);
+#ifdef AMREX_USE_EB
+               Array4<EBCellFlag const> const& f = flags.const_array(mfi);
+#endif
+               ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+               {
+#ifdef AMREX_USE_EB
+                   if (f(i,j,k).isCovered()) { nu_a(i,j,k) = Real(0.0); return; }
+#endif
+                   Real rinv = Real(1.0)/r(i,j,k);
+                   // velocity and temperature diffuse with eta/rho and mu_T/(rho cp);
+                   // a non-conservative tracer diffuses with mu_s itself.
+                   nu_a(i,j,k) = amrex::max(amrex::max(nu_a(i,j,k), mu_T_eff)*rinv,
+                                            mu_s_max*amrex::max(Real(1.0), rinv));
+               });
+           }
+       }
+
 #ifdef AMREX_USE_EB
         if (!vel.isAllRegular()) {
             auto const& flag = EBFactory(lev).getMultiEBCellFlagFab();
@@ -71,22 +115,7 @@
                            return mx;
                        });
             if (explicit_diffusion) {
-                diff_lev = amrex::ReduceMax(rho, flag, 0,
-                           [=] AMREX_GPU_HOST_DEVICE (Box const& b,
-                                                      Array4<Real const> const& r,
-                                                      Array4<EBCellFlag const> const& f) -> Real
-                          {
-                              Real mx = Real(-1.0);
-                              amrex::Loop(b, [=,&mx] (int i, int j, int k) noexcept
-                              {
-                                  if (!f(i,j,k).isCovered()) {
-                                      Real rho_inv = Real(1.0)/r(i,j,k);
-                                      mx = amrex::max(rho_inv, mx);
-                                  }
-                              });
-                              return mx;
-                          });
-                diff_lev *= m_mu;
+                diff_lev = nu.max(0, 0, true);
             }
 
             // Forcing term -- old way of computing
@@ -130,19 +159,7 @@
                        });
 
             if (explicit_diffusion) {
-                diff_lev = amrex::ReduceMax(rho, 0,
-                           [=] AMREX_GPU_HOST_DEVICE (Box const& b,
-                                                      Array4<Real const> const& r) -> Real
-                           {
-                               Real mx = Real(-1.0);
-                               amrex::Loop(b, [=,&mx] (int i, int j, int k) noexcept
-                               {
-                                   Real rho_inv = Real(1.0)/r(i,j,k);
-                                   mx = amrex::max(rho_inv, mx);
-                               });
-                               return mx;
-                           });
-                diff_lev *= m_mu;
+                diff_lev = nu.max(0, 0, true);
             }
 
             // Forcing term -- old way of computing

(audit/notes/T2-scratch/005.diff, checked with git apply --check; not compiled.) nu.max(0,0,true) is a
local max like the ReduceMax it replaces; the existing ParallelAllReduce::Max below still does the MPI
reduction. With explicit_diffusion == false the block is skipped and nothing changes.

[asalmgren]

closed by #265