User function to patch flow variables after restart from plot file (#398

) * Added function declaration for patching ignition kernel * Corrected for running on HIP * Made changes to EB Backward facing step flame case * Updated Pelephysics * Clang Tidyied * Compiled and check EB_BFSF * Added documentation * Modified ignition patching to flow variable patching * Formatted * Corrected input file for EBBFS * cleaned up some of the print statements * Changed input file EBBF to be reactive * Incorporating Bruce's suggestions * Making changes suggested by Bruce * clang tidy --------- Co-authored-by: Bruce Perry <[email protected]>
AMReX-Combustion · Jul 20, 2024 · e79f505 · e79f505
1 parent 5c77f2e
commit e79f505
Show file tree

Hide file tree

Showing 16 changed files with 86 additions and 6 deletions.
diff --git a/CMake/BuildPeleExe.cmake b/CMake/BuildPeleExe.cmake
@@ -46,6 +46,8 @@ function(build_pele_exe pele_exe_name pele_physics_lib_name)
        ${SRC_DIR}/PeleLMeX_Evolve.cpp
        ${SRC_DIR}/PeleLMeX_FlowController.cpp
        ${SRC_DIR}/PeleLMeX_Forces.cpp
+       ${SRC_DIR}/PeleLMeX_PatchFlowVariables.H
+       ${SRC_DIR}/PeleLMeX_PatchFlowVariables.cpp
        ${SRC_DIR}/PeleLMeX_Init.cpp
        ${SRC_DIR}/PeleLMeX_Plot.cpp
        ${SRC_DIR}/PeleLMeX_Projection.cpp

diff --git a/Docs/sphinx/manual/LMeXControls.rst b/Docs/sphinx/manual/LMeXControls.rst
@@ -255,6 +255,7 @@ PeleLMeX algorithm
     peleLM.deltaT_iterMax = 5              # [OPT, DEF=10] Maximum number of deltaT iterations
     peleLM.deltaT_tol = 1e-10              # [OPT, DEF=1.e-10] Tolerance of the deltaT solve
     peleLM.evaluate_vars =...              # [OPT, DEF=""] In evaluate mode, list unitTest: diffTerm, divU, instRR, transportCC
+    peleLM.do_patch_flow_variables = false # [OPT, DEF=false] Enable user-defined flow variable patching after reading a plot solution file
 
 Transport coefficients and LES
 ------------------------------

diff --git a/Docs/sphinx/manual/Tutorials_BFSFlame.rst b/Docs/sphinx/manual/Tutorials_BFSFlame.rst
@@ -129,6 +129,7 @@ The user parameters are gathered in the struct defined in  ``pelelmex_prob_parm.
         amrex::Real T_hot  = 1800.0_rt;
         amrex::Real Twall  = 300.0_rt;
         amrex::Real meanFlowMag = 0.0;
+
     };
 
 * ``T_mean`` : inlet and initial gas temperature
@@ -145,10 +146,12 @@ The user parameters are gathered in the struct defined in  ``pelelmex_prob_parm.
 
 * ``meanFlowMag`` : inlet :math:`x` velocity
 
+
 The initial solution consists of a premixed methane/air mixture in the upper part of the domain
 and pure hot air in the wake of the step. The default parameters provided above are overwritten
 using AMReX ParmParse in ``pelelmex_prob.cpp`` and the initial/boundary conditions implemented in
-``pelelmex_prob.H``.
+``pelelmex_prob.H``. Alternatively, the user can write a custom function to enforce an ignition kernel through the ``patchFlowVariables`` function in the problem-specific ``PeleLMeX_PatchFlowVariables.cpp`` file. 
+It should be kept in mind that the ``patchFlowVariables`` function can be used if the user wants to patch certain flow variables after reading an existing solution from a plot file ( ``peleLM.do_patch_flow_variables`` should be set to true).
 
 In addition to these three C++ files, an extra header is needed in the present case compared to
 :doc:`Tutorials_FlameSheet` : ``EBUserDefined.H``. This file is necessary to specify more complex EB

diff --git a/Exec/RegTests/EB_BackwardStepFlame/GNUmakefile b/Exec/RegTests/EB_BackwardStepFlame/GNUmakefile
@@ -27,6 +27,7 @@ THREAD_SANITIZER = FALSE
 
 # PeleLMeX
 CEXE_headers += EBUserDefined.H
+CEXE_sources += PeleLMeX_PatchFlowVariables.cpp
 
 # PelePhysics
 Chemistry_Model = drm19

