dissipation rate source terms for manifold variable variances

Source/PeleLMeX.H

@@ -827,6 +827,8 @@ public:
 
   void addSpark(const PeleLM::TimeStamp& a_time);
 
+  void addScalarDissipation(const PeleLM::TimeStamp& a_time);
+
   //-----------------------------------------------------------------------------
 
 #ifdef PELE_USE_SPRAY
@@ -1802,6 +1804,7 @@ public:
   int m_plotChemDiag = 0;
   int m_plotHeatRelease = 0;
   int m_plotStateSpec = 0;
+  int m_plotExtSource = 0;
   int m_plot_int = 0;
   bool m_plot_overwrite = false;
   int m_plot_zeroEBcovered = 1;
@@ -1911,7 +1914,7 @@ public:
   amrex::Real m_les_cs_smag = 0.18;
   amrex::Real m_les_cm_wale = 0.60;
   amrex::Real m_les_cs_sigma = 1.35;
-  amrex::Vector<amrex::Real> m_turb_visc_time;
+  amrex::Real m_les_c_chi = 20.0;
 
   // switch on/off different processes
   // TODO: actually implement that

Source/PeleLMeX_Forces.cpp

@@ -237,6 +237,82 @@ PeleLM::addSpark(const TimeStamp& a_timestamp)
   }
 }
 
+// Manifold model - dissipation rate sources for variances
+void
+PeleLM::addScalarDissipation(const TimeStamp& a_timestamp)
+{
+  BL_PROFILE("PeleLM::addScalarDissipationRate");
+  // no scalar dissipation sources if not using a manifold model
+#ifndef USE_MANIFOLD_EOS
+  return;
+#else
+
+  // Dissipation source term for subfilter variances
+  for (int lev = 0; lev <= finest_level; lev++) {
+
+    auto* ldata_p = getLevelDataPtr(lev, a_timestamp);
+    auto const& leosparm = eos_parms.host_parm();
+
+    // For now - we require the turbulent viscosity to be pre-computed
+    // it always is stored at AmrOldTime, so we just use that
+    // We need cell-centered mu_t but have it at faces
+    // The simple interpolation below probably isn't valid for EB
+#ifdef AMREX_USE_EB
+    amrex::Abort(
+      "PeleLM::addScalarDissipation(): this is not supported with EB");
+#endif
+
+    for (int n = 0; n < MANIFOLD_DIM; ++n) {
+      if (leosparm.is_variance_of[n] >= 0) {
+        if (!m_do_les) {
+          amrex::Abort("PeleLM::addScalarDissipation(): cannot add a "
+                       "scalarDissipation without an active LES model");
+        }
+
+        constexpr amrex::Real fact = 0.5 / AMREX_SPACEDIM;
+        const amrex::Real C_chi = m_les_c_chi;
+
+        AMREX_D_TERM(auto const& mut_arr_x =
+                       m_leveldata_old[lev]->visc_turb_fc[0].const_arrays();
+                     , auto const& mut_arr_y =
+                         m_leveldata_old[lev]->visc_turb_fc[1].const_arrays();
+                     , auto const& mut_arr_z =
+                         m_leveldata_old[lev]->visc_turb_fc[2].const_arrays();)
+        auto extma = m_extSource[lev]->arrays();
+        auto statema = ldata_p->state.const_arrays();
+
+        // l_scale will also need modification for EB
+        const amrex::Real vol = AMREX_D_TERM(
+          geom[lev].CellSize(0), *geom[lev].CellSize(1),
+          *geom[lev].CellSize(2));
+        const amrex::Real l_scale =
+          (AMREX_SPACEDIM == 2) ? std::sqrt(vol) : std::cbrt(vol);
+        const amrex::Real inv_l_scale2 = 1.0 / (l_scale * l_scale);
+
+        amrex::ParallelFor(
+          *m_extSource[lev],
+          [=] AMREX_GPU_DEVICE(int box_no, int i, int j, int k) noexcept {
+            amrex::Real mu_t =
+              fact *
+              (AMREX_D_TERM(
+                mut_arr_x[box_no](i, j, k) + mut_arr_x[box_no](i + 1, j, k),
+                +mut_arr_y[box_no](i, j, k) + mut_arr_y[box_no](i, j + 1, k),
+                +mut_arr_z[box_no](i, j, k) + mut_arr_z[box_no](i, j, k + 1)));
+
+            // Linear Relaxation model
+            // rho chi_sgs = C_chi * mu_t / Delta^2 * Variance
+            extma[box_no](i, j, k, FIRSTSPEC + n) -=
+              C_chi * mu_t * inv_l_scale2 *
+              statema[box_no](i, j, k, FIRSTSPEC + n);
+          });
+      }
+    }
+    Gpu::streamSynchronize();
+  }
+
+#endif
+}
+
 // Calculate additional external sources (soot, radiation, user defined, etc.)
 void
 PeleLM::getExternalSources(
@@ -268,6 +344,8 @@ PeleLM::getExternalSources(
   }
 #endif
 
+  addScalarDissipation(a_timestamp_old);
+
   // User defined external sources
   if (m_user_defined_ext_sources) {
     for (int lev = 0; lev <= finest_level; lev++) {
@@ -279,4 +357,4 @@ PeleLM::getExternalSources(
         ldata_p_new->state, ext_src, geom[lev].data(), *prob_parm_d);
     }
   }
-}
+}

