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;
+  bool m_plotExtSource = false;
   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_Advance.cpp

@@ -105,6 +105,9 @@ PeleLM::Advance(int is_initIter)
     poissonSolveEF(AmrOldTime);
 #endif
   }
+  if (m_do_les) {
+    calcTurbViscosity(AmrOldTime);
+  }
 
   //----------------------------------------------------------------
   BL_PROFILE_VAR_STOP(PLM_SETUP);

Source/PeleLMeX_Evaluate.cpp

@@ -171,6 +171,10 @@ PeleLM::MLevaluate(
       finest_level, grids, dmap, m_factory, m_nGrowAdv, m_use_wbar,
       m_use_soret);
     calcDiffusivity(AmrNewTime);
+    // If doing LES, need to be able to add in turbulent component
+    if (m_do_les) {
+      calcTurbViscosity(AmrNewTime);
+    }
     computeDifferentialDiffusionTerms(AmrNewTime, diffData);
     for (int lev = 0; lev <= finest_level; ++lev) {
       MultiFab::Copy(
@@ -215,6 +219,7 @@ PeleLM::MLevaluate(
   } else if (a_var == "transportCC") {
     // Cell-centered transport coefficients functions go through the level
     // data container. Simply copy once the later has been filled.
+    // No LES component here - this is just moledular transport coeffs
     calcViscosity(AmrNewTime);
     calcDiffusivity(AmrNewTime);
     for (int lev = 0; lev <= finest_level; ++lev) {
@@ -290,6 +295,9 @@ PeleLM::evaluateChemExtForces(
   // compute t^{n} data
   calcViscosity(AmrOldTime);
   calcDiffusivity(AmrOldTime);
+  if (m_do_les) {
+    calcTurbViscosity(AmrOldTime);
+  }
 
   floorSpecies(AmrOldTime);
   setThermoPress(AmrOldTime);
@@ -392,6 +400,9 @@ PeleLM::evaluateAdvectionTerms(
   // compute t^{n} data
   calcViscosity(AmrOldTime);
   calcDiffusivity(AmrOldTime);
+  if (m_do_les) {
+    calcTurbViscosity(AmrOldTime);
+  }
 
   floorSpecies(AmrOldTime);
   setThermoPress(AmrOldTime);

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
+  amrex::ignore_unused(a_timestamp);
+#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 @@
     ncomp += 1;
   }
 
+  if (m_plotExtSource) {
+    ncomp += NVAR;
+  }
+
   //----------------------------------------------------------------
   // Plot MultiFabs
   Vector<MultiFab> mf_plt(finest_level + 1);
@@ -282,6 +286,12 @@
     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 @@
         });
       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);

Source/PeleLMeX_TransportProp.cpp

@@ -414,24 +414,9 @@ PeleLM::getDiffusivity(
   // supported
   if ((addTurbContrib != 0) and m_do_les) {
 
-    // If initializing the simulation, always recompute turbulent viscosity
-    // otherwise, only recompute once per level per timestep (at old time)
-    // calcTurbViscosity computes for all levels, so only call from the base
-    // level
-    TimeStamp tstamp;
-    if (getTime(lev, AmrNewTime) == 0.0) {
-      tstamp = AmrNewTime;
-      if (lev == 0) {
-        calcTurbViscosity(tstamp);
-      }
-    } else if (lev == 0 and getTime(lev, AmrOldTime) > m_turb_visc_time[lev]) {
-      tstamp = AmrOldTime;
-      calcTurbViscosity(tstamp);
-      m_turb_visc_time[lev] = getTime(lev, AmrOldTime);
-    } else {
-      tstamp = AmrOldTime;
-    }
-    auto* ldata_p = getLevelDataPtr(lev, tstamp);
+    // Turbulent viscosity - always in old data
+    // (not updated because velocity doesn't get updated until end of advance)
+    auto* ldata_p = getLevelDataPtr(lev, AmrOldTime);
 
     // Identify and add the correct turbulent contribution
     for (int idim = 0; idim < AMREX_SPACEDIM; idim++) {