use runtime specnames for manifold EOS, also add more physics checks

CMake/BuildPelePhysicsLib.cmake

@@ -113,7 +113,6 @@ function(build_pele_physics_lib pele_physics_lib_name)
 
     target_sources(${pele_physics_lib_name}
       PRIVATE
-        ${PELE_PHYSICS_EOS_DIR}/EOS.cpp
         ${PELE_PHYSICS_EOS_DIR}/EOS.H
         ${PELE_PHYSICS_EOS_DIR}/GammaLaw.H
         ${PELE_PHYSICS_EOS_DIR}/Fuego.H

Exec/RegTests/FlameSheet/input.manifold

@@ -51,7 +51,12 @@ 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
+amr.derive_plot_vars = avg_pressure mag_vort mass_fractions maniout progress_variable
+
+peleLM.progressVariable.format = Cantera
+peleLM.progressVariable.weights = "PROG:1"
+peleLM.progressVariable.coldState = "PROG:0"
+peleLM.progressVariable.hotState = "PROG:1"
 
 # ------------------- INPUTS DERIVED DIAGS ------------------
 peleLM.fuel_name = CH4
@@ -82,4 +87,4 @@ manifold.density_lookup_type = log
 peleLM.use_wbar = 0
 peleLM.chi_correction_type = DivuFirstIter
 #peleLM.chi_correction_type = NoDivu
-peleLM.print_chi_convergence = 1
+peleLM.print_chi_convergence = 1

Source/PeleLMeX_Plot.cpp

@@ -175,7 +175,8 @@ PeleLM::WritePlotFile()
   //----------------------------------------------------------------
   // Components names
   Vector<std::string> names;
-  pele::physics::eos::speciesNames<pele::physics::PhysicsType::eos_type>(names);
+  pele::physics::eos::speciesNames<pele::physics::PhysicsType::eos_type>(
+    names, &(eos_parms.host_parm()));
 
   Vector<std::string> plt_VarsName;
   AMREX_D_TERM(plt_VarsName.push_back("x_velocity");
@@ -841,7 +842,7 @@ PeleLM::initLevelDataFromPlt(int a_lev, const std::string& a_dataPltFile)
   // Find required data in pltfile
   Vector<std::string> spec_names;
   pele::physics::eos::speciesNames<pele::physics::PhysicsType::eos_type>(
-    spec_names);
+    spec_names, &(eos_parms.host_parm()));
   int idT = -1, idV = -1, idY = -1, nSpecPlt = 0;
 #ifdef PELE_USE_EFIELD
   int inE = -1, iPhiV = -1;

Source/PeleLMeX_Setup.cpp

@@ -75,14 +75,17 @@ PeleLM::Setup()
   makeEBGeometry();
 #endif
 
-  // Setup the state variables
-  variablesSetup();
-
-  // Initialize EOS and others
+  // Initialize EOS
   if (m_incompressible == 0) {
     amrex::Print() << " Initialization of Eos ... \n";
     eos_parms.initialize();
+  }
+
+  // Setup the state variables
+  variablesSetup();
 
+  // Initialize Transport and others
+  if (m_incompressible == 0) {
     amrex::Print() << " Initialization of Transport ... \n";
 #ifdef USE_MANIFOLD_TRANSPORT
     trans_parms.host_only_parm().manfunc_par =
@@ -130,7 +133,8 @@ PeleLM::Setup()
       m_reactor =
         pele::physics::reactions::ReactorBase::create(m_chem_integrator);
       m_reactor->init(reactor_type, ncells_chem);
-      m_reactor->set_eos_parm(eos_parms.device_parm());
+      m_reactor->set_eos_parm(
+        &(eos_parms.host_parm()), eos_parms.device_parm());
       // For ReactorNull, we need to also skip instantaneous RR used in divU
       if (m_chem_integrator == "ReactorNull") {
         m_skipInstantRR = 1;
@@ -179,10 +183,8 @@ PeleLM::Setup()
   resizeArray();
 
   // Mixture fraction & Progress variable
-  if (pele::physics::PhysicsType::eos_type::identifier() != "Manifold") {
-    initMixtureFraction();
-    initProgressVariable();
-  }
+  initMixtureFraction();
+  initProgressVariable();
 
   // Initialize turbulence injection
   turb_inflow.init(Geom(0));
@@ -249,12 +251,6 @@ PeleLM::readParameters()
   if ((verbose != 0) && (m_closed_chamber != 0)) {
     Print() << " Simulation performed with the closed chamber algorithm \n";
   }
-  if (
-    (m_closed_chamber != 0) &&
-    (pele::physics::PhysicsType::eos_type::identifier() == "Manifold")) {
-    amrex::Abort(
-      "Simulation with closed chamber not supported for Manifold EOS");
-  }
 
 #ifdef PELE_USE_EFIELD
   ParmParse ppef("ef");
@@ -490,11 +486,6 @@ PeleLM::readParameters()
                       "fixed_Pr or fixed_Le is true"
                    << std::endl;
   }
-  if (
-    (m_use_wbar != 0) &&
-    (pele::physics::PhysicsType::eos_type::identifier() == "Manifold")) {
-    amrex::Abort("Use of Wbar fluxes is not compatible with Manifold EOS");
-  }
 
   pp.query("deltaT_verbose", m_deltaT_verbose);
   pp.query("deltaT_iterMax", m_deltaTIterMax);
@@ -757,6 +748,33 @@ PeleLM::readParameters()
   }
 #endif
 
+  // Ensure unsupported physics is not used with manifiold models
+  if (pele::physics::PhysicsType::eos_type::identifier() == "Manifold") {
+    if (m_closed_chamber != 0) {
+      amrex::Abort(
+        "Simulation with closed chamber is not yet supported for Manifold EOS");
+    }
+    if (m_use_wbar != 0) {
+      amrex::Abort("Use of Wbar fluxes is not compatible with Manifold EOS");
+    }
+#ifdef PELE_USE_RADIATION
+    if (do_rad_solve) {
+      amrex::Abort("Radiation models are not yet supported for Manifold EOS");
+    }
+#endif
+#ifdef PELE_USE_SOOT
+    if (do_soot_solve) {
+      amrex::Abort("Soot models are not yet supported for Manifold EOS");
+    }
+#endif
+#ifdef PELE_USE_SPRAY
+    amrex::Abort("Spray models are not yet supported for Manifold EOS");
+#endif
+#ifdef PELE_USE_EFIELD
+    amrex::Abort("Efield models are not yet supported for Manifold EOS");
+#endif
+  }
+
   // -----------------------------------------
   // External Sources
   // -----------------------------------------
@@ -843,7 +861,7 @@ PeleLM::variablesSetup()
     Print() << " First species: " << FIRSTSPEC << "\n";
     Vector<std::string> names;
     pele::physics::eos::speciesNames<pele::physics::PhysicsType::eos_type>(
-      names);
+      names, &(eos_parms.host_parm()));
     for (int n = 0; n < NUM_SPECIES; n++) {
       stateComponents.emplace_back(FIRSTSPEC + n, "rho.Y(" + names[n] + ")");
       reactComponents.emplace_back(n, "I_R(" + names[n] + ")");
@@ -1024,7 +1042,7 @@ PeleLM::derivedSetup()
     // Get species names
     Vector<std::string> spec_names;
     pele::physics::eos::speciesNames<pele::physics::PhysicsType::eos_type>(
-      spec_names);
+      spec_names, &(eos_parms.host_parm()));
 
     // Set species mass fractions
     Vector<std::string> var_names_massfrac(NUM_SPECIES);
