Add full manifold model capability

Exec/Make.PeleLMeX

@@ -77,6 +77,19 @@ endif
 ifeq ($(Eos_Model),$(filter $(Eos_Model),Soave-Redlich-Kwong))
    DEFINES += -DUSE_SRK_EOS
 endif
+ifeq ($(Eos_Model),$(filter $(Eos_Model),Manifold))
+   DEFINES += -DUSE_MANIFOLD_EOS
+   DEFINES += -DMANIFOLD_DIM=$(Manifold_Dim)
+   ifeq ($(Manifold_Type),Table)
+      DEFINES += -DMANIFOLD_EOS_TYPE=1
+   else
+      ifeq ($(Manifold_Type),Network)
+         DEFINES += -DMANIFOLD_EOS_TYPE=2
+      else
+         DEFINES += -DMANIFOLD_EOS_TYPE=0
+      endif
+   endif
+endif
 
 # Transport model switches
 ifeq ($(Transport_Model), Simple)
@@ -92,6 +105,9 @@ endif
 ifeq ($(Transport_Model), Sutherland)
   DEFINES += -DUSE_SUTHERLAND_TRANSPORT
 endif
+ifeq ($(Transport_Model), Manifold)
+   DEFINES += -DUSE_MANIFOLD_TRANSPORT
+endif
 
 ifeq ($(PELE_USE_KLU), TRUE)
   DEFINES += -DPELE_USE_KLU

Exec/RegTests/FlameSheet/GNUmakefile

@@ -29,9 +29,18 @@ THREAD_SANITIZER = FALSE
 USE_EFIELD = FALSE
 
 # PelePhysics
-Chemistry_Model = drm19
-Eos_Model = Fuego
-Transport_Model = Simple
+USE_MANIFOLD = FALSE
+ifeq ($(USE_MANIFOLD), TRUE)
+  Eos_Model       = Manifold
+  Chemistry_Model = Null
+  Transport_Model = Manifold
+  Manifold_Dim = 1
+ # Manifold_Type = Table
+else
+  Chemistry_Model = drm19
+  Eos_Model = Fuego
+  Transport_Model = Simple
+endif
 
 PELE_HOME ?= ../../..
 include $(PELE_HOME)/Exec/Make.PeleLMeX

Exec/RegTests/FlameSheet/input.manifold

