From 746448330005d64ad96225b739b559b3cc3272ad Mon Sep 17 00:00:00 2001 From: Mukul Dave Date: Sun, 21 Jul 2024 16:47:06 -0700 Subject: [PATCH] Enforcing inflow-outflow boundary solvability (#129) Introduces routines to calculate outflux and influx at the in-out boundaries and then correct the outflux to match with influx. Three steps: 1. set masks for tagging cells at the boundary for inflow or outflow 2. calculate influx and outflux using reductions 3. correct outflow velocities using the correction factor --------- Co-authored-by: Ann Almgren --- .gitignore | 2 + Docs/source/InOutSolvability.rst | 28 ++ Docs/source/Utilities.rst | 1 + Utils/CMakeLists.txt | 1 + Utils/Make.package | 1 + Utils/hydro_enforce_inout_solvability.cpp | 300 ++++++++++++++++++++++ Utils/hydro_utils.H | 11 + 7 files changed, 344 insertions(+) create mode 100644 Docs/source/InOutSolvability.rst create mode 100644 Utils/hydro_enforce_inout_solvability.cpp diff --git a/.gitignore b/.gitignore index 632a8ac17..05485ceb2 100644 --- a/.gitignore +++ b/.gitignore @@ -54,3 +54,5 @@ Tutorials/EB/Donut/Exec/ls_plt* # local config Tools/GNUMake/Make.local + +.vscode* diff --git a/Docs/source/InOutSolvability.rst b/Docs/source/InOutSolvability.rst new file mode 100644 index 000000000..c74caa515 --- /dev/null +++ b/Docs/source/InOutSolvability.rst @@ -0,0 +1,28 @@ +.. include:: CustomCommands.rst + +Enforcing inflow-outflow solvability +------------------------------------ + +This routine enforces solvability for inflow-outflow boundaries, +which have both inflow and outflow cells. +A flux correction factor, :math:`\alpha_\text{fcf}` is introduced +which scales the outflow velocities to match the inflow: + +.. math:: + + \sum_{\partial\Omega_\text{in}} {\bf u} \cdot {\bf \area} + \alpha_\text{fcf} \sum_{\partial\Omega_\text{out}} {\bf u} \cdot {\bf \area} = 0 + + +The new flux-conserving velocities to be used for the MAC/nodal projections, +:math:`{\bf u}_\text{fc}`, are calculated from the correction factor as: + +.. math:: + \alpha_\text{fcf} = \frac{-\sum_{\partial\Omega_\text{in}} {\bf u} \cdot {\bf \area}}{\sum_{\partial\Omega_\text{out}} {\bf u} \cdot {\bf \area}}, + +.. math:: + {\bf u}_\text{fc} = \alpha_\text{fcf} \cdot {\bf u}, \ \forall {\bf x} \in \partial\Omega_\text{out}. + +It must be noted that this routine currently only accounts for boundaries +with the math BC ``BCType::direction_dependent``, which is to be used for +an inflow-outflow boundary. It does not compute or correct the boundary velocities +over other math BC types such as those representing pure inflow or pure outflow. \ No newline at end of file diff --git a/Docs/source/Utilities.rst b/Docs/source/Utilities.rst index 59f12aa7f..1380724f5 100644 --- a/Docs/source/Utilities.rst +++ b/Docs/source/Utilities.rst @@ -38,3 +38,4 @@ Note that EB is only an option for Cartesian geometries as of this writing. .. include:: AdvectiveTerm.rst +.. include:: InOutSolvability.rst diff --git a/Utils/CMakeLists.txt b/Utils/CMakeLists.txt index 0f3c7302f..40d3b0375 100644 --- a/Utils/CMakeLists.txt +++ b/Utils/CMakeLists.txt @@ -14,4 +14,5 @@ target_sources( hydro_utils.cpp hydro_constants.H hydro_bcs_K.H + hydro_enforce_inout_solvability.cpp ) diff --git a/Utils/Make.package b/Utils/Make.package index bdafec9f3..e0e59b5f6 100644 --- a/Utils/Make.package +++ b/Utils/Make.package @@ -1,6 +1,7 @@ CEXE_sources += hydro_utils.cpp CEXE_sources += hydro_compute_edgestate_and_flux.cpp CEXE_sources += hydro_extrap_vel_to_faces.cpp +CEXE_sources += hydro_enforce_inout_solvability.cpp CEXE_headers += hydro_bcs_K.H CEXE_headers += hydro_utils.H diff --git a/Utils/hydro_enforce_inout_solvability.cpp b/Utils/hydro_enforce_inout_solvability.cpp new file mode 100644 index 000000000..3a73f0e0c --- /dev/null +++ b/Utils/hydro_enforce_inout_solvability.cpp @@ -0,0 +1,300 @@ +/** \addtogroup Utilities + * @{ + */ + +#include + +using namespace amrex; + +namespace HydroUtils { + +namespace { + +void set_inout_masks( + const int lev, + const Vector>& vels_vec, + Array& inout_masks, + const BCRec* bc_type, + const Box& domain, + const bool corners) +{ + for (OrientationIter oit; oit != nullptr; ++oit) { + const auto ori = oit(); + const auto side = ori.faceDir(); + const int dir = ori.coordDir(); + const auto islow = ori.isLow(); + const auto ishigh = ori.isHigh(); + + // Multifab for normal velocity + const auto& vel_mf = vels_vec[lev][dir]; + + // mask iMF for the respective velocity direction + auto& inout_mask = inout_masks[dir]; + + IndexType::CellIndex dir_index_type = (vel_mf->ixType()).ixType(dir); + // domain extent indices for the velocities + int dlo; + if (dir_index_type == IndexType::CellIndex::CELL) { + // lower boundary is at -1 for cell-centered velocity + dlo = domain.smallEnd(dir) - 1; + } else { + // lower boundary is at 0 for face-centered velocity + dlo = domain.smallEnd(dir); + } + int dhi = domain.bigEnd(dir) + 1; + + // get BCs for the normal velocity and set the boundary index + // based on low or high side + const BCRec ibcrec = bc_type[dir]; + int bc, bndry; + if (side == Orientation::low) { + bc = ibcrec.lo(dir); + bndry = dlo; + } else { + bc = ibcrec.hi(dir); + bndry = dhi; + } + + // limit influx/outflux calculations to the in-out boundaries only + if (bc == BCType::direction_dependent) { + for (MFIter mfi(*vel_mf, TilingIfNotGPU()); mfi.isValid(); ++mfi) { + + Box box = mfi.validbox(); + + // include ghost cells along normal velocity for cell-centered + // not for face-centered as boundary lies in valid region + if (dir_index_type == IndexType::CellIndex::CELL) { + box.grow(dir, 1); + } + + // include boundary corners if specified + if (corners) { + box.grow((dir+1)%AMREX_SPACEDIM, 1); + box.grow((dir+2)%AMREX_SPACEDIM, 1); + } + + // Enter further only if the box bndry is at the domain bndry + if ((islow && (box.smallEnd(dir) == dlo)) + || (ishigh && (box.bigEnd(dir) == dhi))) { + + // create a 2D box normal to dir at the low/high bndry + Box box2d(box); box2d.setRange(dir, bndry); + + auto vel_arr = vel_mf->array(mfi); + auto inout_mask_arr = inout_mask.array(mfi); + + // tag cells as inflow or outflow by checking vel direction + ParallelFor(box2d, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + if ((side == Orientation::low && vel_arr(i,j,k) >= 0) + || (side == Orientation::high && vel_arr(i,j,k) <= 0)) { + inout_mask_arr(i,j,k) = -1; + } else { + inout_mask_arr(i,j,k) = +1; + } + }); + } + } + } + } +} + +void compute_influx_outflux( + const int lev, + const Vector>& vels_vec, + const Array& inout_masks, + const Real* a_dx, + Real& influx, + Real& outflux, + const bool corners) +{ + influx = 0.0, outflux = 0.0; + + for (int idim = 0; idim < AMREX_SPACEDIM; idim++) { + + // normal face area + const Real ds = + a_dx[(idim+1) % AMREX_SPACEDIM] * a_dx[(idim+2) % AMREX_SPACEDIM]; + + // Multifab for normal velocity + const auto& vel_mf = vels_vec[lev][idim]; + + // grow in the respective direction if vel is cell-centered + IndexType index_type = vel_mf->ixType(); + index_type.flip(idim); IntVect ngrow = index_type.ixType(); + + // grow in the transverse direction to include boundary corners + if (corners) { + ngrow[(idim+1)%AMREX_SPACEDIM] = 1; + ngrow[(idim+2)%AMREX_SPACEDIM] = 1; + } + + // mask iMF for the respective velocity direction + const auto& inout_mask = inout_masks[idim]; + + // define "multi-arrays" and perform reduction using the mask + auto const& vel_ma = vel_mf->const_arrays(); + auto const& inout_mask_ma = inout_mask.const_arrays(); + + influx += ds * + ParReduce(TypeList{}, + TypeList{}, + *vel_mf, ngrow, + [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) + noexcept -> GpuTuple + { + if (inout_mask_ma[box_no](i,j,k) == -1) { + return { std::abs(vel_ma[box_no](i,j,k)) }; + } else { + return { 0. }; + } + }); + + outflux += ds * + ParReduce(TypeList{}, + TypeList{}, + *vel_mf, ngrow, + [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) + noexcept -> GpuTuple + { + if (inout_mask_ma[box_no](i,j,k) == 1) { + return { std::abs(vel_ma[box_no](i,j,k)) }; + } else { + return { 0. }; + } + }); + } + ParallelDescriptor::ReduceRealSum(influx); + ParallelDescriptor::ReduceRealSum(outflux); +} + +void correct_outflow( + const int lev, + const Vector>& vels_vec, + const Array& inout_masks, + const BCRec* bc_type, + const Box& domain, + const Real alpha_fcf, + const bool corners) +{ + for (OrientationIter oit; oit != nullptr; ++oit) { + const auto ori = oit(); + const auto side = ori.faceDir(); + const int dir = ori.coordDir(); + const auto islow = ori.isLow(); + const auto ishigh = ori.isHigh(); + + // Multifab for normal velocity + const auto& vel_mf = vels_vec[lev][dir]; + + // mask iMF for the respective velocity direction + const auto& inout_mask = inout_masks[dir]; + + IndexType::CellIndex dir_index_type = (vel_mf->ixType()).ixType(dir); + // domain extent indices for the velocities + int dlo; + if (dir_index_type == IndexType::CellIndex::CELL) { + dlo = domain.smallEnd(dir) - 1; // cell-centered boundary + } else { + dlo = domain.smallEnd(dir); // face-centered boundary + } + int dhi = domain.bigEnd(dir) + 1; + + // get BCs for the normal velocity and set the boundary index + const BCRec ibcrec = bc_type[dir]; + int bc, bndry; + if (side == Orientation::low) { + bc = ibcrec.lo(dir); + bndry = dlo; + } else { + bc = ibcrec.hi(dir); + bndry = dhi; + } + + if (bc == BCType::direction_dependent) { + for (MFIter mfi(*vel_mf, TilingIfNotGPU()); mfi.isValid(); ++mfi) { + + Box box = mfi.validbox(); + if (dir_index_type == IndexType::CellIndex::CELL) { + box.grow(dir, 1); + } + if (corners) { + box.grow((dir+1)%AMREX_SPACEDIM, 1); + box.grow((dir+2)%AMREX_SPACEDIM, 1); + } + + // Enter further only if the box boundary is at the domain boundary + if ((islow && (box.smallEnd(dir) == dlo)) + || (ishigh && (box.bigEnd(dir) == dhi))) { + + // create a 2D box normal to dir at the low/high boundary + Box box2d(box); box2d.setRange(dir, bndry); + + auto vel_arr = vel_mf->array(mfi); + auto inout_mask_arr = inout_mask.array(mfi); + + ParallelFor(box2d, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + if (inout_mask_arr(i,j,k) == 1) { + vel_arr(i,j,k) *= alpha_fcf; + } + }); + } + } + } + } +} + +} // file-local namespace + +void enforceInOutSolvability ( + const Vector>& vels_vec, + const BCRec* bc_type, + const Vector& geom, + bool include_bndry_corners +) +{ + const Box domain = geom[0].Domain(); + + const auto nlevs = int(vels_vec.size()); + for (int lev = 0; lev < nlevs; ++lev) { + + // masks to tag inflow/outflow at the boundaries + // separate iMultifab for each velocity direction + // 0 for interior cells, + // -1 for inflow bndry cells, +1 for outflow bndry cells + Array inout_masks; + + // defining the mask iMultifabs in each direction + for (int idim = 0; idim < AMREX_SPACEDIM; idim++) + { + const auto& vel_mf = vels_vec[lev][idim]; // normal velocity multifab + + // grow in the respective direction if vel is cell-centered + // to include the boundary cells + IndexType index_type = vel_mf->ixType(); + index_type.flip(idim); IntVect ngrow = index_type.ixType(); + + // grow in the transverse direction to include boundary corners + if (include_bndry_corners) { + ngrow[(idim+1)%AMREX_SPACEDIM] = 1; + ngrow[(idim+2)%AMREX_SPACEDIM] = 1; + } + + inout_masks[idim].define(vel_mf->boxArray(), vel_mf->DistributionMap(), 1, ngrow); + inout_masks[idim].setVal(0); + } + + set_inout_masks(lev, vels_vec, inout_masks, bc_type, domain, include_bndry_corners); + + const Real* a_dx = geom[lev].CellSize(); + Real influx = 0.0, outflux = 0.0; + compute_influx_outflux(lev, vels_vec, inout_masks, a_dx, influx, outflux, include_bndry_corners); + + const Real alpha_fcf = influx/outflux; // flux correction factor + correct_outflow(lev, vels_vec, inout_masks, bc_type, domain, alpha_fcf, include_bndry_corners); + + } // levels loop +} + +} diff --git a/Utils/hydro_utils.H b/Utils/hydro_utils.H index 8df7772fb..ac3699bf0 100644 --- a/Utils/hydro_utils.H +++ b/Utils/hydro_utils.H @@ -20,6 +20,17 @@ */ namespace HydroUtils { + +/** + * \brief Enforces solvablity by scaling outflow to match with inflow. + * + */ +void enforceInOutSolvability ( + const amrex::Vector>& vels_vec, + amrex::BCRec const* bc_type, + const amrex::Vector& geom, + bool include_bndry_corners = false); + /** * \brief Compute edge state and flux. Most general version for use with multilevel synchronization. * All other versions ultimately call this one.