@@ -1219,7 +1237,7 @@ PeleLM::evaluateSetup()
   // Get species names
   Vector<std::string> spec_names;
   pele::physics::eos::speciesNames<pele::physics::PhysicsType::eos_type>(
-    spec_names);
+    spec_names, &(eos_parms.host_parm()));
 
   // divU
   evaluate_lst.add("divU", IndexType::TheCellType(), 1, the_same_box);

Source/PeleLMeX_Utils.cpp

@@ -1019,7 +1019,7 @@ PeleLM::initProgressVariable()
 {
   Vector<std::string> varNames;
   pele::physics::eos::speciesNames<pele::physics::PhysicsType::eos_type>(
-    varNames);
+    varNames, &(eos_parms.host_parm()));
   varNames.push_back("temp");
 
   ParmParse pp("peleLM");
@@ -1416,7 +1416,7 @@ PeleLM::setTypicalValues(const TimeStamp& a_time, int is_init)
       Print() << "\tH:        " << typical_values[RHOH] << '\n';
       Vector<std::string> spec_names;
       pele::physics::eos::speciesNames<pele::physics::PhysicsType::eos_type>(
-        spec_names);
+        spec_names, &(eos_parms.host_parm()));
       for (int n = 0; n < NUM_SPECIES; n++) {
         Print() << "\tY_" << spec_names[n]
                 << std::setw(
@@ -1722,10 +1722,56 @@ PeleLM::checkMemory(const std::string& a_message) const
 void
 PeleLM::initMixtureFraction()
 {
-  // Get default fuel and oxy tank composition: pure fuel vs air
+  // set up a few variables
+  auto const* leosparm = &eos_parms.host_parm();
+  auto eos = pele::physics::PhysicsType::eos(leosparm);
   Vector<std::string> specNames;
   pele::physics::eos::speciesNames<pele::physics::PhysicsType::eos_type>(
-    specNames);
+    specNames, leosparm);
+  ParmParse pp("peleLM");
+
+  // Do simpler things for some EOS
+  if (pele::physics::PhysicsType::eos_type::identifier() == "GammaLaw") {
+    // Do nothing - Bilger mixture fraction has no meaning here
+    // error will be raised if user tries to derive mixfrac because Zfu is
+    // negative
+    return;
+  }
+  if (pele::physics::PhysicsType::eos_type::identifier() == "Manifold") {
+    // Just take a mixture fraction if it is a manifold parameter
+    // otherwise do nothing (an error will later be raised if the user tries to
+    // derive it) also raise an error if the user requests it and it is not
+    // found
+    std::string mixfrac_name = "ZMIX";
+    bool requested_mixfrac =
+      pp.contains("mixtureFraction.manifoldParameterName");
+    pp.query("mixtureFraction.manifoldParameterName", mixfrac_name);
+    bool found = false;
+    for (int n = 0; n < NUM_SPECIES; ++n) {
+      if (specNames[n] == mixfrac_name) {
+        if (!found) {
+          found = true;
+          spec_Bilger_fact[n] = 1.0;
+        } else {
+          amrex::Abort("initMixtureFraction: requested manifold parameter "
+                       "found multiple times");
+        }
+      } else {
+        spec_Bilger_fact[n] = 0.0;
+      }
+    }
+    if (found) {
+      Zfu = 1.0;
+      Zox = 0.0;
+    } else if (requested_mixfrac) {
+      amrex::Abort(
+        "initMixtureFraction: requested manifold parameter not found");
+    }
+    return;
+  }
+
+  // Otherwise - compute Bilger weights for detailed chemistry
+
   amrex::Real YF[NUM_SPECIES], YO[NUM_SPECIES];
   for (int i = 0; i < NUM_SPECIES; ++i) {
     YF[i] = 0.0;
@@ -1741,10 +1787,7 @@ PeleLM::initMixtureFraction()
     }
   }
 
-  auto const* leosparm = &eos_parms.host_parm();
-  auto eos = pele::physics::PhysicsType::eos(leosparm);
   // Overwrite with user-defined value if provided in input file
-  ParmParse pp("peleLM");
   std::string MFformat;
   int hasUserMF = static_cast<int>(pp.contains("mixtureFraction.format"));
   if (hasUserMF != 0) {
@@ -1796,7 +1839,6 @@ PeleLM::initMixtureFraction()
         for (int i = 0; i < NUM_SPECIES; ++i) {
           XF[i] = compositionIn[i];
         }
-
         eos.X2Y(XO, YO);
         eos.X2Y(XF, YF);
       } else {
@@ -1812,6 +1854,7 @@ PeleLM::initMixtureFraction()
                "peleLM.fuel_name keyword \n";
   }
 
+  // Detailed chem - compute Bilger coefficients
   // Only interested in CHON -in that order. Compute Bilger weights
   amrex::Real atwCHON[4] = {0.0};
   pele::physics::eos::atomic_weightsCHON<pele::physics::PhysicsType::eos_type>(
@@ -1852,7 +1895,7 @@ PeleLM::parseComposition(
   // Get species names
   Vector<std::string> specNames;
   pele::physics::eos::speciesNames<pele::physics::PhysicsType::eos_type>(
-    specNames);
+    specNames, &(eos_parms.host_parm()));
 
   // For each entry in the user-provided composition, parse name and value
   std::string delimiter = ":";