@@ -0,0 +1,85 @@
+#----------------------DOMAIN DEFINITION------------------------
+geometry.is_periodic = 0 0               # For each dir, 0: non-perio, 1: periodic
+geometry.coord_sys   = 0                  # 0 => cart, 1 => RZ
+geometry.prob_lo     = 0.0 0.0 0.0        # x_lo y_lo (z_lo)
+geometry.prob_hi     = 0.003 0.012 0.016        # x_hi y_hi (z_hi)
+
+# >>>>>>>>>>>>>  BC FLAGS <<<<<<<<<<<<<<<<
+# Interior, Inflow, Outflow, Symmetry,
+# SlipWallAdiab, NoSlipWallAdiab, SlipWallIsotherm, NoSlipWallIsotherm
+peleLM.lo_bc = Interior Inflow
+peleLM.hi_bc = Interior Outflow
+peleLM.lo_bc =  SlipWallAdiab Inflow
+peleLM.hi_bc =  SlipWallAdiab Outflow
+
+
+#-------------------------AMR CONTROL----------------------------
+amr.n_cell          = 64 256 32      # Level 0 number of cells in each direction
+amr.v               = 1                # AMR verbose
+amr.max_level       = 0                # maximum level number allowed
+amr.ref_ratio       = 2 2 2 2          # refinement ratio
+amr.regrid_int      = 5                # how often to regrid
+amr.n_error_buf     = 1 1 2 2          # number of buffer cells in error est
+amr.grid_eff        = 0.7              # what constitutes an efficient grid
+amr.blocking_factor = 16               # block factor in grid generation (min box size)
+amr.max_grid_size   = 64               # max box size
+
+
+#--------------------------- Problem -------------------------------
+prob.P_mean = 101325.0
+prob.standoff = -.03
+prob.pertmag = 0.0000
+pmf.datafile = "prem_drm19_phi1_p1_t298_mani.dat"
+
+#-------------------------PeleLM CONTROL----------------------------
+peleLM.v = 3
+peleLM.incompressible = 0
+peleLM.rho = 1.17
+peleLM.mu = 0.0
+peleLM.sdc_iterMax = 1
+peleLM.floor_species = 0
+
+peleLM.do_temporals = 1
+peleLM.temporal_int = 2
+peleLM.mass_balance = 1
+
+#amr.restart = chk00005
+#amr.check_int = 20
+amr.plot_int = 1
+amr.max_step = 10
+amr.dt_shrink = 0.01
+amr.stop_time = 0.01
+#amr.stop_time = 1.00
+amr.cfl = 0.5
+amr.derive_plot_vars = avg_pressure mag_vort mass_fractions maniout
+
+# ------------------- INPUTS DERIVED DIAGS ------------------
+peleLM.fuel_name = CH4
+
+# --------------- INPUTS TO CHEMISTRY REACTOR ---------------
+peleLM.chem_integrator = "ReactorRK64"
+
+mac_proj.verbose = 0
+nodal_proj.verbose = 0
+
+#--------------------REFINEMENT CONTROL------------------------
+amr.refinement_indicators = temp
+amr.temp.max_level     = 3
+amr.temp.adjacent_difference_greater = 10
+amr.temp.field_name    = temp
+
+#amrex.fpe_trap_invalid = 1
+#amrex.fpe_trap_zero = 1
+#amrex.fpe_trap_overflow = 1
+
+#----------------------- Manifold Stuff ----------------------------
+manifold.model = Table
+manifold.table.v = 1
+manifold.table.filename = prem_drm19_phi1_p1_t298.ctb
+manifold.nominal_pressure_cgs = 1013250.0
+manifold.compute_temperature = 1
+manifold.density_lookup_type = log
+peleLM.use_wbar = 0
+peleLM.chi_correction_type = DivuFirstIter
+#peleLM.chi_correction_type = NoDivu
+peleLM.print_chi_convergence = 1

Exec/RegTests/FlameSheet/pelelmex_prob.H

