Boundary patch implementation (AMReX-Combustion#369)

* Added class BPatch and associated functions to print species fluxes through user-defined patches based on simple geometric shapes

* Formatted using clang-format

* Removed declaration of an old and obsolete function specific to A74 case

* Modified gitignore file

* Compiled using GNU and previous errors resolved. Need to try running tests on Github.

* Clang formatted

* Reverted back to old clang-format

* Reverted Init file to clan format

* Corrected warnings. Corrected to clang-format 18

* Corrected 2 more warning on gcc compiler.

* Corrected clang tidy warnings

* correct host/device data access for bpatches

* add block for bpatches in RZ coordinates

* Added documentation for bpatch

* Clang Tidy ied Temporals.cpp

---------

Co-authored-by: Bruce Perry <[email protected]>

.gitignore

@@ -16,3 +16,10 @@ temporals
 # OSX files
 .DS_Store
 .AppleDouble
+
+#Eclipse file
+.cproject
+.project
+Exec/**/*.cproject
+Exec/**/*.project
+

CMake/BuildPeleExe.cmake

@@ -53,6 +53,8 @@ function(build_pele_exe pele_exe_name pele_physics_lib_name)
        ${SRC_DIR}/PeleLMeX_Regrid.cpp
        ${SRC_DIR}/PeleLMeX_Setup.cpp
        ${SRC_DIR}/PeleLMeX_Tagging.cpp
+       ${SRC_DIR}/PeleLMeX_BPatch.H
+       ${SRC_DIR}/PeleLMeX_BPatch.cpp
        ${SRC_DIR}/PeleLMeX_Temporals.cpp
        ${SRC_DIR}/PeleLMeX_Timestep.cpp
        ${SRC_DIR}/PeleLMeX_TransportProp.cpp

Docs/sphinx/manual/LMeXControls.rst

@@ -426,6 +426,44 @@ to activate `temporal` diagnostics performing these reductions at given interval
     peleLM.do_extremas = 1                      # [OPT, DEF=0] Trigger extremas, if temporals activated
     peleLM.do_mass_balance = 1                  # [OPT, DEF=0] Compute mass balance, if temporals activated
     peleLM.do_species_balance = 1               # [OPT, DEF=0] Compute species mass balance, if temporals activated
+    peleLM.do_patch_mfr=1                       # [OPT, DEF=0] Activate patch based species flux diagbostics
+    peleLM.bpatch.patchnames= <patch_name1 patch_name2 ..> # List of patchnames
+    
+    bpatch.patch_name1.patchtype=full-boundary             # patchtype one of "full-boundary", "circle, "rectangle", "circle-annular" or "rectangle-annular"	 
+    bpatch.patch_name1.boundary_direction=2                # patch normal direction 
+    bpatch.patch_name1.boundary_lo_or_hi=0                 # patch in low or high side of boundary    
+    bpatch.patch_name1.species= O2 N2                      # list of species names
+    
+    bpatch.patch_name2.patchtype=circle                    # patchtype one of "full-boundary", "circle, "rectangle", "circle-annular" or "rectangle-annular"	 
+    bpatch.patch_name2.boundary_direction=2                # patch normal direction 
+    bpatch.patch_name2.boundary_lo_or_hi=0                 # patch in low or high side of boundary   
+    bpatch.patch_name2.patch_circle_radius=0.1             # radius of the patch
+    bpatch.patch_name2.patch_circle_center=0.0 0.0 0.0     # coordinates of patch center 
+    bpatch.patch_name2.species= O2 N2                      # list of species names
+    
+    bpatch.patch_name3.patchtype=rectangle      	 
+    bpatch.patch_name3.boundary_direction=2                # patch normal direction 
+    bpatch.patch_name3.boundary_lo_or_hi=0                 # patch in low or high side of boundary   
+    bpatch.patch_name3.patch_rectangle_lo=0.0 0.0 0.0      # coordinates of low corner of rectangle
+    bpatch.patch_name3.patch_rectangle_hi=1.0 1.0 1.0      # coordinates of high corner of rectangle
+    bpatch.patch_name3.species= O2 N2                      # list of species names
+    
+    bpatch.patch_name4.patchtype=circle-annular  	 
+    bpatch.patch_name4.boundary_direction=2                # patch normal direction 
+    bpatch.patch_name4.boundary_lo_or_hi=0                 # patch in low or high side of boundary   
+    bpatch.patch_name4.patch_circ_ann_center= 0.0 0.0 0.0  # center of annular circle
+    bpatch.patch_name4.patch_circ_ann_inner_radius=0.1     # coordinates of patch center
+    bpatch.patch_name4.patch_circ_ann_outer_radius=0.2     # coordinates of patch center 
+    bpatch.patch_name4.species= O2 N2                      # list of species names
+    
+    bpatch.patch_name5.patchtype=rectangle-annular         
+    bpatch.patch_name5.boundary_direction=2                     # patch normal direction 
+    bpatch.patch_name5.boundary_lo_or_hi=0                      # patch in low or high side of boundary   
+    bpatch.patch_name5.patch_rect_ann_outer_lo = -1.0 -1.0 -1.0 # coordinates of low corner of outer rectangle 
+    bpatch.patch_name5.patch_rect_ann_outer_hi =  1.0  1.0  1.0 # coordinates of high corner of outer rectangle  
+    bpatch.patch_name5.patch_rect_ann_inner_lo = -0.5 -0.5 -0.5 # coordinates of low corner of inner rectangle 
+    bpatch.patch_name5.patch_rect_ann_inner_hi =  0.5  0.5  0.5 # coordinates of high corner of inner rectangle 
+    bpatch.patch_name5.species= O2 N2                           # list of species names
 
 The `do_temporal` flag will trigger the creation of a `temporals` folder in your run directory and the following entries
 will be appended to an ASCII `temporals/tempState` file: step, time, dt, kin. energy integral, enstrophy integral, mean pressure
@@ -434,6 +472,10 @@ turn on state extremas (stored in `temporals/tempExtremas` as min/max for each s
 `temporals/tempMass`) computing the total mass, dMdt and advective mass fluxes across the domain boundaries as well as the error in
 the balance (dMdt - sum of fluxes), and species balance (stored in `temporals/tempSpec`) computing each species total mass, dM_Ydt,
 advective \& diffusive fluxes across the domain boundaries, consumption rate integral and the error (dMdt - sum of fluxes - reaction).
+Users can also monitor species advective fluxes through specific regions of the domain boundaries (called as boundary patches).
+Patches can be defined on the low or high sides of non-embedded boundaries through the use of pre-defined shapes such as `circle`,
+`rectangle`,`circle-annular`, `rectangle-annular` and `full-boundary`. The zero AMR level, advective fluxes of each of the user-specified species will be 
+reported in the ASCII `temppatchmfr` file in the temporals folder. 
 
 Combustion diagnostics often involve the use of a mixture fraction and/or a progress variable, both of which can be defined
 at run time and added to the derived variables included in the plotfile. If `mixture_fraction` or `progress_variable` is

Exec/RegTests/TaylorGreen/input.3d_PLM

@@ -58,3 +58,11 @@ amr.derive_plot_vars = avg_pressure mag_vort
 #amrex.fpe_trap_invalid = 1
 #amrex.fpe_trap_zero = 1
 #amrex.fpe_trap_overflow = 1
+
+peleLM.do_patch_mfr=1
+peleLM.bpatch.patchnames= full_z
+
+bpatch.full_z.patchtype=full-boundary  #(full-boundary,circle,rectangle,circle-annular,rectangle-annular)
+bpatch.full_z.boundary_direction=2
+bpatch.full_z.boundary_lo_or_hi=0
+bpatch.full_z.species= O2 N2 

Exec/RegTests/TripleFlame/input.2d-regt

@@ -107,3 +107,16 @@ amr.HR.field_name = HeatRelease
 #amrex.fpe_trap_invalid = 1
 #amrex.fpe_trap_zero = 1
 #amrex.fpe_trap_overflow = 1
+
+peleLM.do_temporals = 1
+peleLM.temporal_int = 1
+peleLM.do_patch_mfr=1
+
+peleLM.bpatch.patchnames= outlet
+
+bpatch.outlet.patchtype=line  #(full-boundary,circle,rectangle,circle-annular,rectangle-annular)
+bpatch.outlet.patch_line_radius = 0.005
+bpatch.outlet.patch_line_center = 0.005 0
+bpatch.outlet.boundary_direction=1
+bpatch.outlet.boundary_lo_or_hi=0
+bpatch.outlet.species= O2 N2 CH4 

Source/Make.package

@@ -8,6 +8,7 @@ CEXE_headers += PeleLMeX_DiffusionOp.H
 CEXE_headers += PeleLMeX_DeriveFunc.H
 CEXE_headers += PeleLMeX_EBUserDefined.H
 CEXE_headers += PeleLMeX_FlowControllerData.H
+CEXE_headers += PeleLMeX_BPatch.H
 
 ## Sources
 CEXE_sources += main.cpp
@@ -41,6 +42,7 @@ CEXE_sources += PeleLMeX_EB.cpp
 CEXE_sources += PeleLMeX_Diagnostics.cpp
 CEXE_sources += PeleLMeX_FlowController.cpp
 CEXE_sources += PeleLMeX_DeriveUserDefined.cpp
+CEXE_sources += PeleLMeX_BPatch.cpp
 
 ifeq ($(USE_SOOT), TRUE)
   CEXE_sources += PeleLMeX_Soot.cpp

Source/PeleLMeX.H

@@ -9,6 +9,8 @@
 #include "PeleLMeX_DiffusionOp.H"
 #include "DiagBase.H"
 #include "PeleLMeX_FlowControllerData.H"
+#include "PeleLMeX_BPatch.H"
+
 #ifdef PELE_USE_EFIELD
 #include "PrecondOp.H"
 #endif
@@ -45,6 +47,8 @@ class SprayParticleContainer;
 class SootModel;
 #endif
 
+class BPatch;
+
 const std::string PrettyLine = " " + std::string(78, '=') + "\n";
 
 class PeleLM : public amrex::AmrCore
@@ -977,6 +981,9 @@ public:
   // Temporal
   void massBalance();
   void speciesBalance();
+  void speciesBalancePatch();
+  void initBPatches(amrex::Geometry& a_geom);
+
   void rhoHBalance();
   void initTemporals(const PeleLM::TimeStamp& a_time = AmrOldTime);
   void writeTemporals();
@@ -990,6 +997,10 @@ public:
     const amrex::Array<const amrex::MultiFab*, AMREX_SPACEDIM>& a_fluxes,
     const amrex::Geometry& a_geom,
     const amrex::Real& a_factor = 1.0);
