diff --git a/.gitignore b/.gitignore index bc345021..80f0e71c 100644 --- a/.gitignore +++ b/.gitignore @@ -16,3 +16,10 @@ temporals # OSX files .DS_Store .AppleDouble + +#Eclipse file +.cproject +.project +Exec/**/*.cproject +Exec/**/*.project + diff --git a/CMake/BuildPeleExe.cmake b/CMake/BuildPeleExe.cmake index 1d45d570..6dc1cc0d 100644 --- a/CMake/BuildPeleExe.cmake +++ b/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 diff --git a/Docs/sphinx/manual/LMeXControls.rst b/Docs/sphinx/manual/LMeXControls.rst index 0c70e5bd..859f65f4 100644 --- a/Docs/sphinx/manual/LMeXControls.rst +++ b/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= # 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 diff --git a/Exec/RegTests/TaylorGreen/input.3d_PLM b/Exec/RegTests/TaylorGreen/input.3d_PLM index 11628df2..6427863d 100644 --- a/Exec/RegTests/TaylorGreen/input.3d_PLM +++ b/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 diff --git a/Exec/RegTests/TripleFlame/input.2d-regt b/Exec/RegTests/TripleFlame/input.2d-regt index a311eb8c..ff00a390 100644 --- a/Exec/RegTests/TripleFlame/input.2d-regt +++ b/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 diff --git a/Source/Make.package b/Source/Make.package index fa8d9f7e..1532ea82 100644 --- a/Source/Make.package +++ b/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 diff --git a/Source/PeleLMeX.H b/Source/PeleLMeX.H index 7e2d736a..549bb228 100644 --- a/Source/PeleLMeX.H +++ b/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& a_fluxes, const amrex::Geometry& a_geom, const amrex::Real& a_factor = 1.0); + void addRhoYFluxesPatch( + const amrex::Array& a_fluxes, + const amrex::Geometry& a_geom, + const amrex::Real& a_factor = 1.0); void addUmacFluxes( std::unique_ptr& 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> 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 diff --git a/Source/PeleLMeX_Advection.cpp b/Source/PeleLMeX_Advection.cpp index 85249a05..45d9d326 100644 --- a/Source/PeleLMeX_Advection.cpp +++ b/Source/PeleLMeX_Advection.cpp @@ -661,6 +661,10 @@ PeleLM::computeScalarAdvTerms(std::unique_ptr& 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]); diff --git a/Source/PeleLMeX_BPatch.H b/Source/PeleLMeX_BPatch.H new file mode 100644 index 00000000..69a099df --- /dev/null +++ b/Source/PeleLMeX_BPatch.H @@ -0,0 +1,239 @@ +#ifndef BPATCH_H +#define BPATCH_H + +#include + +class BPatch +{ + +public: + BPatch() = default; + ~BPatch() = default; + + BPatch(const std::string& patch_name, const amrex::Geometry& geom); + + struct BpatchDataContainer + { + int m_boundary_dir; + int m_boundary_lo_hi; + int num_species; + int m_patchtype_num; + int* speciesIndex = nullptr; + amrex::Real* speciesFlux = nullptr; + + // In 2D, a boundary patch is a line + amrex::GpuArray m_patch_line_center; + amrex::Real m_patch_line_radius; + + // circle variables + amrex::GpuArray m_patch_circle_center; + amrex::Real m_patch_circle_radius; + + // rectangle variables + amrex::GpuArray m_patch_rectangle_lo; + amrex::GpuArray m_patch_rectangle_hi; + + // circle annular variables + amrex::GpuArray m_patch_circ_ann_center; + amrex::Real m_patch_circ_ann_inner_radius; + amrex::Real m_patch_circ_ann_outer_radius; + + // Rectangular annular variables + amrex::GpuArray m_patch_rect_ann_outer_lo; + amrex::GpuArray m_patch_rect_ann_outer_hi; + amrex::GpuArray m_patch_rect_ann_inner_lo; + amrex::GpuArray m_patch_rect_ann_inner_hi; + + AMREX_GPU_DEVICE + AMREX_FORCE_INLINE + bool CheckifPointInside( + amrex::GpuArray point_coordinate, + amrex::Real dx) const + { + bool inside = false; + + // auto dx = geom.CellSizeArray(); + amrex::Real xp, yp; + + point_coordinate[m_boundary_dir] = 0.0; +#if (AMREX_SPACEDIM == 2) + xp = point_coordinate[0]; + yp = point_coordinate[1]; + if (m_patchtype_num == 0) { + inside = true; + } else if (m_patchtype_num == 1) { + amrex::Real patch_line_radius_touse = m_patch_circle_radius + dx / 2.0; + amrex::Real rad = sqrt( + (xp - m_patch_circle_center[0]) * (xp - m_patch_circle_center[0]) + + (yp - m_patch_circle_center[1]) * (yp - m_patch_circle_center[1])); + if (rad <= patch_line_radius_touse) { + inside = true; + } + } +#elif (AMREX_SPACEDIM == 3) + amrex::Real zp; + const amrex::Real sqrt2 = sqrt(2.0); + + xp = point_coordinate[0]; + yp = point_coordinate[1]; + zp = point_coordinate[2]; + + if (m_patchtype_num == 0) { + inside = true; + } else if (m_patchtype_num == 1) { + + amrex::Real patch_circle_radius_touse = + m_patch_circle_radius + dx * sqrt2 / 2.0; + amrex::Real rad = sqrt( + (xp - m_patch_circle_center[0]) * (xp - m_patch_circle_center[0]) + + (yp - m_patch_circle_center[1]) * (yp - m_patch_circle_center[1]) + + (zp - m_patch_circle_center[2]) * (zp - m_patch_circle_center[2])); + if (rad <= patch_circle_radius_touse) { + inside = true; + } + + } else if (m_patchtype_num == 2) { + amrex::GpuArray patch_rectangle_lo_touse = + {AMREX_D_DECL( + m_patch_rectangle_lo[0] - dx / 2.0, + m_patch_rectangle_lo[1] - dx / 2.0, + m_patch_rectangle_lo[2] - dx / 2.0)}; + amrex::GpuArray patch_rectangle_hi_touse = + {AMREX_D_DECL( + m_patch_rectangle_hi[0] + dx / 2.0, + m_patch_rectangle_hi[1] + dx / 2.0, + m_patch_rectangle_hi[2] + dx / 2.0)}; + patch_rectangle_hi_touse[m_boundary_dir] = 0.0; + patch_rectangle_lo_touse[m_boundary_dir] = 0.0; + + if ( + (xp >= patch_rectangle_lo_touse[0] && + xp <= patch_rectangle_hi_touse[0]) && + (yp >= patch_rectangle_lo_touse[1] && + yp <= patch_rectangle_hi_touse[1]) && + (zp >= patch_rectangle_lo_touse[2] && + zp <= patch_rectangle_hi_touse[2])) { + inside = true; + } + } else if (m_patchtype_num == 3) { + amrex::Real patch_circ_ann_outer_radius_touse = + m_patch_circ_ann_outer_radius + dx * sqrt2 / 2.0; + amrex::Real patch_circ_ann_inner_radius_touse = + m_patch_circ_ann_inner_radius - dx * sqrt2 / 2.0; + + amrex::Real rad = sqrt( + (xp - m_patch_circ_ann_center[0]) * + (xp - m_patch_circ_ann_center[0]) + + (yp - m_patch_circ_ann_center[1]) * + (yp - m_patch_circ_ann_center[1]) + + (zp - m_patch_circ_ann_center[2]) * + (zp - m_patch_circ_ann_center[2])); + if ( + rad <= patch_circ_ann_outer_radius_touse && + rad >= patch_circ_ann_inner_radius_touse) { + inside = true; + } + } else if (m_patchtype_num == 4) { + bool inside_outer = false; + bool inside_inner = false; + amrex::GpuArray + patch_rect_ann_inner_lo_touse = {AMREX_D_DECL( + m_patch_rect_ann_inner_lo[0] + dx / 2.0, + m_patch_rect_ann_inner_lo[1] + dx / 2.0, + m_patch_rect_ann_inner_lo[2] + dx / 2.0)}; + amrex::GpuArray + patch_rect_ann_inner_hi_touse = {AMREX_D_DECL( + m_patch_rect_ann_inner_hi[0] - dx / 2.0, + m_patch_rect_ann_inner_hi[1] - dx / 2.0, + m_patch_rect_ann_inner_hi[2] - dx / 2.0)}; + amrex::GpuArray + patch_rect_ann_outer_lo_touse = {AMREX_D_DECL( + m_patch_rect_ann_outer_lo[0] - dx / 2.0, + m_patch_rect_ann_outer_lo[1] - dx / 2.0, + m_patch_rect_ann_outer_lo[2] - dx / 2.0)}; + amrex::GpuArray + patch_rect_ann_outer_hi_touse = {AMREX_D_DECL( + m_patch_rect_ann_outer_hi[0] + dx / 2.0, + m_patch_rect_ann_outer_hi[1] + dx / 2.0, + m_patch_rect_ann_outer_hi[2] + dx / 2.0)}; + patch_rect_ann_inner_lo_touse[m_boundary_dir] = 0.0; + patch_rect_ann_inner_hi_touse[m_boundary_dir] = 0.0; + patch_rect_ann_outer_lo_touse[m_boundary_dir] = 0.0; + patch_rect_ann_outer_hi_touse[m_boundary_dir] = 0.0; + + if (((xp >= patch_rect_ann_outer_lo_touse[0] && + xp <= patch_rect_ann_outer_hi_touse[0]) && + (yp >= patch_rect_ann_outer_lo_touse[1] && + yp <= patch_rect_ann_outer_hi_touse[1]) && + (zp >= patch_rect_ann_outer_lo_touse[2] && + zp <= patch_rect_ann_outer_hi_touse[2]))) { + inside_outer = true; + } + if (((xp >= patch_rect_ann_inner_lo_touse[0] && + xp <= patch_rect_ann_inner_hi_touse[0]) && + (yp >= patch_rect_ann_inner_lo_touse[1] && + yp <= patch_rect_ann_inner_hi_touse[1]) && + (zp >= patch_rect_ann_inner_lo_touse[2] && + zp <= patch_rect_ann_inner_hi_touse[2]))) { + inside_inner = true; + } + if (inside_outer && !inside_inner) { + inside = true; + } + } +#endif + + return inside; + } + }; + + void allocate() + { + if (!m_device_allocated) { + m_bpdata_d = + (BpatchDataContainer*)amrex::The_Arena()->alloc(sizeof(m_bpdata_h)); + m_device_allocated = true; + sync_to_device(); + } + } + + void deallocate() + { + if (m_host_allocated) { + amrex::The_Pinned_Arena()->free(m_bpdata_h.speciesIndex); + amrex::The_Pinned_Arena()->free(m_bpdata_h.speciesFlux); + } + if (m_device_allocated) { + amrex::The_Arena()->free(m_bpdata_d); + } + } + + void sync_to_device() + { + if (!m_device_allocated) { + amrex::Abort("Device params not allocated yet"); + } else { + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, &m_bpdata_h, &m_bpdata_h + 1, m_bpdata_d); + } + } + + // Accessors + BpatchDataContainer getHostData() { return m_bpdata_h; } + BpatchDataContainer* getHostDataPtr() { return &m_bpdata_h; } + BpatchDataContainer* getDeviceData() { return m_bpdata_d; } + + std::string m_patchname; + std::string + m_patchtype; // values: + // "fullboundary","circle","rectangle","annular-circle","annular-rectangle" + amrex::Vector speciesList; + +protected: + BpatchDataContainer m_bpdata_h; + BpatchDataContainer* m_bpdata_d = nullptr; + bool m_device_allocated{false}; + bool m_host_allocated{false}; +}; + +#endif diff --git a/Source/PeleLMeX_BPatch.cpp b/Source/PeleLMeX_BPatch.cpp new file mode 100644 index 00000000..527e611b --- /dev/null +++ b/Source/PeleLMeX_BPatch.cpp @@ -0,0 +1,181 @@ +#include "PeleLMeX_BPatch.H" +#include "PelePhysics.H" + +BPatch::BPatch(const std::string& patch_name, const amrex::Geometry& geom) + : m_patchname(std::move(patch_name)) +{ + + std::string ppbpatch = "bpatch." + m_patchname; + amrex::ParmParse ps(ppbpatch); + + auto prob_lo = geom.ProbLoArray(); + auto prob_hi = geom.ProbHiArray(); + + // Reading patch related info + ps.get("patchtype", m_patchtype); + ps.get("boundary_direction", m_bpdata_h.m_boundary_dir); + ps.get("boundary_lo_or_hi", m_bpdata_h.m_boundary_lo_hi); + + // Checking if patch type matches defined ones +#if (AMREX_SPACEDIM == 1) + if (m_patchtype != "full-boundary") { + amrex::Abort( + "\nPatch type for 1-D domains should be specified as full-boundary.\n"); + } +#elif (AMREX_SPACEDIM == 2) + if (m_patchtype != "line" && m_patchtype != "full-boundary") { + amrex::Abort( + "\nPatch type for 2-D domains should be line or full-boundary\n"); + } +#else + if ( + m_patchtype != "full-boundary" && m_patchtype != "circle" && + m_patchtype != "rectangle" && m_patchtype != "circle-annular" && + m_patchtype != "rectangle-annular") { + amrex::Abort( + "\nPatch type do not match allowed values " + "(full-boundary,circle,rectangle,circle-annular,rectangle-annular) \n"); + } +#endif + + if ( + m_bpdata_h.m_boundary_dir != 0 && m_bpdata_h.m_boundary_dir != 1 && + m_bpdata_h.m_boundary_dir != 2) { + amrex::Abort("\nBoundary direction should be 0,1 or 2"); + } + + if (m_bpdata_h.m_boundary_lo_hi != 0 && m_bpdata_h.m_boundary_lo_hi != 1) { + amrex::Abort("\nBoundary high low should be 0 or 1"); + } + +#if (AMREX_SPACEDIM == 2) + if (m_bpdata_h.m_boundary_dir == 2) { + amrex::Abort("\nFor 2D problems, patch boundary direction cannot be 2"); + } +#endif + + // Define patch_num + if (m_patchtype == "full-boundary") { + m_bpdata_h.m_patchtype_num = 0; + } else if (m_patchtype == "line" || m_patchtype == "circle") { + m_bpdata_h.m_patchtype_num = 1; + } else if (m_patchtype == "rectangle") { + m_bpdata_h.m_patchtype_num = 2; + } else if (m_patchtype == "circle-annular") { + m_bpdata_h.m_patchtype_num = 3; + } else if (m_patchtype == "rectangle-annular") { + m_bpdata_h.m_patchtype_num = 4; + } else { + amrex::Abort("\nError! Unknown patch type"); + } + + // Define patch variables + if (m_patchtype == "line") { + for (int n = 0; n < 2; ++n) { + ps.get("patch_line_center", m_bpdata_h.m_patch_line_center[n], n); + } + m_bpdata_h.m_patch_line_center[m_bpdata_h.m_boundary_dir] = 0.0; + ps.get("patch_line_radius", m_bpdata_h.m_patch_line_radius); + if (m_bpdata_h.m_patch_line_radius <= 0.0) { + amrex::Abort("\nPatch radius for line should be greater than 0"); + } + } else if (m_patchtype == "full-boundary") { + for (int n = 0; n < AMREX_SPACEDIM; ++n) { + m_bpdata_h.m_patch_rectangle_lo[n] = prob_lo[n]; + m_bpdata_h.m_patch_rectangle_hi[n] = prob_hi[n]; + } + m_bpdata_h.m_patch_rectangle_lo[m_bpdata_h.m_boundary_dir] = 0.0; + m_bpdata_h.m_patch_rectangle_hi[m_bpdata_h.m_boundary_dir] = 0.0; + } else if (m_patchtype == "circle") { + for (int n = 0; n < 3; ++n) { + ps.get("patch_circle_center", m_bpdata_h.m_patch_circle_center[n], n); + } + m_bpdata_h.m_patch_circle_center[m_bpdata_h.m_boundary_dir] = 0.0; + ps.get("patch_circle_radius", m_bpdata_h.m_patch_circle_radius); + if (m_bpdata_h.m_patch_circle_radius <= 0.0) { + amrex::Abort("\nPatch radius for circle should be greater than 0"); + } + } else if (m_patchtype == "rectangle") { + for (int n = 0; n < 3; ++n) { + ps.get("patch_rectangle_lo", m_bpdata_h.m_patch_rectangle_lo[n], n); + } + for (int n = 0; n < 3; ++n) { + ps.get("patch_rectangle_hi", m_bpdata_h.m_patch_rectangle_hi[n], n); + } + m_bpdata_h.m_patch_rectangle_lo[m_bpdata_h.m_boundary_dir] = 0.0; + m_bpdata_h.m_patch_rectangle_hi[m_bpdata_h.m_boundary_dir] = 0.0; + } else if (m_patchtype == "circle-annular") { + for (int n = 0; n < 3; ++n) { + ps.get("patch_circ_ann_center", m_bpdata_h.m_patch_circ_ann_center[n], n); + } + m_bpdata_h.m_patch_circ_ann_center[m_bpdata_h.m_boundary_dir] = 0.0; + ps.get( + "patch_circ_ann_inner_radius", m_bpdata_h.m_patch_circ_ann_inner_radius); + ps.get( + "patch_circ_ann_outer_radius", m_bpdata_h.m_patch_circ_ann_outer_radius); + } else if (m_patchtype == "rectangle-annular") { + for (int n = 0; n < 3; ++n) { + ps.get( + "patch_rect_ann_outer_lo", m_bpdata_h.m_patch_rect_ann_outer_lo[n], n); + } + for (int n = 0; n < 3; ++n) { + ps.get( + "patch_rect_ann_outer_hi", m_bpdata_h.m_patch_rect_ann_outer_hi[n], n); + } + for (int n = 0; n < 3; ++n) { + ps.get( + "patch_rect_ann_inner_lo", m_bpdata_h.m_patch_rect_ann_inner_lo[n], n); + } + for (int n = 0; n < 3; ++n) { + ps.get( + "patch_rect_ann_inner_hi", m_bpdata_h.m_patch_rect_ann_inner_hi[n], n); + } + + m_bpdata_h.m_patch_rect_ann_inner_hi[m_bpdata_h.m_boundary_dir] = 0.0; + m_bpdata_h.m_patch_rect_ann_outer_hi[m_bpdata_h.m_boundary_dir] = 0.0; + m_bpdata_h.m_patch_rect_ann_inner_lo[m_bpdata_h.m_boundary_dir] = 0.0; + m_bpdata_h.m_patch_rect_ann_outer_lo[m_bpdata_h.m_boundary_dir] = 0.0; + } else { + amrex::Abort("\nError! Unknown patch type"); + } + + m_bpdata_h.num_species = ps.countval("species"); + ps.getarr("species", speciesList); + + if (m_bpdata_h.num_species > 0) { + speciesList.resize(m_bpdata_h.num_species); + m_bpdata_h.speciesIndex = (int*)amrex::The_Pinned_Arena()->alloc( + m_bpdata_h.num_species * sizeof(int)); + m_bpdata_h.speciesFlux = (amrex::Real*)amrex::The_Pinned_Arena()->alloc( + m_bpdata_h.num_species * sizeof(amrex::Real)); + } else { + amrex::Abort("\nError! No species provided to plot flux at boundary patch"); + } + + amrex::Vector names; + pele::physics::eos::speciesNames(names); + names.resize(names.size()); + + for (int n = 0; n < names.size(); n++) { + m_bpdata_h.speciesIndex[n] = -1; + } + + for (int m = 0; m < m_bpdata_h.num_species; m++) { + for (int n = 0; n < names.size(); n++) { + if (speciesList[m] == names[n]) { + m_bpdata_h.speciesIndex[m] = n; + } + } + } + + for (int n = 0; n < m_bpdata_h.num_species; n++) { + if (m_bpdata_h.speciesIndex[n] == -1) { + std::string msg = + "\nError! Unable to find species index " + std::to_string(n); + amrex::Abort(msg); + } + } + + allocate(); + amrex::Gpu::streamSynchronize(); +} diff --git a/Source/PeleLMeX_Setup.cpp b/Source/PeleLMeX_Setup.cpp index 2f5bf570..deb511d1 100644 --- a/Source/PeleLMeX_Setup.cpp +++ b/Source/PeleLMeX_Setup.cpp @@ -1,6 +1,7 @@ #include #include #include +#include #include "PelePhysics.H" #include #ifdef PELE_USE_EFIELD @@ -44,8 +45,10 @@ 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]), &&amrex::almostEqual(dx[1], dx[2]))); + , amrex::almostEqual(dx[0], dx[1], 10), + &&amrex::almostEqual(dx[1], dx[2], 10))); } // Print build info to screen const char* githash1 = buildInfoGetGitHash(1); @@ -93,6 +96,11 @@ PeleLM::Setup() // Diagnostics setup createDiagnostics(); + // Boundary Patch Setup + if (m_do_patch_mfr != 0) { + initBPatches(Geom(0)); + } + // Initialize Level Hierarchy data resizeArray(); @@ -540,6 +548,7 @@ PeleLM::readParameters() pp.query("do_extremas", m_do_extremas); pp.query("do_mass_balance", m_do_massBalance); pp.query("do_species_balance", m_do_speciesBalance); + pp.query("do_patch_mfr", m_do_patch_mfr); } // ----------------------------------------- diff --git a/Source/PeleLMeX_Temporals.cpp b/Source/PeleLMeX_Temporals.cpp index e738fb08..1f5fd377 100644 --- a/Source/PeleLMeX_Temporals.cpp +++ b/Source/PeleLMeX_Temporals.cpp @@ -1,4 +1,5 @@ #include +#include using namespace amrex; @@ -55,6 +56,20 @@ PeleLM::massBalance() tmpMassFile.flush(); } +void +PeleLM::speciesBalancePatch() +{ + tmppatchmfrFile << m_nstep << " " << m_cur_time; // Time info + for (int n = 0; n < m_bPatches.size(); n++) { + BPatch::BpatchDataContainer* bphost = m_bPatches[n]->getHostDataPtr(); + for (int i = 0; i < bphost->num_species; i++) { + tmppatchmfrFile << " " << bphost->speciesFlux[i]; + } + } + tmppatchmfrFile << "\n"; + tmppatchmfrFile.flush(); +} + void PeleLM::speciesBalance() { @@ -501,6 +516,118 @@ PeleLM::addRhoYFluxes( } } +void +PeleLM::initBPatches(Geometry& a_geom) +{ + std::string pele_prefix = "peleLM.bpatch"; + ParmParse pp(pele_prefix); + int num_bPatches = 0; + num_bPatches = pp.countval("patchnames"); + + Vector bpatch_name; + if (num_bPatches > 0) { + + m_bPatches.resize(num_bPatches); + bpatch_name.resize(num_bPatches); + } + for (int n = 0; n < num_bPatches; ++n) { + pp.get("patchnames", bpatch_name[n], n); + m_bPatches[n] = std::make_unique(bpatch_name[n], a_geom); + if (m_verbose > 0) { + Print() << " Initializing boundary patch: " << bpatch_name[n] + << std::endl; + } + } +} + +void +PeleLM::addRhoYFluxesPatch( + const Array& a_fluxes, + const Geometry& a_geom, + const Real& a_factor) +{ + + if (!(m_nstep % m_temp_int == m_temp_int - 1)) { + return; + } + + const auto dx = a_geom.CellSizeArray(); + auto prob_lo = a_geom.ProbLoArray(); + amrex::GpuArray area; + +#if (AMREX_SPACEDIM == 1) + area[0] = 1.0; +#elif (AMREX_SPACEDIM == 2) + if (geom[0].IsRZ() && m_bPatches.size() > 0) { + Abort("Bpatches not supported in RZ coordinates"); + } + area[0] = dx[1]; + area[1] = dx[0]; +#else + area[0] = dx[1] * dx[2]; + area[1] = dx[0] * dx[2]; + area[2] = dx[0] * dx[1]; +#endif + + // Loop through all patches + for (int n = 0; n < m_bPatches.size(); n++) { + + BPatch* patch = m_bPatches[n].get(); + BPatch::BpatchDataContainer const* bpdevice = patch->getDeviceData(); + BPatch::BpatchDataContainer const* bphost = patch->getHostDataPtr(); + const int idim = bphost->m_boundary_dir; + + auto faceDomain = + amrex::convert(a_geom.Domain(), IntVect::TheDimensionVector(idim)); + auto const& fma = a_fluxes[idim]->const_arrays(); + + // Loop through species specified by user + for (int m = 0; m < bphost->num_species; m++) { + + Real sum_species_flux_global = 0.0; + + { + auto r = amrex::ParReduce( + TypeList{}, TypeList{}, + *a_fluxes[idim], IntVect(0), + [=] AMREX_GPU_DEVICE( + int box_no, int i, int j, int k) noexcept -> GpuTuple { + Array4 const& flux = fma[box_no]; + int idx = + (bpdevice->m_boundary_dir == 0 + ? i + : (bpdevice->m_boundary_dir == 1 ? j : k)); + int idx_lo_hi = + (bpdevice->m_boundary_lo_hi == 0 ? faceDomain.smallEnd(idim) + : faceDomain.bigEnd(idim)); + + amrex::GpuArray point_coordinates{ + AMREX_D_DECL( + prob_lo[0] + (i + 0.5) * dx[0], prob_lo[1] + (j + 0.5) * dx[1], + prob_lo[2] + (k + 0.5) * dx[2])}; + + Real sum_species_flux = 0.0; + Real dummy = 0.0; + const bool ifinside = + bpdevice->CheckifPointInside(point_coordinates, dx[0]); + + if (idx == idx_lo_hi and ifinside) { + int species_idx = bpdevice->speciesIndex[m]; + sum_species_flux += flux(i, j, k, species_idx) * area[idim]; + } + return {sum_species_flux, dummy}; + }); + sum_species_flux_global = amrex::get<0>(r); + + ParallelAllReduce::Sum( + {sum_species_flux_global}, ParallelContext::CommunicatorSub()); + bphost->speciesFlux[m] = a_factor * sum_species_flux_global; + // amrex::Print()<<"\nNew func = "<getHostData(); + for (int i = 0; i < bphost.num_species; i++) { + tmppatchmfrFile << "," + << patch->m_patchname + "_" + patch->speciesList[i]; + } + } + tmppatchmfrFile << "\n"; + } #ifdef PELE_USE_EFIELD if (m_do_ionsBalance) { tempFileName = "temporals/tempIons"; @@ -644,6 +793,10 @@ PeleLM::closeTempFile() tmpExtremasFile.flush(); tmpExtremasFile.close(); } + if (m_do_patch_mfr != 0) { + tmppatchmfrFile.flush(); + tmppatchmfrFile.close(); + } #ifdef PELE_USE_EFIELD if (m_do_ionsBalance) { tmpIonsFile.flush();