@@ -45,7 +45,6 @@ pelelmex_initdata(
   amrex::Real pert = 0.0;
 
   if (prob_parm.pertmag > 0.0) {
-
 #if (AMREX_SPACEDIM == 2)
     amrex::GpuArray<amrex::Real, AMREX_SPACEDIM - 1> lpert{Lx};
     if (prob_parm.pertlength > 0.0) {
@@ -58,10 +57,8 @@ pelelmex_initdata(
             1.017 * std::sin(2 * Pi * 5 * (x - 0.0033) / lpert[0]) +
             0.982 * std::sin(2 * Pi * 5 * (x - 0.014234) / lpert[0]));
   }
-
   amrex::Real y1 = (y - prob_parm.standoff - 0.5 * dx[1] + pert) * 100;
   amrex::Real y2 = (y - prob_parm.standoff + 0.5 * dx[1] + pert) * 100;
-
 #elif (AMREX_SPACEDIM == 3)
     amrex::GpuArray<amrex::Real, AMREX_SPACEDIM - 1> lpert{Lx, Ly};
     if (prob_parm.pertlength > 0.0) {
@@ -80,29 +77,36 @@ pelelmex_initdata(
               std::sin(2 * Pi * 6 * (y - 0.018) / lpert[1]) +
             0.982 * std::sin(2 * Pi * 5 * (x - 0.014234) / lpert[0]));
   }
-
   amrex::Real y1 = (z - prob_parm.standoff - 0.5 * dx[2] + pert) * 100;
   amrex::Real y2 = (z - prob_parm.standoff + 0.5 * dx[2] + pert) * 100;
 #endif
 
   pele::physics::PMF::pmf(pmf_data, y1, y2, pmf_vals);
-
   state(i, j, k, TEMP) = pmf_vals[0];
-  ;
+  for (int dim = 0; dim < AMREX_SPACEDIM - 1; dim++) {
+    state(i, j, k, VELX + dim) = 0.0;
+  }
+  state(i, j, k, VELX + AMREX_SPACEDIM - 1) = pmf_vals[1] * 1e-2;
 
+#ifdef USE_MANIFOLD_EOS
+  for (int n = 0; n < NUM_SPECIES - 1; n++) {
+    massfrac[n] = pmf_vals[3 + n];
+  }
+  amrex::Real rho_temp, tmp1, tmp2;
+  eos.PYT2R(tmp1, massfrac, tmp2, rho_temp);
+  rho_temp *= 1.0e3;
+  state(i, j, k, DENSITY) = rho_temp; // CGS -> MKS conversion
+  state(i, j, k, RHOH) = 0.0;         // No RhoH transport in manifold model
+  for (int n = 0; n < NUM_SPECIES - 1; n++) {
+    state(i, j, k, FIRSTSPEC + n) = massfrac[n] * rho_temp;
+  }
+  state(i, j, k, FIRSTSPEC + NUM_SPECIES - 1) = rho_temp;
+#else
   for (int n = 0; n < NUM_SPECIES; n++) {
     molefrac[n] = pmf_vals[3 + n];
   }
   eos.X2Y(molefrac, massfrac);
 
-  state(i, j, k, VELX) = 0.0;
-#if (AMREX_SPACEDIM == 2)
-  state(i, j, k, VELY) = pmf_vals[1] * 1e-2;
-#elif (AMREX_SPACEDIM == 3)
-  state(i, j, k, VELY) = 0.0;
-  state(i, j, k, VELZ) = pmf_vals[1] * 1e-2;
-#endif
-
   amrex::Real P_cgs = prob_parm.P_mean * 10.0;
 
   // Density
@@ -119,6 +123,7 @@ pelelmex_initdata(
   for (int n = 0; n < NUM_SPECIES; n++) {
     state(i, j, k, FIRSTSPEC + n) = massfrac[n] * state(i, j, k, DENSITY);
   }
+#endif
 }
 
 AMREX_GPU_DEVICE
@@ -144,15 +149,25 @@ bcnormal(
   if (sgn == 1) {
     pele::physics::PMF::pmf(pmf_data, prob_lo[idir], prob_lo[idir], pmf_vals);
 
-    s_ext[VELX] = 0.0;
-#if (AMREX_SPACEDIM == 2)
-    s_ext[VELY] = pmf_vals[1] * 1e-2;
-#elif (AMREX_SPACEDIM == 3)
-    s_ext[VELY] = 0.0;
-    s_ext[VELZ] = pmf_vals[1] * 1e-2;
-#endif
-
+    for (int dim = 0; dim < AMREX_SPACEDIM - 1; dim++) {
+      s_ext[VELX + dim] = 0.0;
+    }
+    s_ext[VELX + AMREX_SPACEDIM - 1] = pmf_vals[1] * 1e-2;
     s_ext[TEMP] = pmf_vals[0];
+
+#ifdef USE_MANIFOLD_EOS
+    for (int n = 0; n < NUM_SPECIES; n++) {
+      massfrac[n] = pmf_vals[3 + n];
+    }
+    amrex::Real P_dummy, rho_cgs;
+    eos.PYT2R(P_dummy, massfrac, s_ext[TEMP], rho_cgs);
+    s_ext[DENSITY] = rho_cgs * 1.0e3;
+    s_ext[RHOH] = 0.0;
+    for (int n = 0; n < NUM_SPECIES - 1; n++) {
+      s_ext[FIRSTSPEC + n] = massfrac[n] * s_ext[DENSITY];
+    }
+    s_ext[FIRSTSPEC + NUM_SPECIES - 1] = s_ext[DENSITY];
+#else
     for (int n = 0; n < NUM_SPECIES; n++) {
       molefrac[n] = pmf_vals[3 + n];
     }
@@ -170,6 +185,7 @@ bcnormal(
     for (int n = 0; n < NUM_SPECIES; n++) {
       s_ext[FIRSTSPEC + n] = massfrac[n] * s_ext[DENSITY];
     }
+#endif
   }
 }
 

Source/PeleLMeX.H

@@ -606,6 +606,7 @@ public:
    * velocity, ensuring they sum up to zero. \param a_fluxes diffusion fluxes to
    * be updated \param a_spec species rhoYs state data on all levels
    */
+  template <typename EOSType>
   void adjustSpeciesFluxes(
     const amrex::Vector<amrex::Array<amrex::MultiFab*, AMREX_SPACEDIM>>&
       a_spfluxes,

Source/PeleLMeX_Advection.cpp

@@ -488,9 +488,8 @@ PeleLM::computeScalarAdvTerms(std::unique_ptr<AdvanceAdvData>& advData)
                   afrac] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
               rho_ed(i, j, k) = 0.0;
               if (afrac(i, j, k) > 0.0) { // Uncovered faces
-                for (int n = 0; n < NUM_SPECIES; n++) {
-                  rho_ed(i, j, k) += rhoY_ed(i, j, k, n);
-                }
+                pele::physics::PhysicsType::eos_type::RY2R(
+                  rhoY_ed.cellData(i, j, k), rho_ed(i, j, k));
               }
             });
         } else // Regular boxes