diff --git a/Exec/RegTests/EB_BackwardStepFlame/PeleLMeX_PatchFlowVariables.cpp b/Exec/RegTests/EB_BackwardStepFlame/PeleLMeX_PatchFlowVariables.cpp
@@ -0,0 +1,31 @@
+
+#include "PeleLMeX.H"
+using namespace amrex;
+
+void
+patchFlowVariables(
+  const amrex::Geometry& geom, ProbParm const& lprobparm, amrex::MultiFab& a_mf)
+{
+
+  amrex::Print() << "\nPatching flow variables..";
+  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> prob_lo =
+    geom.ProbLoArray();
+  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> dx = geom.CellSizeArray();
+
+  for (amrex::MFIter mfi(a_mf, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
+    const amrex::Box& bx = mfi.tilebox();
+    auto const& temp_arr = a_mf.array(mfi, TEMP);
+    amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+      amrex::Real x[3] = {
+        prob_lo[0] + static_cast<amrex::Real>(i + 0.5) * dx[0],
+        prob_lo[1] + static_cast<amrex::Real>(j + 0.5) * dx[1],
+        prob_lo[2] + static_cast<amrex::Real>(k + 0.5) * dx[2]};
+
+      amrex::ignore_unused(x);
+      /*User can define how to patch flow variables here.*/
+      temp_arr(i, j, k) = lprobparm.T_mean;
+    });
+  }
+
+  amrex::Print() << "Done\n";
+}
diff --git a/Exec/RegTests/EB_BackwardStepFlame/input.2d b/Exec/RegTests/EB_BackwardStepFlame/input.2d
@@ -29,14 +29,13 @@ prob.T_wall = 300
 prob.meanFlowMag = 10.0
 prob.Y_fuel = 0.0445269
 prob.Y_oxid = 0.2226345
-
 #---------------------- PeleLM CONTROL ---------------------------
 peleLM.v = 1                                # PeleLMeX verbose
 peleLM.use_wbar = 1                         # Include Wbar term in species diffusion fluxes
 peleLM.sdc_iterMax = 2                      # Number of SDC iterations
 peleLM.num_init_iter = 3                    # Number of initial iterations
 peleLM.floor_species = 0                    # Enforce species positivity (non-conservative !)
-
+peleLM.do_patch_flow_variables = false      # Enable flow variable patching after reading a plot solution file
 #---------------------- Temporal CONTROL -------------------------
 peleLM.do_temporals = 1                     # Turn temporals ON/OFF
 peleLM.temporal_int = 10                    # Frequency of temporals
@@ -52,6 +51,7 @@ amr.dt_shrink = 0.001                       # Scale back initial timestep
 amr.dt_change_max = 1.1                     # Maximum dt increase btw successive steps
 
 #---------------------- IO CONTROL -------------------------------
+#amr.initDataPlt = plt_cold
 #amr.restart = chk00250                     # Restart checkpoint file
 amr.check_int = 100                         # Frequency of checkpoint output
 amr.plot_int = 50                           # Frequency of pltfile output

diff --git a/Exec/RegTests/EB_BackwardStepFlame/pelelmex_prob.cpp b/Exec/RegTests/EB_BackwardStepFlame/pelelmex_prob.cpp
@@ -13,4 +13,9 @@ PeleLM::readProbParm() // NOLINT(readability-make-member-function-const)
   pp.query("T_hot", prob_parm->T_hot);
   pp.query("T_wall", prob_parm->Twall);
   pp.query("meanFlowMag", prob_parm->meanFlowMag);
+  pp.query("do_ignite", prob_parm->ignite_flow);
+
+  if (prob_parm->ignite_flow) {
+    pp.query("ignition_temperature", prob_parm->ignition_temperature);
+  }
 }
diff --git a/Exec/RegTests/EB_BackwardStepFlame/pelelmex_prob_parm.H b/Exec/RegTests/EB_BackwardStepFlame/pelelmex_prob_parm.H
@@ -15,5 +15,8 @@ struct ProbParm
   amrex::Real T_hot = 1800.0_rt;
   amrex::Real Twall = 300.0_rt;
   amrex::Real meanFlowMag = 0.0;
+
+  bool ignite_flow = false;
+  amrex::Real ignition_temperature = 1700.0;
 };
 #endif
diff --git a/Source/Make.package b/Source/Make.package
@@ -9,6 +9,7 @@ CEXE_headers += PeleLMeX_DeriveFunc.H
 CEXE_headers += PeleLMeX_EBUserDefined.H
 CEXE_headers += PeleLMeX_FlowControllerData.H
 CEXE_headers += PeleLMeX_BPatch.H
+CEXE_headers += PeleLMeX_PatchFlowVariables.H
 
 ## Sources
 CEXE_sources += main.cpp