Source/PeleLMeX_Plot.cpp

@@ -165,6 +165,10 @@ PeleLM::WritePlotFile()
     ncomp += 1;
   }
 
+  if (m_plotExtSource) {
+    ncomp += NVAR;
+  }
+
   //----------------------------------------------------------------
   // Plot MultiFabs
   Vector<MultiFab> mf_plt(finest_level + 1);
@@ -282,6 +286,12 @@ PeleLM::WritePlotFile()
     plt_VarsName.push_back("viscturb");
   }
 
+  if (m_plotExtSource) {
+    for (int ivar = 0; ivar < NVAR; ++ivar) {
+      plt_VarsName.push_back("extsource_" + stateVariableName(ivar));
+    }
+  }
+
 #if NUM_ODE > 0
   for (int n = 0; n < NUM_ODE; n++) {
     plt_VarsName.push_back(m_ode_names[n]);
@@ -437,6 +447,12 @@ PeleLM::WritePlotFile()
         });
       Gpu::streamSynchronize();
     }
+    cnt += 1;
+
+    if (m_plotExtSource) {
+      MultiFab::Copy(mf_plt[lev], *m_extSource[lev], 0, cnt, NVAR, 0);
+      cnt += NVAR;
+    }
 
 #ifdef AMREX_USE_EB
     if (m_plot_zeroEBcovered != 0) {

Source/PeleLMeX_Setup.cpp

@@ -413,9 +413,7 @@ PeleLM::readParameters()
     m_les_verbose = m_verbose;
     pp.query("plot_les", m_plot_les);
     pp.query("les_v", m_les_verbose);
-    for (int lev = 0; lev <= max_level; ++lev) {
-      m_turb_visc_time.push_back(-1.0E200);
-    }
+    pp.query("les_c_chi", m_les_c_chi);
 #ifdef PELE_USE_EFIELD
     amrex::Abort("LES implementation is not yet compatible with efield/ions");
 #endif
@@ -837,6 +835,7 @@ PeleLM::readIOParameters()
   }
   pp.query("plot_zeroEBcovered", m_plot_zeroEBcovered);
   pp.query("plot_speciesState", m_plotStateSpec);
+  pp.query("plot_extSource", m_plotExtSource);
   m_initial_grid_file = "";
   m_regrid_file = "";
   pp.query("initial_grid_file", m_initial_grid_file);