+  void addRhoYFluxesPatch(
+    const amrex::Array<const amrex::MultiFab*, AMREX_SPACEDIM>& a_fluxes,
+    const amrex::Geometry& a_geom,
+    const amrex::Real& a_factor = 1.0);
   void addUmacFluxes(
     std::unique_ptr<AdvanceAdvData>& advData, const amrex::Geometry& a_geom);
   void openTempFile();
@@ -1685,6 +1696,10 @@ public:
   ProbParm* prob_parm = nullptr;
   ProbParm* prob_parm_d = nullptr;
 
+  // Patch parameters
+  BPatch* boundary_patches = nullptr;
+  BPatch* boundary_patches_d = nullptr;
+
 #ifdef PELE_USE_SOOT
   // Soot parameters
   SootModel* soot_model = nullptr;
@@ -1900,6 +1915,8 @@ public:
 
   // Temporals
   int m_do_temporals = 0;
+  int m_do_patch_mfr = 0;
+  amrex::Vector<std::unique_ptr<BPatch>> m_bPatches;
   int m_temp_int = 5;
   int m_do_extremas = 0;
   int m_do_massBalance = 0;
@@ -1920,6 +1937,7 @@ public:
   std::ofstream tmpExtremasFile;
   std::ofstream tmpMassFile;
   std::ofstream tmpSpecFile;
+  std::ofstream tmppatchmfrFile;
 
   // Number of ghost cells
 #ifdef AMREX_USE_EB

Source/PeleLMeX_Advection.cpp

@@ -661,6 +661,10 @@ PeleLM::computeScalarAdvTerms(std::unique_ptr<AdvanceAdvData>& advData)
     if (m_do_speciesBalance != 0) {
       addRhoYFluxes(GetArrOfConstPtrs(fluxes[0]), geom[0]);
     }
+
+    if (m_do_patch_mfr != 0) {
+      addRhoYFluxesPatch(GetArrOfConstPtrs(fluxes[0]), geom[0]);
+    }
   }
   // Compute face domain integral for U at every SDC iteration
   addUmacFluxes(advData, geom[0]);