@@ -43,6 +44,7 @@ CEXE_sources += PeleLMeX_Diagnostics.cpp
 CEXE_sources += PeleLMeX_FlowController.cpp
 CEXE_sources += PeleLMeX_DeriveUserDefined.cpp
 CEXE_sources += PeleLMeX_BPatch.cpp
+CEXE_sources += PeleLMeX_PatchFlowVariables.cpp
 
 ifeq ($(USE_SOOT), TRUE)
   CEXE_sources += PeleLMeX_Soot.cpp

diff --git a/Source/PeleLMeX.H b/Source/PeleLMeX.H
@@ -1787,6 +1787,7 @@ public:
   int m_ioDigits = 5;
   amrex::Vector<std::string> m_evaluatePlotVars;
   bool m_write_hdf5_pltfile = false;
+  bool m_do_patch_flow_variables = false;
 
   //-----------------------------------------------------------------------------
   // ALGORITHM

diff --git a/Source/PeleLMeX_PatchFlowVariables.H b/Source/PeleLMeX_PatchFlowVariables.H
@@ -0,0 +1,10 @@
+#ifndef PELELM_PATCHFLOWVARIABLES_H
+#define PELELM_PATCHFLOWVARIABLES_H
+
+class PeleLM;
+
+void patchFlowVariables(
+  const amrex::Geometry& geom,
+  ProbParm const& prob_parm,
+  amrex::MultiFab& a_mf);
+#endif
diff --git a/Source/PeleLMeX_PatchFlowVariables.cpp b/Source/PeleLMeX_PatchFlowVariables.cpp
@@ -0,0 +1,15 @@
+#include <PeleLMeX.H>
+#include <PeleLMeX_K.H>
+#include <PeleLMeX_PatchFlowVariables.H>
+
+using namespace amrex;
+
+void
+patchFlowVariables(
+  const amrex::Geometry& /*geom*/,
+  ProbParm const& /*prob_parm*/,
+  amrex::MultiFab& /*a_mf*/)
+{
+  Abort("Using patchFlowVariables requires providing a definition in local "
+        "PatchFlowVariables.cpp");
+}
diff --git a/Source/PeleLMeX_Plot.cpp b/Source/PeleLMeX_Plot.cpp
@@ -6,6 +6,7 @@
 #include <AMReX_ParmParse.H>
 #include <PeleLMeX_BCfill.H>
 #include <AMReX_FillPatchUtil.H>
+#include <PeleLMeX_PatchFlowVariables.H>
 #include <memory>
 #ifdef AMREX_USE_EB
 #include <AMReX_EBInterpolater.H>
@@ -951,6 +952,12 @@ PeleLM::initLevelDataFromPlt(int a_lev, const std::string& a_dataPltFile)
 
   ProbParm const* lprobparm = prob_parm_d;
 
+  // If m_do_patch_flow_variables is set as true, call user-defined function to
+  // patch flow variables
+  if (m_do_patch_flow_variables) {
+    patchFlowVariables(geom[a_lev], *lprobparm, ldata_p->state);
+  }
+
   // Enforce rho and rhoH consistent with temperature and mixture
   // The above handles species mapping (to some extent), but nothing enforce
   // sum of Ys = 1 -> use N2 in the following if N2 is present

diff --git a/Source/PeleLMeX_Setup.cpp b/Source/PeleLMeX_Setup.cpp
@@ -45,7 +45,6 @@ PeleLM::Setup()
   // Ensure grid is isotropic
   {
     auto const dx = geom[0].CellSizeArray();
-    amrex::Print() << "\n Dx = " << dx[0] << " " << dx[1] << " " << dx[2];
     AMREX_ALWAYS_ASSERT(AMREX_D_TERM(
       , amrex::almostEqual(dx[0], dx[1], 10),
       &&amrex::almostEqual(dx[1], dx[2], 10)));
@@ -419,6 +418,7 @@ PeleLM::readParameters()
   pp.query("num_divu_iter", m_numDivuIter);
   pp.query("do_init_proj", m_do_init_proj);
   pp.query("num_init_iter", m_init_iter);
+  pp.query("do_patch_flow_variables", m_do_patch_flow_variables);
 
   // -----------------------------------------
   // advance

diff --git a/Source/PeleLMeX_Temporals.cpp b/Source/PeleLMeX_Temporals.cpp
@@ -622,7 +622,6 @@ PeleLM::addRhoYFluxesPatch(
         ParallelAllReduce::Sum<Real>(
           {sum_species_flux_global}, ParallelContext::CommunicatorSub());
         bphost->speciesFlux[m] = a_factor * sum_species_flux_global;
-        // amrex::Print()<<"\nNew func = "<<a_factor * sum_species_flux_global;
       }
     }
   }

diff --git a/Submodules/PelePhysics b/Submodules/PelePhysics