@@ -499,10 +498,8 @@ PeleLM::computeScalarAdvTerms(std::unique_ptr<AdvanceAdvData>& advData)
           amrex::ParallelFor(
             ebx,
             [rho_ed, rhoY_ed] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
-              rho_ed(i, j, k) = 0.0;
-              for (int n = 0; n < NUM_SPECIES; n++) {
-                rho_ed(i, j, k) += rhoY_ed(i, j, k, n);
-              }
+              pele::physics::PhysicsType::eos_type::RY2R(
+                rhoY_ed.cellData(i, j, k), rho_ed(i, j, k));
             });
         }
       }
@@ -751,11 +748,9 @@ PeleLM::computeScalarAdvTerms(std::unique_ptr<AdvanceAdvData>& advData)
     amrex::ParallelFor(
       advData->AofS[lev],
       [=] AMREX_GPU_DEVICE(int box_no, int i, int j, int k) noexcept {
-        aofsma[box_no](i, j, k, DENSITY) = 0.0;
-        for (int n = 0; n < NUM_SPECIES; n++) {
-          aofsma[box_no](i, j, k, DENSITY) +=
-            aofsma[box_no](i, j, k, FIRSTSPEC + n);
-        }
+        pele::physics::PhysicsType::eos_type::RY2R(
+          aofsma[box_no].cellData(i, j, k), aofsma[box_no](i, j, k, DENSITY),
+          FIRSTSPEC);
       });
   }
   Gpu::streamSynchronize();

Source/PeleLMeX_DeriveFunc.H

@@ -290,6 +290,22 @@ void pelelmex_deruserdef(
   const amrex::Vector<amrex::BCRec>& bcrec,
   int level);
 
+#ifdef USE_MANIFOLD_EOS
+void pelelmex_dermaniout(
+  PeleLM* a_pelelm,
+  const amrex::Box& bx,
+  amrex::FArrayBox& derfab,
+  int dcomp,
+  int ncomp,
+  const amrex::FArrayBox& statefab,
+  const amrex::FArrayBox& reactfab,
+  const amrex::FArrayBox& pressfab,
+  const amrex::Geometry& geom,
+  amrex::Real time,
+  const amrex::Vector<amrex::BCRec>& bcrec,
+  int level);
+#endif
+
 #ifdef PELE_USE_EFIELD
 #include <PeleLMeX_EFDeriveFunc.H>
 #endif