From e97b7fa4f9b0456cefbbd1acf2994e050b45f11e Mon Sep 17 00:00:00 2001 From: Mukul Dave Date: Tue, 27 Feb 2024 13:59:59 -0800 Subject: [PATCH 01/22] adding vscode config file to gitignore --- .gitignore | 2 ++ 1 file changed, 2 insertions(+) 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* From a0a1b7a851e757dcb572f55e01cdc595ef0b8b91 Mon Sep 17 00:00:00 2001 From: Mukul Dave Date: Wed, 28 Feb 2024 18:14:54 -0800 Subject: [PATCH 02/22] initial commit for the outflow correction --- Projections/hydro_MacProjector.H | 5 +++ Projections/hydro_MacProjector.cpp | 66 ++++++++++++++++++++++++++++++ 2 files changed, 71 insertions(+) diff --git a/Projections/hydro_MacProjector.H b/Projections/hydro_MacProjector.H index df6c7bf85..72b8f3d6f 100644 --- a/Projections/hydro_MacProjector.H +++ b/Projections/hydro_MacProjector.H @@ -81,6 +81,11 @@ public: const amrex::Vector& a_divu = {}); #endif + void enforceSolvability ( + const amrex::Vector>& a_umac, + amrex::BCRec const* bc_type, + const amrex::Vector& geom); + /** Initialize the underlying linear operator and MLMG instances */ void initProjector ( diff --git a/Projections/hydro_MacProjector.cpp b/Projections/hydro_MacProjector.cpp index 8dd8fafc3..62422a048 100644 --- a/Projections/hydro_MacProjector.cpp +++ b/Projections/hydro_MacProjector.cpp @@ -4,6 +4,7 @@ #include #include +#include #include @@ -50,6 +51,71 @@ MacProjector::MacProjector (const Vector >& a_um setDivU(a_divu); } +void MacProjector::enforceSolvability ( + const Vector>& a_umac, + const BCRec* bc_type, + const Vector& geom +) +{ + // get the level zero domain + const Box domain = geom[0].Domain(); + + int lev = 0; // change this into a loop later ******** + //for (int lev = 0; lev < m_repo.num_active_levels(); ++lev) { + + const Real* a_dx = geom[lev].CellSize(); + + Real influx, outflux; + + // loop over the six orientations + for (OrientationIter oit; oit != nullptr; ++oit) { + const auto ori = oit(); + const auto side = ori.faceDir(); + const int dir = ori.coordDir(); + + // domain extent indices for the mac velocities + const int dlo = domain.smallEnd(dir); + const int dhi = domain.bigEnd(dir) + 1; // because face-centered + + // get BCs for the normal velocity and set the boundary index + const BCRec ibcrec = bc_type[dir]; + if (side == Orientation::low) { + const int bc = ibcrec.lo(dir); + const int bndry = dlo; + } else { + const int bc = ibcrec.hi(dir); + const int bndry = dhi; + } + + + if (bc == BCType::user_1) { + + // normal face area + const Real ds = a_dx[(dir+1) % AMREX_SPACEDIM] * a_dx[(dir+2) % AMREX_SPACEDIM]; + // Multifab reference for normal mac velocity + auto& mac_vel_mf = a_umac[lev][dir]; + + for (MFIter mfi(*mac_vel_mf, TilingIfNotGPU()); mfi.isValid(); ++mfi) { + + const Box& box = mfi.validbox(); + // create a 2D box normal to dir at the low/high boundary + const Box box2d = box.setRange(dir, bndry); + + auto mac_vel = mac_vel_mf->array(mfi); + + ParallelFor(box2d, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + a(i,j,k) += b(i,j,k) * c(i,j,k); + }); + + + + + } + } + } +} + void MacProjector::initProjector ( const LPInfo& a_lpinfo, const Vector >& a_beta, From 510711528c873c2dac4e982aa3bd8b281936efb2 Mon Sep 17 00:00:00 2001 From: Mukul Dave Date: Thu, 29 Feb 2024 16:56:14 -0800 Subject: [PATCH 03/22] function for setting the mask values for flux calculation at the boundaries --- Projections/hydro_MacProjector.cpp | 139 +++++++++++++++++++---------- 1 file changed, 90 insertions(+), 49 deletions(-) diff --git a/Projections/hydro_MacProjector.cpp b/Projections/hydro_MacProjector.cpp index 62422a048..fbfeaf528 100644 --- a/Projections/hydro_MacProjector.cpp +++ b/Projections/hydro_MacProjector.cpp @@ -12,6 +12,75 @@ using namespace amrex; namespace Hydro { +namespace { + +void set_masks( + const int lev, + const Vector>& a_umac, + const Vector& inflow_masks, + const Vector& outflow_masks, + const BCRec* bc_type, + const Box& domain) +{ + // loop over the six orientations + for (OrientationIter oit; oit != nullptr; ++oit) { + const auto ori = oit(); + const auto side = ori.faceDir(); + const int dir = ori.coordDir(); + + // domain extent indices for the mac velocities + const int dlo = domain.smallEnd(dir); + const int dhi = domain.bigEnd(dir) + 1; // because face-centered(?) + + // 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; + } + + // Multifab for normal mac velocity + auto& mac_vel_mf = a_umac[lev][dir]; + + // mask iMFs for the respective velocity direction + auto inflow_mask = inflow_masks[dir]; + auto outflow_mask = outflow_masks[dir]; + inflow_mask->setVal(0); outflow_mask->setVal(0); + + for (MFIter mfi(*mac_vel_mf, TilingIfNotGPU()); mfi.isValid(); ++mfi) { + + const Box box = mfi.validbox(); + // create a 2D box normal to dir at the low/high boundary + Box box2d(box); box2d.setRange(dir, bndry); + + auto mac_vel = mac_vel_mf->array(mfi); + auto in_mask = inflow_mask->array(mfi); + auto out_mask = outflow_mask->array(mfi); + + ParallelFor(box2d, [=] AMREX_GPU_DEVICE (int i, int j, int k) // need a reduction! + { + if (bc == BCType::user_1) + { + if ((side == Orientation::low && mac_vel(i,j,k) >= 0) + || (side == Orientation::high && mac_vel(i,j,k) <= 0)) { + in_mask(i,j,k) = 1; + } else { + out_mask(i,j,k) = 1; + } + } + }); + + } + } + +} + +} // file-local namespace + MacProjector::MacProjector( const Vector& a_geom, MLMG::Location a_umac_loc, @@ -63,57 +132,29 @@ void MacProjector::enforceSolvability ( int lev = 0; // change this into a loop later ******** //for (int lev = 0; lev < m_repo.num_active_levels(); ++lev) { - const Real* a_dx = geom[lev].CellSize(); + // masks to tag in/out flow at in-out boundaries + // separate iMultifab for each direction + Vector inflow_masks(AMREX_SPACEDIM); + Vector outflow_masks(AMREX_SPACEDIM); - Real influx, outflux; - - // loop over the six orientations - for (OrientationIter oit; oit != nullptr; ++oit) { - const auto ori = oit(); - const auto side = ori.faceDir(); - const int dir = ori.coordDir(); - - // domain extent indices for the mac velocities - const int dlo = domain.smallEnd(dir); - const int dhi = domain.bigEnd(dir) + 1; // because face-centered - - // get BCs for the normal velocity and set the boundary index - const BCRec ibcrec = bc_type[dir]; - if (side == Orientation::low) { - const int bc = ibcrec.lo(dir); - const int bndry = dlo; - } else { - const int bc = ibcrec.hi(dir); - const int bndry = dhi; - } - - - if (bc == BCType::user_1) { - - // normal face area - const Real ds = a_dx[(dir+1) % AMREX_SPACEDIM] * a_dx[(dir+2) % AMREX_SPACEDIM]; - // Multifab reference for normal mac velocity - auto& mac_vel_mf = a_umac[lev][dir]; - - for (MFIter mfi(*mac_vel_mf, TilingIfNotGPU()); mfi.isValid(); ++mfi) { - - const Box& box = mfi.validbox(); - // create a 2D box normal to dir at the low/high boundary - const Box box2d = box.setRange(dir, bndry); - - auto mac_vel = mac_vel_mf->array(mfi); - - ParallelFor(box2d, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - a(i,j,k) += b(i,j,k) * c(i,j,k); - }); - - - - - } - } + for (int idim = 0; idim < AMREX_SPACEDIM; idim++) + { + auto& mac_vel_mf = a_umac[lev][idim]; // normal mac velocity multifab + inflow_masks[idim]->define(mac_vel_mf->boxArray(), mac_vel_mf->DistributionMap(), 1, 0); + outflow_masks[idim]->define(mac_vel_mf->boxArray(), mac_vel_mf->DistributionMap(), 1, 0); } + set_masks(lev, a_umac, inflow_masks, outflow_masks, bc_type, domain); + + const Real* a_dx = geom[lev].CellSize(); + const Real influx = 0.0, outflux = 0.0; + // now calculate the influx and outflux separately + // compute_influx_outflux(lev, a_umac, inflow_masks, outflow_masks, a_dx); + // normal face area + //const Real ds = a_dx[(dir+1) % AMREX_SPACEDIM] * a_dx[(dir+2) % AMREX_SPACEDIM]; + + // correctionFactor + // const Real alpha = + // correct_outflow(lev, a_umac, outflow_masks, alpha); } void MacProjector::initProjector ( From 4d5c9396362471cd7e57264427a99ed8264896d8 Mon Sep 17 00:00:00 2001 From: Mukul Dave Date: Fri, 19 Apr 2024 17:11:44 -0700 Subject: [PATCH 04/22] working on the routine to set inflow and outflow masks --- Projections/hydro_MacProjector.cpp | 61 +++++++++++++++--------------- 1 file changed, 31 insertions(+), 30 deletions(-) diff --git a/Projections/hydro_MacProjector.cpp b/Projections/hydro_MacProjector.cpp index fbfeaf528..66a188e8e 100644 --- a/Projections/hydro_MacProjector.cpp +++ b/Projections/hydro_MacProjector.cpp @@ -17,8 +17,8 @@ namespace { void set_masks( const int lev, const Vector>& a_umac, - const Vector& inflow_masks, - const Vector& outflow_masks, + Array& inflow_masks, + Array& outflow_masks, const BCRec* bc_type, const Box& domain) { @@ -30,7 +30,7 @@ void set_masks( // domain extent indices for the mac velocities const int dlo = domain.smallEnd(dir); - const int dhi = domain.bigEnd(dir) + 1; // because face-centered(?) + const int dhi = domain.bigEnd(dir) + 1; // because face-centered // get BCs for the normal velocity and set the boundary index const BCRec ibcrec = bc_type[dir]; @@ -45,35 +45,34 @@ void set_masks( // Multifab for normal mac velocity auto& mac_vel_mf = a_umac[lev][dir]; - + //Print() << mac_vel_mf->boxArray() << std::endl; // mask iMFs for the respective velocity direction - auto inflow_mask = inflow_masks[dir]; - auto outflow_mask = outflow_masks[dir]; - inflow_mask->setVal(0); outflow_mask->setVal(0); - - for (MFIter mfi(*mac_vel_mf, TilingIfNotGPU()); mfi.isValid(); ++mfi) { - - const Box box = mfi.validbox(); - // create a 2D box normal to dir at the low/high boundary - Box box2d(box); box2d.setRange(dir, bndry); - - auto mac_vel = mac_vel_mf->array(mfi); - auto in_mask = inflow_mask->array(mfi); - auto out_mask = outflow_mask->array(mfi); - - ParallelFor(box2d, [=] AMREX_GPU_DEVICE (int i, int j, int k) // need a reduction! - { - if (bc == BCType::user_1) + auto& inflow_mask = inflow_masks[dir]; + auto& outflow_mask = outflow_masks[dir]; + inflow_mask.setVal(0); outflow_mask.setVal(0); + + if (bc == BCType::user_1) { + for (MFIter mfi(*mac_vel_mf, TilingIfNotGPU()); mfi.isValid(); ++mfi) { + + const Box box = mfi.validbox(); + // create a 2D box normal to dir at the low/high boundary + Box box2d(box); box2d.setRange(dir, bndry); + + auto mac_vel = mac_vel_mf->array(mfi); + auto in_mask = inflow_mask.array(mfi); + auto out_mask = outflow_mask.array(mfi); + Print() << "looping over 2d box: " << box2d << std::endl; + ParallelFor(box2d, [=] AMREX_GPU_DEVICE (int i, int j, int k) // need a reduction! { if ((side == Orientation::low && mac_vel(i,j,k) >= 0) || (side == Orientation::high && mac_vel(i,j,k) <= 0)) { + //Print() << i << " " << j << " " << k in_mask(i,j,k) = 1; } else { out_mask(i,j,k) = 1; } - } - }); - + }); + } } } @@ -131,18 +130,20 @@ void MacProjector::enforceSolvability ( int lev = 0; // change this into a loop later ******** //for (int lev = 0; lev < m_repo.num_active_levels(); ++lev) { - + Print() << "Declaring mask MF pointer vectors" << std::endl; // masks to tag in/out flow at in-out boundaries - // separate iMultifab for each direction - Vector inflow_masks(AMREX_SPACEDIM); - Vector outflow_masks(AMREX_SPACEDIM); + // separate iMultifab for each velocity direction + Array inflow_masks; + Array outflow_masks; + Print() << "Defining mask MF pointer vectors" << std::endl; for (int idim = 0; idim < AMREX_SPACEDIM; idim++) { auto& mac_vel_mf = a_umac[lev][idim]; // normal mac velocity multifab - inflow_masks[idim]->define(mac_vel_mf->boxArray(), mac_vel_mf->DistributionMap(), 1, 0); - outflow_masks[idim]->define(mac_vel_mf->boxArray(), mac_vel_mf->DistributionMap(), 1, 0); + inflow_masks[idim].define(mac_vel_mf->boxArray(), mac_vel_mf->DistributionMap(), 1, 0); + outflow_masks[idim].define(mac_vel_mf->boxArray(), mac_vel_mf->DistributionMap(), 1, 0); } + Print() << "Setting inflow and ouflow cell masks" << std::endl; set_masks(lev, a_umac, inflow_masks, outflow_masks, bc_type, domain); const Real* a_dx = geom[lev].CellSize(); From 5e83b8d1e5b017eab08820624c38350195225fa3 Mon Sep 17 00:00:00 2001 From: Mukul Dave Date: Tue, 23 Apr 2024 15:28:31 -0700 Subject: [PATCH 05/22] added the in- and out-flux computation --- Projections/hydro_MacProjector.cpp | 78 +++++++++++++++++++++++++++--- 1 file changed, 71 insertions(+), 7 deletions(-) diff --git a/Projections/hydro_MacProjector.cpp b/Projections/hydro_MacProjector.cpp index 66a188e8e..a0419e929 100644 --- a/Projections/hydro_MacProjector.cpp +++ b/Projections/hydro_MacProjector.cpp @@ -49,32 +49,96 @@ void set_masks( // mask iMFs for the respective velocity direction auto& inflow_mask = inflow_masks[dir]; auto& outflow_mask = outflow_masks[dir]; - inflow_mask.setVal(0); outflow_mask.setVal(0); if (bc == BCType::user_1) { for (MFIter mfi(*mac_vel_mf, TilingIfNotGPU()); mfi.isValid(); ++mfi) { const Box box = mfi.validbox(); +//Print() << "validbox = " << box << std::endl; // create a 2D box normal to dir at the low/high boundary Box box2d(box); box2d.setRange(dir, bndry); auto mac_vel = mac_vel_mf->array(mfi); auto in_mask = inflow_mask.array(mfi); auto out_mask = outflow_mask.array(mfi); - Print() << "looping over 2d box: " << box2d << std::endl; +//Print() << "looping over 2d box: " << box2d << std::endl; ParallelFor(box2d, [=] AMREX_GPU_DEVICE (int i, int j, int k) // need a reduction! { if ((side == Orientation::low && mac_vel(i,j,k) >= 0) || (side == Orientation::high && mac_vel(i,j,k) <= 0)) { - //Print() << i << " " << j << " " << k +//Print() << "inflow at: " << i << " " << j << " " << k +// << " mac_vel = " << mac_vel(i,j,k) << std::endl; in_mask(i,j,k) = 1; } else { +//Print() << "outflow at: " << i << " " << j << " " << k +// << " mac_vel = " << mac_vel(i,j,k) << std::endl; out_mask(i,j,k) = 1; } }); } } } +} + +void compute_influx_outflux( + const int lev, + const Vector>& a_umac, + const Array& inflow_masks, + const Array& outflow_masks, + const Real* a_dx, + Real& influx, + Real& outflux) +{ + //Array influx_a; + + // loop over the three dimensions + 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]; + //Print() << "ds is " << ds << std::endl; + // Multifab for normal mac velocity + auto& mac_vel_mf = a_umac[lev][idim]; + //Print() << mac_vel_mf->boxArray() << std::endl; + // mask iMFs for the respective velocity direction + auto& inflow_mask = inflow_masks[idim]; + auto& outflow_mask = outflow_masks[idim]; + + auto const& mac_vel_ma = mac_vel_mf->const_arrays(); + auto const& inflow_mask_ma = inflow_mask.const_arrays(); + auto const& outflow_mask_ma = outflow_mask.const_arrays(); + + influx += ds * + ParReduce(TypeList{}, + TypeList{}, + *mac_vel_mf, IntVect(0), // zero ghost cells + [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) + noexcept -> GpuTuple + { + if (inflow_mask_ma[box_no](i,j,k)) { + return { std::abs(mac_vel_ma[box_no](i,j,k)) }; + } else { + return { 0. }; + } + }); + + outflux += ds * + ParReduce(TypeList{}, + TypeList{}, + *mac_vel_mf, IntVect(0), // zero ghost cells + [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) + noexcept -> GpuTuple + { + if (outflow_mask_ma[box_no](i,j,k)) { + return { std::abs(mac_vel_ma[box_no](i,j,k)) }; + } else { + return { 0. }; + } + }); + } + Print() << "total influx is " << influx << std::endl; + Print() << "total outflux is " << outflux << std::endl; } @@ -141,17 +205,17 @@ void MacProjector::enforceSolvability ( { auto& mac_vel_mf = a_umac[lev][idim]; // normal mac velocity multifab inflow_masks[idim].define(mac_vel_mf->boxArray(), mac_vel_mf->DistributionMap(), 1, 0); + inflow_masks[idim].setVal(0); outflow_masks[idim].define(mac_vel_mf->boxArray(), mac_vel_mf->DistributionMap(), 1, 0); + outflow_masks[idim].setVal(0); } Print() << "Setting inflow and ouflow cell masks" << std::endl; set_masks(lev, a_umac, inflow_masks, outflow_masks, bc_type, domain); const Real* a_dx = geom[lev].CellSize(); - const Real influx = 0.0, outflux = 0.0; + Real influx = 0.0, outflux = 0.0; // now calculate the influx and outflux separately - // compute_influx_outflux(lev, a_umac, inflow_masks, outflow_masks, a_dx); - // normal face area - //const Real ds = a_dx[(dir+1) % AMREX_SPACEDIM] * a_dx[(dir+2) % AMREX_SPACEDIM]; + compute_influx_outflux(lev, a_umac, inflow_masks, outflow_masks, a_dx, influx, outflux); // correctionFactor // const Real alpha = From 0d64ad3d82302a666c33b305d0dc60c99b72d5b7 Mon Sep 17 00:00:00 2001 From: Mukul Dave Date: Tue, 23 Apr 2024 18:23:22 -0700 Subject: [PATCH 06/22] added the outflow velocity correction routine --- Projections/hydro_MacProjector.cpp | 79 ++++++++++++++++++++++++++---- 1 file changed, 70 insertions(+), 9 deletions(-) diff --git a/Projections/hydro_MacProjector.cpp b/Projections/hydro_MacProjector.cpp index a0419e929..085a4846e 100644 --- a/Projections/hydro_MacProjector.cpp +++ b/Projections/hydro_MacProjector.cpp @@ -62,7 +62,7 @@ void set_masks( auto in_mask = inflow_mask.array(mfi); auto out_mask = outflow_mask.array(mfi); //Print() << "looping over 2d box: " << box2d << std::endl; - ParallelFor(box2d, [=] AMREX_GPU_DEVICE (int i, int j, int k) // need a reduction! + ParallelFor(box2d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { if ((side == Orientation::low && mac_vel(i,j,k) >= 0) || (side == Orientation::high && mac_vel(i,j,k) <= 0)) { @@ -89,6 +89,7 @@ void compute_influx_outflux( Real& influx, Real& outflux) { + influx = 0.0, outflux = 0.0; //Array influx_a; // loop over the three dimensions @@ -137,9 +138,66 @@ void compute_influx_outflux( } }); } - Print() << "total influx is " << influx << std::endl; - Print() << "total outflux is " << outflux << std::endl; + Print() << "##### total influx is " << influx << std::endl; + Print() << "##### total outflux is " << outflux << std::endl; + // !!!!!!!!!! need to reduce these over MPI +} + +void correct_outflow( + const int lev, + const Vector>& a_umac, + const Array& outflow_masks, + const BCRec* bc_type, + const Box& domain, + const Real alpha) +{ + // loop over the six orientations + for (OrientationIter oit; oit != nullptr; ++oit) { + const auto ori = oit(); + const auto side = ori.faceDir(); + const int dir = ori.coordDir(); + + // domain extent indices for the mac velocities + const int dlo = domain.smallEnd(dir); + const int dhi = domain.bigEnd(dir) + 1; // because face-centered + + // 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; + } + // Multifab for normal mac velocity + auto& mac_vel_mf = a_umac[lev][dir]; + //Print() << mac_vel_mf->boxArray() << std::endl; + // mask iMFs for the respective velocity direction + auto& outflow_mask = outflow_masks[dir]; + + if (bc == BCType::user_1) { + for (MFIter mfi(*mac_vel_mf, TilingIfNotGPU()); mfi.isValid(); ++mfi) { + + const Box box = mfi.validbox(); +//Print() << "validbox = " << box << std::endl; + // create a 2D box normal to dir at the low/high boundary + Box box2d(box); box2d.setRange(dir, bndry); + + auto mac_vel = mac_vel_mf->array(mfi); + auto out_mask = outflow_mask.array(mfi); +//Print() << "looping over 2d box: " << box2d << std::endl; + ParallelFor(box2d, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + if (out_mask(i,j,k)) { + mac_vel(i,j,k) *= alpha; + } + }); + } + } + } } } // file-local namespace @@ -194,13 +252,12 @@ void MacProjector::enforceSolvability ( int lev = 0; // change this into a loop later ******** //for (int lev = 0; lev < m_repo.num_active_levels(); ++lev) { - Print() << "Declaring mask MF pointer vectors" << std::endl; + // masks to tag in/out flow at in-out boundaries // separate iMultifab for each velocity direction Array inflow_masks; Array outflow_masks; - Print() << "Defining mask MF pointer vectors" << std::endl; for (int idim = 0; idim < AMREX_SPACEDIM; idim++) { auto& mac_vel_mf = a_umac[lev][idim]; // normal mac velocity multifab @@ -209,7 +266,7 @@ void MacProjector::enforceSolvability ( outflow_masks[idim].define(mac_vel_mf->boxArray(), mac_vel_mf->DistributionMap(), 1, 0); outflow_masks[idim].setVal(0); } - Print() << "Setting inflow and ouflow cell masks" << std::endl; + Print() << "##### Setting inflow and ouflow cell masks" << std::endl; set_masks(lev, a_umac, inflow_masks, outflow_masks, bc_type, domain); const Real* a_dx = geom[lev].CellSize(); @@ -217,9 +274,13 @@ void MacProjector::enforceSolvability ( // now calculate the influx and outflux separately compute_influx_outflux(lev, a_umac, inflow_masks, outflow_masks, a_dx, influx, outflux); - // correctionFactor - // const Real alpha = - // correct_outflow(lev, a_umac, outflow_masks, alpha); + // apply correction factor to outflow + Print() << "##### Correcting outflow to match with inflow" << std::endl; + const Real alpha = influx/outflux; + correct_outflow(lev, a_umac, outflow_masks, bc_type, domain, alpha); + + // verify flux balance + compute_influx_outflux(lev, a_umac, inflow_masks, outflow_masks, a_dx, influx, outflux); } void MacProjector::initProjector ( From c43f0757f9b63cfd503f57e0ed7a13590ee9520a Mon Sep 17 00:00:00 2001 From: Mukul Dave Date: Thu, 6 Jun 2024 11:49:31 -0700 Subject: [PATCH 07/22] replaced user_1 with new direction_dependent BCType in amrex --- Projections/hydro_MacProjector.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Projections/hydro_MacProjector.cpp b/Projections/hydro_MacProjector.cpp index 085a4846e..b8304d95d 100644 --- a/Projections/hydro_MacProjector.cpp +++ b/Projections/hydro_MacProjector.cpp @@ -50,7 +50,7 @@ void set_masks( auto& inflow_mask = inflow_masks[dir]; auto& outflow_mask = outflow_masks[dir]; - if (bc == BCType::user_1) { + if (bc == BCType::direction_dependent) { for (MFIter mfi(*mac_vel_mf, TilingIfNotGPU()); mfi.isValid(); ++mfi) { const Box box = mfi.validbox(); @@ -178,7 +178,7 @@ void correct_outflow( // mask iMFs for the respective velocity direction auto& outflow_mask = outflow_masks[dir]; - if (bc == BCType::user_1) { + if (bc == BCType::direction_dependent) { for (MFIter mfi(*mac_vel_mf, TilingIfNotGPU()); mfi.isValid(); ++mfi) { const Box box = mfi.validbox(); From a6721078ac89f15f57bcb87854b6e498aa35d027 Mon Sep 17 00:00:00 2001 From: Mukul Dave Date: Thu, 6 Jun 2024 16:26:21 -0700 Subject: [PATCH 08/22] put in MPI reduction and corrected a bug in set_masks --- Projections/hydro_MacProjector.cpp | 51 ++++++++++++++++++++---------- 1 file changed, 35 insertions(+), 16 deletions(-) diff --git a/Projections/hydro_MacProjector.cpp b/Projections/hydro_MacProjector.cpp index b8304d95d..c35ffded0 100644 --- a/Projections/hydro_MacProjector.cpp +++ b/Projections/hydro_MacProjector.cpp @@ -27,6 +27,8 @@ void set_masks( const auto ori = oit(); const auto side = ori.faceDir(); const int dir = ori.coordDir(); + const auto islow = ori.isLow(); + const auto ishigh = ori.isHigh(); // domain extent indices for the mac velocities const int dlo = domain.smallEnd(dir); @@ -55,26 +57,32 @@ void set_masks( const Box box = mfi.validbox(); //Print() << "validbox = " << box << std::endl; - // create a 2D box normal to dir at the low/high boundary - Box box2d(box); box2d.setRange(dir, bndry); - auto mac_vel = mac_vel_mf->array(mfi); - auto in_mask = inflow_mask.array(mfi); - auto out_mask = outflow_mask.array(mfi); + // 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 mac_vel = mac_vel_mf->array(mfi); + auto in_mask = inflow_mask.array(mfi); + auto out_mask = outflow_mask.array(mfi); //Print() << "looping over 2d box: " << box2d << std::endl; - ParallelFor(box2d, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - if ((side == Orientation::low && mac_vel(i,j,k) >= 0) - || (side == Orientation::high && mac_vel(i,j,k) <= 0)) { + ParallelFor(box2d, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + if ((side == Orientation::low && mac_vel(i,j,k) >= 0) + || (side == Orientation::high && mac_vel(i,j,k) <= 0)) { //Print() << "inflow at: " << i << " " << j << " " << k // << " mac_vel = " << mac_vel(i,j,k) << std::endl; - in_mask(i,j,k) = 1; - } else { + in_mask(i,j,k) = 1; + } else { //Print() << "outflow at: " << i << " " << j << " " << k // << " mac_vel = " << mac_vel(i,j,k) << std::endl; - out_mask(i,j,k) = 1; - } - }); + out_mask(i,j,k) = 1; + } + }); + } } } } @@ -138,9 +146,17 @@ void compute_influx_outflux( } }); } - Print() << "##### total influx is " << influx << std::endl; - Print() << "##### total outflux is " << outflux << std::endl; + AllPrint() << "##### influx on rank " << ParallelDescriptor::MyProc() << " is " << influx << std::endl; + AllPrint() << "##### outflux on rank " << ParallelDescriptor::MyProc() << " is " << outflux << std::endl; // !!!!!!!!!! need to reduce these over MPI + ParallelDescriptor::ReduceRealSum(influx); + ParallelDescriptor::ReduceRealSum(outflux); + //ParallelDescriptor::Barrier(); + AllPrint() << "##### total influx on rank " << ParallelDescriptor::MyProc() << " is " << influx << std::endl; + AllPrint() << "##### total outflux on rank " << ParallelDescriptor::MyProc() << " is " << outflux << std::endl; + //Print() << "##### total influx is " << influx << std::endl; + //Print() << "##### total outflux is " << outflux << std::endl; + } void correct_outflow( @@ -274,6 +290,9 @@ void MacProjector::enforceSolvability ( // now calculate the influx and outflux separately compute_influx_outflux(lev, a_umac, inflow_masks, outflow_masks, a_dx, influx, outflux); + //ParallelDescriptor::Barrier(); + //AllPrint() << "!!!!! total influx on rank " << ParallelDescriptor::MyProc() << " is " << influx << std::endl; + //AllPrint() << "!!!!! total outflux on rank " << ParallelDescriptor::MyProc() << " is " << outflux << std::endl; // apply correction factor to outflow Print() << "##### Correcting outflow to match with inflow" << std::endl; const Real alpha = influx/outflux; From b8ea32b7a7de636a4c6f7a4ce49408ffb2d0a6a8 Mon Sep 17 00:00:00 2001 From: Mukul Dave Date: Fri, 7 Jun 2024 10:20:31 -0700 Subject: [PATCH 09/22] clean up for draft PR --- Projections/hydro_MacProjector.cpp | 31 +++++++++++------------------- 1 file changed, 11 insertions(+), 20 deletions(-) diff --git a/Projections/hydro_MacProjector.cpp b/Projections/hydro_MacProjector.cpp index c35ffded0..3f27a8027 100644 --- a/Projections/hydro_MacProjector.cpp +++ b/Projections/hydro_MacProjector.cpp @@ -47,11 +47,13 @@ void set_masks( // Multifab for normal mac velocity auto& mac_vel_mf = a_umac[lev][dir]; - //Print() << mac_vel_mf->boxArray() << std::endl; +//Print() << mac_vel_mf->boxArray() << std::endl; // mask iMFs for the respective velocity direction auto& inflow_mask = inflow_masks[dir]; auto& outflow_mask = outflow_masks[dir]; + // limit influx/outflux calculations to the in-out boundaries only + // needs to change later? if (bc == BCType::direction_dependent) { for (MFIter mfi(*mac_vel_mf, TilingIfNotGPU()); mfi.isValid(); ++mfi) { @@ -98,7 +100,6 @@ void compute_influx_outflux( Real& outflux) { influx = 0.0, outflux = 0.0; - //Array influx_a; // loop over the three dimensions for (int idim = 0; idim < AMREX_SPACEDIM; idim++) { @@ -106,10 +107,10 @@ void compute_influx_outflux( // normal face area const Real ds = a_dx[(idim+1) % AMREX_SPACEDIM] * a_dx[(idim+2) % AMREX_SPACEDIM]; - //Print() << "ds is " << ds << std::endl; +//Print() << "ds is " << ds << std::endl; // Multifab for normal mac velocity auto& mac_vel_mf = a_umac[lev][idim]; - //Print() << mac_vel_mf->boxArray() << std::endl; +//Print() << mac_vel_mf->boxArray() << std::endl; // mask iMFs for the respective velocity direction auto& inflow_mask = inflow_masks[idim]; auto& outflow_mask = outflow_masks[idim]; @@ -135,7 +136,7 @@ void compute_influx_outflux( outflux += ds * ParReduce(TypeList{}, TypeList{}, - *mac_vel_mf, IntVect(0), // zero ghost cells + *mac_vel_mf, IntVect(0), [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept -> GpuTuple { @@ -146,17 +147,11 @@ void compute_influx_outflux( } }); } - AllPrint() << "##### influx on rank " << ParallelDescriptor::MyProc() << " is " << influx << std::endl; - AllPrint() << "##### outflux on rank " << ParallelDescriptor::MyProc() << " is " << outflux << std::endl; - // !!!!!!!!!! need to reduce these over MPI + ParallelDescriptor::ReduceRealSum(influx); ParallelDescriptor::ReduceRealSum(outflux); - //ParallelDescriptor::Barrier(); - AllPrint() << "##### total influx on rank " << ParallelDescriptor::MyProc() << " is " << influx << std::endl; - AllPrint() << "##### total outflux on rank " << ParallelDescriptor::MyProc() << " is " << outflux << std::endl; - //Print() << "##### total influx is " << influx << std::endl; - //Print() << "##### total outflux is " << outflux << std::endl; - +Print() << "##### total influx is " << influx << std::endl; +Print() << "##### total outflux is " << outflux << std::endl; } void correct_outflow( @@ -190,7 +185,7 @@ void correct_outflow( // Multifab for normal mac velocity auto& mac_vel_mf = a_umac[lev][dir]; - //Print() << mac_vel_mf->boxArray() << std::endl; +//Print() << mac_vel_mf->boxArray() << std::endl; // mask iMFs for the respective velocity direction auto& outflow_mask = outflow_masks[dir]; @@ -282,7 +277,6 @@ void MacProjector::enforceSolvability ( outflow_masks[idim].define(mac_vel_mf->boxArray(), mac_vel_mf->DistributionMap(), 1, 0); outflow_masks[idim].setVal(0); } - Print() << "##### Setting inflow and ouflow cell masks" << std::endl; set_masks(lev, a_umac, inflow_masks, outflow_masks, bc_type, domain); const Real* a_dx = geom[lev].CellSize(); @@ -290,11 +284,8 @@ void MacProjector::enforceSolvability ( // now calculate the influx and outflux separately compute_influx_outflux(lev, a_umac, inflow_masks, outflow_masks, a_dx, influx, outflux); - //ParallelDescriptor::Barrier(); - //AllPrint() << "!!!!! total influx on rank " << ParallelDescriptor::MyProc() << " is " << influx << std::endl; - //AllPrint() << "!!!!! total outflux on rank " << ParallelDescriptor::MyProc() << " is " << outflux << std::endl; // apply correction factor to outflow - Print() << "##### Correcting outflow to match with inflow" << std::endl; +Print() << "##### Correcting outflow to match with inflow" << std::endl; const Real alpha = influx/outflux; correct_outflow(lev, a_umac, outflow_masks, bc_type, domain, alpha); From 3115f03323e965477286b0ec1b6574c981e80d52 Mon Sep 17 00:00:00 2001 From: Mukul Dave Date: Sun, 23 Jun 2024 08:32:18 -0700 Subject: [PATCH 10/22] Added ability to enforce inflow-outflow solvability for cell-centered velocities And moved the functions to their own file under Utils --- Projections/hydro_MacProjector.H | 5 - Projections/hydro_MacProjector.cpp | 242 ----------------- Utils/CMakeLists.txt | 1 + Utils/hydro_enforce_inout_solvability.cpp | 312 ++++++++++++++++++++++ Utils/hydro_utils.H | 7 + 5 files changed, 320 insertions(+), 247 deletions(-) create mode 100644 Utils/hydro_enforce_inout_solvability.cpp diff --git a/Projections/hydro_MacProjector.H b/Projections/hydro_MacProjector.H index 72b8f3d6f..df6c7bf85 100644 --- a/Projections/hydro_MacProjector.H +++ b/Projections/hydro_MacProjector.H @@ -81,11 +81,6 @@ public: const amrex::Vector& a_divu = {}); #endif - void enforceSolvability ( - const amrex::Vector>& a_umac, - amrex::BCRec const* bc_type, - const amrex::Vector& geom); - /** Initialize the underlying linear operator and MLMG instances */ void initProjector ( diff --git a/Projections/hydro_MacProjector.cpp b/Projections/hydro_MacProjector.cpp index 3f27a8027..103fb1c9f 100644 --- a/Projections/hydro_MacProjector.cpp +++ b/Projections/hydro_MacProjector.cpp @@ -12,207 +12,6 @@ using namespace amrex; namespace Hydro { -namespace { - -void set_masks( - const int lev, - const Vector>& a_umac, - Array& inflow_masks, - Array& outflow_masks, - const BCRec* bc_type, - const Box& domain) -{ - // loop over the six orientations - 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(); - - // domain extent indices for the mac velocities - const int dlo = domain.smallEnd(dir); - const int dhi = domain.bigEnd(dir) + 1; // because face-centered - - // 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; - } - - // Multifab for normal mac velocity - auto& mac_vel_mf = a_umac[lev][dir]; -//Print() << mac_vel_mf->boxArray() << std::endl; - // mask iMFs for the respective velocity direction - auto& inflow_mask = inflow_masks[dir]; - auto& outflow_mask = outflow_masks[dir]; - - // limit influx/outflux calculations to the in-out boundaries only - // needs to change later? - if (bc == BCType::direction_dependent) { - for (MFIter mfi(*mac_vel_mf, TilingIfNotGPU()); mfi.isValid(); ++mfi) { - - const Box box = mfi.validbox(); -//Print() << "validbox = " << box << std::endl; - - // 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 mac_vel = mac_vel_mf->array(mfi); - auto in_mask = inflow_mask.array(mfi); - auto out_mask = outflow_mask.array(mfi); -//Print() << "looping over 2d box: " << box2d << std::endl; - ParallelFor(box2d, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - if ((side == Orientation::low && mac_vel(i,j,k) >= 0) - || (side == Orientation::high && mac_vel(i,j,k) <= 0)) { -//Print() << "inflow at: " << i << " " << j << " " << k -// << " mac_vel = " << mac_vel(i,j,k) << std::endl; - in_mask(i,j,k) = 1; - } else { -//Print() << "outflow at: " << i << " " << j << " " << k -// << " mac_vel = " << mac_vel(i,j,k) << std::endl; - out_mask(i,j,k) = 1; - } - }); - } - } - } - } -} - -void compute_influx_outflux( - const int lev, - const Vector>& a_umac, - const Array& inflow_masks, - const Array& outflow_masks, - const Real* a_dx, - Real& influx, - Real& outflux) -{ - influx = 0.0, outflux = 0.0; - - // loop over the three dimensions - 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]; -//Print() << "ds is " << ds << std::endl; - // Multifab for normal mac velocity - auto& mac_vel_mf = a_umac[lev][idim]; -//Print() << mac_vel_mf->boxArray() << std::endl; - // mask iMFs for the respective velocity direction - auto& inflow_mask = inflow_masks[idim]; - auto& outflow_mask = outflow_masks[idim]; - - auto const& mac_vel_ma = mac_vel_mf->const_arrays(); - auto const& inflow_mask_ma = inflow_mask.const_arrays(); - auto const& outflow_mask_ma = outflow_mask.const_arrays(); - - influx += ds * - ParReduce(TypeList{}, - TypeList{}, - *mac_vel_mf, IntVect(0), // zero ghost cells - [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) - noexcept -> GpuTuple - { - if (inflow_mask_ma[box_no](i,j,k)) { - return { std::abs(mac_vel_ma[box_no](i,j,k)) }; - } else { - return { 0. }; - } - }); - - outflux += ds * - ParReduce(TypeList{}, - TypeList{}, - *mac_vel_mf, IntVect(0), - [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) - noexcept -> GpuTuple - { - if (outflow_mask_ma[box_no](i,j,k)) { - return { std::abs(mac_vel_ma[box_no](i,j,k)) }; - } else { - return { 0. }; - } - }); - } - - ParallelDescriptor::ReduceRealSum(influx); - ParallelDescriptor::ReduceRealSum(outflux); -Print() << "##### total influx is " << influx << std::endl; -Print() << "##### total outflux is " << outflux << std::endl; -} - -void correct_outflow( - const int lev, - const Vector>& a_umac, - const Array& outflow_masks, - const BCRec* bc_type, - const Box& domain, - const Real alpha) -{ - // loop over the six orientations - for (OrientationIter oit; oit != nullptr; ++oit) { - const auto ori = oit(); - const auto side = ori.faceDir(); - const int dir = ori.coordDir(); - - // domain extent indices for the mac velocities - const int dlo = domain.smallEnd(dir); - const int dhi = domain.bigEnd(dir) + 1; // because face-centered - - // 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; - } - - // Multifab for normal mac velocity - auto& mac_vel_mf = a_umac[lev][dir]; -//Print() << mac_vel_mf->boxArray() << std::endl; - // mask iMFs for the respective velocity direction - auto& outflow_mask = outflow_masks[dir]; - - if (bc == BCType::direction_dependent) { - for (MFIter mfi(*mac_vel_mf, TilingIfNotGPU()); mfi.isValid(); ++mfi) { - - const Box box = mfi.validbox(); -//Print() << "validbox = " << box << std::endl; - // create a 2D box normal to dir at the low/high boundary - Box box2d(box); box2d.setRange(dir, bndry); - - auto mac_vel = mac_vel_mf->array(mfi); - auto out_mask = outflow_mask.array(mfi); -//Print() << "looping over 2d box: " << box2d << std::endl; - ParallelFor(box2d, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - if (out_mask(i,j,k)) { - mac_vel(i,j,k) *= alpha; - } - }); - } - } - } -} - -} // file-local namespace - MacProjector::MacProjector( const Vector& a_geom, MLMG::Location a_umac_loc, @@ -252,47 +51,6 @@ MacProjector::MacProjector (const Vector >& a_um setDivU(a_divu); } -void MacProjector::enforceSolvability ( - const Vector>& a_umac, - const BCRec* bc_type, - const Vector& geom -) -{ - // get the level zero domain - const Box domain = geom[0].Domain(); - - int lev = 0; // change this into a loop later ******** - //for (int lev = 0; lev < m_repo.num_active_levels(); ++lev) { - - // masks to tag in/out flow at in-out boundaries - // separate iMultifab for each velocity direction - Array inflow_masks; - Array outflow_masks; - - for (int idim = 0; idim < AMREX_SPACEDIM; idim++) - { - auto& mac_vel_mf = a_umac[lev][idim]; // normal mac velocity multifab - inflow_masks[idim].define(mac_vel_mf->boxArray(), mac_vel_mf->DistributionMap(), 1, 0); - inflow_masks[idim].setVal(0); - outflow_masks[idim].define(mac_vel_mf->boxArray(), mac_vel_mf->DistributionMap(), 1, 0); - outflow_masks[idim].setVal(0); - } - set_masks(lev, a_umac, inflow_masks, outflow_masks, bc_type, domain); - - const Real* a_dx = geom[lev].CellSize(); - Real influx = 0.0, outflux = 0.0; - // now calculate the influx and outflux separately - compute_influx_outflux(lev, a_umac, inflow_masks, outflow_masks, a_dx, influx, outflux); - - // apply correction factor to outflow -Print() << "##### Correcting outflow to match with inflow" << std::endl; - const Real alpha = influx/outflux; - correct_outflow(lev, a_umac, outflow_masks, bc_type, domain, alpha); - - // verify flux balance - compute_influx_outflux(lev, a_umac, inflow_masks, outflow_masks, a_dx, influx, outflux); -} - void MacProjector::initProjector ( const LPInfo& a_lpinfo, const Vector >& a_beta, 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/hydro_enforce_inout_solvability.cpp b/Utils/hydro_enforce_inout_solvability.cpp new file mode 100644 index 000000000..aed5c58f4 --- /dev/null +++ b/Utils/hydro_enforce_inout_solvability.cpp @@ -0,0 +1,312 @@ +/** \addtogroup Utilities + * @{ + */ + +#include + +using namespace amrex; + +namespace HydroUtils { + +namespace { + +void set_inout_masks( + const int lev, + const Vector>& a_umac, + Array& inflow_masks, + Array& outflow_masks, + const BCRec* bc_type, + const Box& domain, + const bool corners) +{ + // loop over the six orientations + 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 mac velocity + auto& mac_vel_mf = a_umac[lev][dir]; +//Print() << mac_vel_mf->boxArray() << std::endl; + // mask iMFs for the respective velocity direction + auto& inflow_mask = inflow_masks[dir]; + auto& outflow_mask = outflow_masks[dir]; + + // domain extent indices for the velocities + IndexType::CellIndex dir_index_type = (mac_vel_mf->ixType()).ixType(dir); + 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; + } + + // limit influx/outflux calculations to the in-out boundaries only + // needs to change later? + if (bc == BCType::direction_dependent) { + for (MFIter mfi(*mac_vel_mf, TilingIfNotGPU()); mfi.isValid(); ++mfi) { + + Box box = mfi.validbox(); +//Print() << "validbox = " << box << std::endl; + 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 mac_vel = mac_vel_mf->array(mfi); + auto in_mask = inflow_mask.array(mfi); + auto out_mask = outflow_mask.array(mfi); +Print() << "looping over 2d box: " << box2d << std::endl; + ParallelFor(box2d, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + if ((side == Orientation::low && mac_vel(i,j,k) >= 0) + || (side == Orientation::high && mac_vel(i,j,k) <= 0)) { +//Print() << "inflow at: " << i << " " << j << " " << k +// << " mac_vel = " << mac_vel(i,j,k) << std::endl; + in_mask(i,j,k) = 1; + } else { +//Print() << "outflow at: " << i << " " << j << " " << k +// << " mac_vel = " << mac_vel(i,j,k) << std::endl; + out_mask(i,j,k) = 1; + } + }); + } + } + } + } +} + +void compute_influx_outflux( + const int lev, + const Vector>& a_umac, + const Array& inflow_masks, + const Array& outflow_masks, + const Real* a_dx, + Real& influx, + Real& outflux, + const bool corners) +{ + influx = 0.0, outflux = 0.0; + + // loop over the three dimensions + 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]; +//Print() << "ds is " << ds << std::endl; + // Multifab for normal mac velocity + auto& mac_vel_mf = a_umac[lev][idim]; +//Print() << mac_vel_mf->boxArray() << std::endl; + + IndexType index_type = mac_vel_mf->ixType(); + index_type.flip(idim); IntVect ngrow = index_type.ixType(); + if (corners) { + ngrow[(idim+1)%AMREX_SPACEDIM] = 1; + ngrow[(idim+2)%AMREX_SPACEDIM] = 1; + } + + // mask iMFs for the respective velocity direction + auto& inflow_mask = inflow_masks[idim]; + auto& outflow_mask = outflow_masks[idim]; + + auto const& mac_vel_ma = mac_vel_mf->const_arrays(); + auto const& inflow_mask_ma = inflow_mask.const_arrays(); + auto const& outflow_mask_ma = outflow_mask.const_arrays(); + + influx += ds * + ParReduce(TypeList{}, + TypeList{}, + *mac_vel_mf, ngrow, + [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) + noexcept -> GpuTuple + { + if (inflow_mask_ma[box_no](i,j,k)) { +//Print() << "counting inflow at: "<< i << " " << j << " " << k +// << " mac_vel = " << mac_vel_ma[box_no](i,j,k) << std::endl; + return { std::abs(mac_vel_ma[box_no](i,j,k)) }; + } else { + return { 0. }; + } + }); + + outflux += ds * + ParReduce(TypeList{}, + TypeList{}, + *mac_vel_mf, ngrow, + [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) + noexcept -> GpuTuple + { + if (outflow_mask_ma[box_no](i,j,k)) { +//Print() << "counting outflow at: "<< i << " " << j << " " << k +// << " mac_vel = " << mac_vel_ma[box_no](i,j,k) << std::endl; + return { std::abs(mac_vel_ma[box_no](i,j,k)) }; + } else { + return { 0. }; + } + }); + } + + ParallelDescriptor::ReduceRealSum(influx); + ParallelDescriptor::ReduceRealSum(outflux); +Print() << "##### total influx is " << influx << std::endl; +Print() << "##### total outflux is " << outflux << std::endl; +} + +void correct_outflow( + const int lev, + const Vector>& a_umac, + const Array& outflow_masks, + const BCRec* bc_type, + const Box& domain, + const Real alpha, + const bool corners) +{ + // loop over the six orientations + 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 mac velocity + auto& mac_vel_mf = a_umac[lev][dir]; +//Print() << mac_vel_mf->boxArray() << std::endl; + // mask iMFs for the respective velocity direction + auto& outflow_mask = outflow_masks[dir]; + + // domain extent indices for the velocities + IndexType::CellIndex dir_index_type = (mac_vel_mf->ixType()).ixType(dir); + 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(*mac_vel_mf, TilingIfNotGPU()); mfi.isValid(); ++mfi) { + + Box box = mfi.validbox(); +//Print() << "validbox = " << box << std::endl; + 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 mac_vel = mac_vel_mf->array(mfi); + auto out_mask = outflow_mask.array(mfi); +//Print() << "looping over 2d box: " << box2d << std::endl; + ParallelFor(box2d, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + if (out_mask(i,j,k)) { + mac_vel(i,j,k) *= alpha; + } + }); + } + } + } + } +} + +} // file-local namespace + + +void enforceInOutSolvability ( + const Vector>& a_umac, + const BCRec* bc_type, + const Vector& geom, + const bool include_bndry_corners +) +{ + // get the level zero domain + const Box domain = geom[0].Domain(); + + int lev = 0; // change this into a loop later ******** + //for (int lev = 0; lev < m_repo.num_active_levels(); ++lev) { + + // masks to tag in/out flow at in-out boundaries + // separate iMultifab for each velocity direction + Array inflow_masks; + Array outflow_masks; + + for (int idim = 0; idim < AMREX_SPACEDIM; idim++) + { + auto& mac_vel_mf = a_umac[lev][idim]; // normal velocity multifab + + IndexType index_type = mac_vel_mf->ixType(); + index_type.flip(idim); IntVect ngrow = index_type.ixType(); + if (include_bndry_corners) { + ngrow[(idim+1)%AMREX_SPACEDIM] = 1; + ngrow[(idim+2)%AMREX_SPACEDIM] = 1; + } + + inflow_masks[idim].define(mac_vel_mf->boxArray(), mac_vel_mf->DistributionMap(), 1, ngrow); + inflow_masks[idim].setVal(0); + outflow_masks[idim].define(mac_vel_mf->boxArray(), mac_vel_mf->DistributionMap(), 1, ngrow); + outflow_masks[idim].setVal(0); + } + set_inout_masks(lev, a_umac, inflow_masks, outflow_masks, bc_type, domain, include_bndry_corners); + + const Real* a_dx = geom[lev].CellSize(); + Real influx = 0.0, outflux = 0.0; + // now calculate the influx and outflux separately + compute_influx_outflux(lev, a_umac, inflow_masks, outflow_masks, a_dx, influx, outflux, include_bndry_corners); + + // apply correction factor to outflow +Print() << "##### Correcting outflow to match with inflow" << std::endl; + const Real alpha = influx/outflux; + correct_outflow(lev, a_umac, outflow_masks, bc_type, domain, alpha, include_bndry_corners); + + // verify flux balance + compute_influx_outflux(lev, a_umac, inflow_masks, outflow_masks, a_dx, influx, outflux, include_bndry_corners); +} + +} \ No newline at end of file diff --git a/Utils/hydro_utils.H b/Utils/hydro_utils.H index f0e1a6652..955359563 100644 --- a/Utils/hydro_utils.H +++ b/Utils/hydro_utils.H @@ -20,6 +20,13 @@ */ namespace HydroUtils { + +void enforceInOutSolvability ( + const amrex::Vector>& a_umac, + amrex::BCRec const* bc_type, + const amrex::Vector& geom, + const bool include_bndry_corners = false); + /** * \brief Compute edge state and flux. Most general version for use with multilevel synchonization. * All other versions ultimately call this one. From 9cb268a5d3548d33c524cb6edf4f7dc9976869fe Mon Sep 17 00:00:00 2001 From: Mukul Dave Date: Sun, 23 Jun 2024 09:23:05 -0700 Subject: [PATCH 11/22] cleanup and comments --- Utils/hydro_enforce_inout_solvability.cpp | 31 +++++++++++++++++++---- 1 file changed, 26 insertions(+), 5 deletions(-) diff --git a/Utils/hydro_enforce_inout_solvability.cpp b/Utils/hydro_enforce_inout_solvability.cpp index aed5c58f4..c2c0d6653 100644 --- a/Utils/hydro_enforce_inout_solvability.cpp +++ b/Utils/hydro_enforce_inout_solvability.cpp @@ -38,13 +38,16 @@ void set_inout_masks( IndexType::CellIndex dir_index_type = (mac_vel_mf->ixType()).ixType(dir); int dlo; if (dir_index_type == IndexType::CellIndex::CELL) { - dlo = domain.smallEnd(dir) - 1; // cell-centered boundary + // lower boundary is at -1 for cell-centered velocity + dlo = domain.smallEnd(dir) - 1; } else { - dlo = domain.smallEnd(dir); // face-centered boundary + // 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) { @@ -62,25 +65,34 @@ void set_inout_masks( Box box = mfi.validbox(); //Print() << "validbox = " << box << std::endl; + + // include ghost cells 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 + // this is relevant for cell-centered vels only + // make this automatic based on cell-centered check ???? 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 + // 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 boundary + // create a 2D box normal to dir at the low/high bndry Box box2d(box); box2d.setRange(dir, bndry); auto mac_vel = mac_vel_mf->array(mfi); auto in_mask = inflow_mask.array(mfi); auto out_mask = outflow_mask.array(mfi); Print() << "looping over 2d box: " << box2d << std::endl; + + // 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 && mac_vel(i,j,k) >= 0) @@ -123,8 +135,12 @@ void compute_influx_outflux( auto& mac_vel_mf = a_umac[lev][idim]; //Print() << mac_vel_mf->boxArray() << std::endl; + // grow in the respective direction if vel is cell-centered IndexType index_type = mac_vel_mf->ixType(); index_type.flip(idim); IntVect ngrow = index_type.ixType(); + + // grow in the transverse direction to include boundary corners + // make this automatic based on cell-centered check ???? if (corners) { ngrow[(idim+1)%AMREX_SPACEDIM] = 1; ngrow[(idim+2)%AMREX_SPACEDIM] = 1; @@ -258,7 +274,7 @@ void correct_outflow( } // file-local namespace - +// !!!!!!! need to change mac-specific variable names void enforceInOutSolvability ( const Vector>& a_umac, const BCRec* bc_type, @@ -281,8 +297,13 @@ void enforceInOutSolvability ( { auto& mac_vel_mf = a_umac[lev][idim]; // normal velocity multifab + // grow in the respective direction if vel is cell-centered + // to include the boundary cells IndexType index_type = mac_vel_mf->ixType(); index_type.flip(idim); IntVect ngrow = index_type.ixType(); + + // grow in the transverse direction to include boundary corners + // make this automatic based on cell-centered check ???? if (include_bndry_corners) { ngrow[(idim+1)%AMREX_SPACEDIM] = 1; ngrow[(idim+2)%AMREX_SPACEDIM] = 1; From 741c157f528c88ca0be2453bda723169196d231a Mon Sep 17 00:00:00 2001 From: Mukul Dave Date: Mon, 1 Jul 2024 10:45:16 -0700 Subject: [PATCH 12/22] changing variable names to make functions generic rather than being specific to mac velocities --- Utils/hydro_enforce_inout_solvability.cpp | 56 +++++++++++------------ Utils/hydro_utils.H | 6 ++- 2 files changed, 33 insertions(+), 29 deletions(-) diff --git a/Utils/hydro_enforce_inout_solvability.cpp b/Utils/hydro_enforce_inout_solvability.cpp index c2c0d6653..0a4a99a74 100644 --- a/Utils/hydro_enforce_inout_solvability.cpp +++ b/Utils/hydro_enforce_inout_solvability.cpp @@ -12,7 +12,7 @@ namespace { void set_inout_masks( const int lev, - const Vector>& a_umac, + const Vector>& vels_vec, Array& inflow_masks, Array& outflow_masks, const BCRec* bc_type, @@ -28,14 +28,14 @@ void set_inout_masks( const auto ishigh = ori.isHigh(); // Multifab for normal mac velocity - auto& mac_vel_mf = a_umac[lev][dir]; -//Print() << mac_vel_mf->boxArray() << std::endl; + auto& vel_mf = vels_vec[lev][dir]; +//Print() << vel_mf->boxArray() << std::endl; // mask iMFs for the respective velocity direction auto& inflow_mask = inflow_masks[dir]; auto& outflow_mask = outflow_masks[dir]; // domain extent indices for the velocities - IndexType::CellIndex dir_index_type = (mac_vel_mf->ixType()).ixType(dir); + IndexType::CellIndex dir_index_type = (vel_mf->ixType()).ixType(dir); int dlo; if (dir_index_type == IndexType::CellIndex::CELL) { // lower boundary is at -1 for cell-centered velocity @@ -61,7 +61,7 @@ void set_inout_masks( // limit influx/outflux calculations to the in-out boundaries only // needs to change later? if (bc == BCType::direction_dependent) { - for (MFIter mfi(*mac_vel_mf, TilingIfNotGPU()); mfi.isValid(); ++mfi) { + for (MFIter mfi(*vel_mf, TilingIfNotGPU()); mfi.isValid(); ++mfi) { Box box = mfi.validbox(); //Print() << "validbox = " << box << std::endl; @@ -87,7 +87,7 @@ void set_inout_masks( // create a 2D box normal to dir at the low/high bndry Box box2d(box); box2d.setRange(dir, bndry); - auto mac_vel = mac_vel_mf->array(mfi); + auto mac_vel = vel_mf->array(mfi); auto in_mask = inflow_mask.array(mfi); auto out_mask = outflow_mask.array(mfi); Print() << "looping over 2d box: " << box2d << std::endl; @@ -114,7 +114,7 @@ Print() << "looping over 2d box: " << box2d << std::endl; void compute_influx_outflux( const int lev, - const Vector>& a_umac, + const Vector>& vels_vec, const Array& inflow_masks, const Array& outflow_masks, const Real* a_dx, @@ -132,11 +132,11 @@ void compute_influx_outflux( a_dx[(idim+1) % AMREX_SPACEDIM] * a_dx[(idim+2) % AMREX_SPACEDIM]; //Print() << "ds is " << ds << std::endl; // Multifab for normal mac velocity - auto& mac_vel_mf = a_umac[lev][idim]; -//Print() << mac_vel_mf->boxArray() << std::endl; + auto& vel_mf = vels_vec[lev][idim]; +//Print() << vel_mf->boxArray() << std::endl; // grow in the respective direction if vel is cell-centered - IndexType index_type = mac_vel_mf->ixType(); + IndexType index_type = vel_mf->ixType(); index_type.flip(idim); IntVect ngrow = index_type.ixType(); // grow in the transverse direction to include boundary corners @@ -150,14 +150,14 @@ void compute_influx_outflux( auto& inflow_mask = inflow_masks[idim]; auto& outflow_mask = outflow_masks[idim]; - auto const& mac_vel_ma = mac_vel_mf->const_arrays(); + auto const& mac_vel_ma = vel_mf->const_arrays(); auto const& inflow_mask_ma = inflow_mask.const_arrays(); auto const& outflow_mask_ma = outflow_mask.const_arrays(); influx += ds * ParReduce(TypeList{}, TypeList{}, - *mac_vel_mf, ngrow, + *vel_mf, ngrow, [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept -> GpuTuple { @@ -173,7 +173,7 @@ void compute_influx_outflux( outflux += ds * ParReduce(TypeList{}, TypeList{}, - *mac_vel_mf, ngrow, + *vel_mf, ngrow, [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept -> GpuTuple { @@ -195,7 +195,7 @@ Print() << "##### total outflux is " << outflux << std::endl; void correct_outflow( const int lev, - const Vector>& a_umac, + const Vector>& vels_vec, const Array& outflow_masks, const BCRec* bc_type, const Box& domain, @@ -211,13 +211,13 @@ void correct_outflow( const auto ishigh = ori.isHigh(); // Multifab for normal mac velocity - auto& mac_vel_mf = a_umac[lev][dir]; -//Print() << mac_vel_mf->boxArray() << std::endl; + auto& vel_mf = vels_vec[lev][dir]; +//Print() << vel_mf->boxArray() << std::endl; // mask iMFs for the respective velocity direction auto& outflow_mask = outflow_masks[dir]; // domain extent indices for the velocities - IndexType::CellIndex dir_index_type = (mac_vel_mf->ixType()).ixType(dir); + IndexType::CellIndex dir_index_type = (vel_mf->ixType()).ixType(dir); int dlo; if (dir_index_type == IndexType::CellIndex::CELL) { dlo = domain.smallEnd(dir) - 1; // cell-centered boundary @@ -238,7 +238,7 @@ void correct_outflow( } if (bc == BCType::direction_dependent) { - for (MFIter mfi(*mac_vel_mf, TilingIfNotGPU()); mfi.isValid(); ++mfi) { + for (MFIter mfi(*vel_mf, TilingIfNotGPU()); mfi.isValid(); ++mfi) { Box box = mfi.validbox(); //Print() << "validbox = " << box << std::endl; @@ -257,7 +257,7 @@ void correct_outflow( // create a 2D box normal to dir at the low/high boundary Box box2d(box); box2d.setRange(dir, bndry); - auto mac_vel = mac_vel_mf->array(mfi); + auto mac_vel = vel_mf->array(mfi); auto out_mask = outflow_mask.array(mfi); //Print() << "looping over 2d box: " << box2d << std::endl; ParallelFor(box2d, [=] AMREX_GPU_DEVICE (int i, int j, int k) @@ -276,7 +276,7 @@ void correct_outflow( // !!!!!!! need to change mac-specific variable names void enforceInOutSolvability ( - const Vector>& a_umac, + const Vector>& vels_vec, const BCRec* bc_type, const Vector& geom, const bool include_bndry_corners @@ -295,11 +295,11 @@ void enforceInOutSolvability ( for (int idim = 0; idim < AMREX_SPACEDIM; idim++) { - auto& mac_vel_mf = a_umac[lev][idim]; // normal velocity multifab + 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 = mac_vel_mf->ixType(); + IndexType index_type = vel_mf->ixType(); index_type.flip(idim); IntVect ngrow = index_type.ixType(); // grow in the transverse direction to include boundary corners @@ -309,25 +309,25 @@ void enforceInOutSolvability ( ngrow[(idim+2)%AMREX_SPACEDIM] = 1; } - inflow_masks[idim].define(mac_vel_mf->boxArray(), mac_vel_mf->DistributionMap(), 1, ngrow); + inflow_masks[idim].define(vel_mf->boxArray(), vel_mf->DistributionMap(), 1, ngrow); inflow_masks[idim].setVal(0); - outflow_masks[idim].define(mac_vel_mf->boxArray(), mac_vel_mf->DistributionMap(), 1, ngrow); + outflow_masks[idim].define(vel_mf->boxArray(), vel_mf->DistributionMap(), 1, ngrow); outflow_masks[idim].setVal(0); } - set_inout_masks(lev, a_umac, inflow_masks, outflow_masks, bc_type, domain, include_bndry_corners); + set_inout_masks(lev, vels_vec, inflow_masks, outflow_masks, bc_type, domain, include_bndry_corners); const Real* a_dx = geom[lev].CellSize(); Real influx = 0.0, outflux = 0.0; // now calculate the influx and outflux separately - compute_influx_outflux(lev, a_umac, inflow_masks, outflow_masks, a_dx, influx, outflux, include_bndry_corners); + compute_influx_outflux(lev, vels_vec, inflow_masks, outflow_masks, a_dx, influx, outflux, include_bndry_corners); // apply correction factor to outflow Print() << "##### Correcting outflow to match with inflow" << std::endl; const Real alpha = influx/outflux; - correct_outflow(lev, a_umac, outflow_masks, bc_type, domain, alpha, include_bndry_corners); + correct_outflow(lev, vels_vec, outflow_masks, bc_type, domain, alpha, include_bndry_corners); // verify flux balance - compute_influx_outflux(lev, a_umac, inflow_masks, outflow_masks, a_dx, influx, outflux, include_bndry_corners); + compute_influx_outflux(lev, vels_vec, inflow_masks, outflow_masks, a_dx, influx, outflux, include_bndry_corners); } } \ No newline at end of file diff --git a/Utils/hydro_utils.H b/Utils/hydro_utils.H index 955359563..64220db91 100644 --- a/Utils/hydro_utils.H +++ b/Utils/hydro_utils.H @@ -21,8 +21,12 @@ namespace HydroUtils { +/** + * \brief Enforces solvablity by scaling outflow to match with inflow. + * + */ void enforceInOutSolvability ( - const amrex::Vector>& a_umac, + const amrex::Vector>& vels_vec, amrex::BCRec const* bc_type, const amrex::Vector& geom, const bool include_bndry_corners = false); From 1f2c5f59a4aa80c836f47d379a00a6ef63b684e5 Mon Sep 17 00:00:00 2001 From: Mukul Dave Date: Mon, 1 Jul 2024 10:54:48 -0700 Subject: [PATCH 13/22] added the loop over levels --- Utils/hydro_enforce_inout_solvability.cpp | 77 ++++++++++++----------- 1 file changed, 39 insertions(+), 38 deletions(-) diff --git a/Utils/hydro_enforce_inout_solvability.cpp b/Utils/hydro_enforce_inout_solvability.cpp index 0a4a99a74..171c993de 100644 --- a/Utils/hydro_enforce_inout_solvability.cpp +++ b/Utils/hydro_enforce_inout_solvability.cpp @@ -285,49 +285,50 @@ void enforceInOutSolvability ( // get the level zero domain const Box domain = geom[0].Domain(); - int lev = 0; // change this into a loop later ******** - //for (int lev = 0; lev < m_repo.num_active_levels(); ++lev) { - - // masks to tag in/out flow at in-out boundaries - // separate iMultifab for each velocity direction - Array inflow_masks; - Array outflow_masks; - - for (int idim = 0; idim < AMREX_SPACEDIM; idim++) - { - 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(); + const auto nlevs = int(vels_vec.size()); + for (int lev = 0; lev < nlevs; ++lev) { + + // masks to tag in/out flow at in-out boundaries + // separate iMultifab for each velocity direction + Array inflow_masks; + Array outflow_masks; + + for (int idim = 0; idim < AMREX_SPACEDIM; idim++) + { + 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; + } - // grow in the transverse direction to include boundary corners - // make this automatic based on cell-centered check ???? - if (include_bndry_corners) { - ngrow[(idim+1)%AMREX_SPACEDIM] = 1; - ngrow[(idim+2)%AMREX_SPACEDIM] = 1; + inflow_masks[idim].define(vel_mf->boxArray(), vel_mf->DistributionMap(), 1, ngrow); + inflow_masks[idim].setVal(0); + outflow_masks[idim].define(vel_mf->boxArray(), vel_mf->DistributionMap(), 1, ngrow); + outflow_masks[idim].setVal(0); } + set_inout_masks(lev, vels_vec, inflow_masks, outflow_masks, bc_type, domain, include_bndry_corners); - inflow_masks[idim].define(vel_mf->boxArray(), vel_mf->DistributionMap(), 1, ngrow); - inflow_masks[idim].setVal(0); - outflow_masks[idim].define(vel_mf->boxArray(), vel_mf->DistributionMap(), 1, ngrow); - outflow_masks[idim].setVal(0); - } - set_inout_masks(lev, vels_vec, inflow_masks, outflow_masks, bc_type, domain, include_bndry_corners); + const Real* a_dx = geom[lev].CellSize(); + Real influx = 0.0, outflux = 0.0; + // now calculate the influx and outflux separately + compute_influx_outflux(lev, vels_vec, inflow_masks, outflow_masks, a_dx, influx, outflux, include_bndry_corners); - const Real* a_dx = geom[lev].CellSize(); - Real influx = 0.0, outflux = 0.0; - // now calculate the influx and outflux separately - compute_influx_outflux(lev, vels_vec, inflow_masks, outflow_masks, a_dx, influx, outflux, include_bndry_corners); - - // apply correction factor to outflow + // apply correction factor to outflow Print() << "##### Correcting outflow to match with inflow" << std::endl; - const Real alpha = influx/outflux; - correct_outflow(lev, vels_vec, outflow_masks, bc_type, domain, alpha, include_bndry_corners); + const Real alpha = influx/outflux; + correct_outflow(lev, vels_vec, outflow_masks, bc_type, domain, alpha, include_bndry_corners); + + // verify flux balance + compute_influx_outflux(lev, vels_vec, inflow_masks, outflow_masks, a_dx, influx, outflux, include_bndry_corners); - // verify flux balance - compute_influx_outflux(lev, vels_vec, inflow_masks, outflow_masks, a_dx, influx, outflux, include_bndry_corners); + } // levels loop } -} \ No newline at end of file +} From e2fec504ea07e84a486257168a5cf9ed3e445e9b Mon Sep 17 00:00:00 2001 From: Mukul Dave Date: Mon, 1 Jul 2024 15:02:34 -0700 Subject: [PATCH 14/22] cleaned up mac-specific variable names, comments, and debug printouts --- Utils/hydro_enforce_inout_solvability.cpp | 77 ++++++++--------------- 1 file changed, 25 insertions(+), 52 deletions(-) diff --git a/Utils/hydro_enforce_inout_solvability.cpp b/Utils/hydro_enforce_inout_solvability.cpp index 171c993de..22c629ddc 100644 --- a/Utils/hydro_enforce_inout_solvability.cpp +++ b/Utils/hydro_enforce_inout_solvability.cpp @@ -19,7 +19,6 @@ void set_inout_masks( const Box& domain, const bool corners) { - // loop over the six orientations for (OrientationIter oit; oit != nullptr; ++oit) { const auto ori = oit(); const auto side = ori.faceDir(); @@ -27,15 +26,15 @@ void set_inout_masks( const auto islow = ori.isLow(); const auto ishigh = ori.isHigh(); - // Multifab for normal mac velocity + // Multifab for normal velocity auto& vel_mf = vels_vec[lev][dir]; -//Print() << vel_mf->boxArray() << std::endl; + // mask iMFs for the respective velocity direction auto& inflow_mask = inflow_masks[dir]; auto& outflow_mask = outflow_masks[dir]; - // domain extent indices for the velocities 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 @@ -59,22 +58,18 @@ void set_inout_masks( } // limit influx/outflux calculations to the in-out boundaries only - // needs to change later? if (bc == BCType::direction_dependent) { for (MFIter mfi(*vel_mf, TilingIfNotGPU()); mfi.isValid(); ++mfi) { Box box = mfi.validbox(); -//Print() << "validbox = " << box << std::endl; - // include ghost cells for cell-centered + // 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 - // this is relevant for cell-centered vels only - // make this automatic based on cell-centered check ???? if (corners) { box.grow((dir+1)%AMREX_SPACEDIM, 1); box.grow((dir+2)%AMREX_SPACEDIM, 1); @@ -87,22 +82,17 @@ void set_inout_masks( // create a 2D box normal to dir at the low/high bndry Box box2d(box); box2d.setRange(dir, bndry); - auto mac_vel = vel_mf->array(mfi); + auto vel_arr = vel_mf->array(mfi); auto in_mask = inflow_mask.array(mfi); auto out_mask = outflow_mask.array(mfi); -Print() << "looping over 2d box: " << box2d << std::endl; // 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 && mac_vel(i,j,k) >= 0) - || (side == Orientation::high && mac_vel(i,j,k) <= 0)) { -//Print() << "inflow at: " << i << " " << j << " " << k -// << " mac_vel = " << mac_vel(i,j,k) << std::endl; + if ((side == Orientation::low && vel_arr(i,j,k) >= 0) + || (side == Orientation::high && vel_arr(i,j,k) <= 0)) { in_mask(i,j,k) = 1; } else { -//Print() << "outflow at: " << i << " " << j << " " << k -// << " mac_vel = " << mac_vel(i,j,k) << std::endl; out_mask(i,j,k) = 1; } }); @@ -124,23 +114,20 @@ void compute_influx_outflux( { influx = 0.0, outflux = 0.0; - // loop over the three dimensions 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]; -//Print() << "ds is " << ds << std::endl; - // Multifab for normal mac velocity + + // Multifab for normal velocity auto& vel_mf = vels_vec[lev][idim]; -//Print() << vel_mf->boxArray() << std::endl; // 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 - // make this automatic based on cell-centered check ???? if (corners) { ngrow[(idim+1)%AMREX_SPACEDIM] = 1; ngrow[(idim+2)%AMREX_SPACEDIM] = 1; @@ -150,7 +137,8 @@ void compute_influx_outflux( auto& inflow_mask = inflow_masks[idim]; auto& outflow_mask = outflow_masks[idim]; - auto const& mac_vel_ma = vel_mf->const_arrays(); + // define "multi-arrays" and perform reduction using masks + auto const& vel_ma = vel_mf->const_arrays(); auto const& inflow_mask_ma = inflow_mask.const_arrays(); auto const& outflow_mask_ma = outflow_mask.const_arrays(); @@ -162,9 +150,7 @@ void compute_influx_outflux( noexcept -> GpuTuple { if (inflow_mask_ma[box_no](i,j,k)) { -//Print() << "counting inflow at: "<< i << " " << j << " " << k -// << " mac_vel = " << mac_vel_ma[box_no](i,j,k) << std::endl; - return { std::abs(mac_vel_ma[box_no](i,j,k)) }; + return { std::abs(vel_ma[box_no](i,j,k)) }; } else { return { 0. }; } @@ -178,19 +164,14 @@ void compute_influx_outflux( noexcept -> GpuTuple { if (outflow_mask_ma[box_no](i,j,k)) { -//Print() << "counting outflow at: "<< i << " " << j << " " << k -// << " mac_vel = " << mac_vel_ma[box_no](i,j,k) << std::endl; - return { std::abs(mac_vel_ma[box_no](i,j,k)) }; + return { std::abs(vel_ma[box_no](i,j,k)) }; } else { return { 0. }; } }); } - ParallelDescriptor::ReduceRealSum(influx); ParallelDescriptor::ReduceRealSum(outflux); -Print() << "##### total influx is " << influx << std::endl; -Print() << "##### total outflux is " << outflux << std::endl; } void correct_outflow( @@ -199,10 +180,9 @@ void correct_outflow( const Array& outflow_masks, const BCRec* bc_type, const Box& domain, - const Real alpha, + const Real alpha_fcf, const bool corners) { - // loop over the six orientations for (OrientationIter oit; oit != nullptr; ++oit) { const auto ori = oit(); const auto side = ori.faceDir(); @@ -210,14 +190,14 @@ void correct_outflow( const auto islow = ori.isLow(); const auto ishigh = ori.isHigh(); - // Multifab for normal mac velocity + // Multifab for normal velocity auto& vel_mf = vels_vec[lev][dir]; -//Print() << vel_mf->boxArray() << std::endl; + // mask iMFs for the respective velocity direction auto& outflow_mask = outflow_masks[dir]; - // domain extent indices for the velocities 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 @@ -241,7 +221,6 @@ void correct_outflow( for (MFIter mfi(*vel_mf, TilingIfNotGPU()); mfi.isValid(); ++mfi) { Box box = mfi.validbox(); -//Print() << "validbox = " << box << std::endl; if (dir_index_type == IndexType::CellIndex::CELL) { box.grow(dir, 1); } @@ -257,13 +236,13 @@ void correct_outflow( // create a 2D box normal to dir at the low/high boundary Box box2d(box); box2d.setRange(dir, bndry); - auto mac_vel = vel_mf->array(mfi); + auto vel_arr = vel_mf->array(mfi); auto out_mask = outflow_mask.array(mfi); -//Print() << "looping over 2d box: " << box2d << std::endl; + ParallelFor(box2d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { if (out_mask(i,j,k)) { - mac_vel(i,j,k) *= alpha; + vel_arr(i,j,k) *= alpha_fcf; } }); } @@ -274,7 +253,6 @@ void correct_outflow( } // file-local namespace -// !!!!!!! need to change mac-specific variable names void enforceInOutSolvability ( const Vector>& vels_vec, const BCRec* bc_type, @@ -282,17 +260,17 @@ void enforceInOutSolvability ( const bool include_bndry_corners ) { - // get the level zero domain const Box domain = geom[0].Domain(); const auto nlevs = int(vels_vec.size()); for (int lev = 0; lev < nlevs; ++lev) { - // masks to tag in/out flow at in-out boundaries + // masks to tag inflow/outflow at the boundaries // separate iMultifab for each velocity direction Array inflow_masks; Array outflow_masks; + // defining the mask iMultifabs in each direction for (int idim = 0; idim < AMREX_SPACEDIM; idim++) { auto& vel_mf = vels_vec[lev][idim]; // normal velocity multifab @@ -313,20 +291,15 @@ void enforceInOutSolvability ( outflow_masks[idim].define(vel_mf->boxArray(), vel_mf->DistributionMap(), 1, ngrow); outflow_masks[idim].setVal(0); } + set_inout_masks(lev, vels_vec, inflow_masks, outflow_masks, bc_type, domain, include_bndry_corners); const Real* a_dx = geom[lev].CellSize(); Real influx = 0.0, outflux = 0.0; - // now calculate the influx and outflux separately compute_influx_outflux(lev, vels_vec, inflow_masks, outflow_masks, a_dx, influx, outflux, include_bndry_corners); - // apply correction factor to outflow -Print() << "##### Correcting outflow to match with inflow" << std::endl; - const Real alpha = influx/outflux; - correct_outflow(lev, vels_vec, outflow_masks, bc_type, domain, alpha, include_bndry_corners); - - // verify flux balance - compute_influx_outflux(lev, vels_vec, inflow_masks, outflow_masks, a_dx, influx, outflux, include_bndry_corners); + const Real alpha_fcf = influx/outflux; // flux correction factor + correct_outflow(lev, vels_vec, outflow_masks, bc_type, domain, alpha_fcf, include_bndry_corners); } // levels loop } From c4bd659a9795a118493634d762c13b726263bd3e Mon Sep 17 00:00:00 2001 From: Mukul Dave Date: Tue, 2 Jul 2024 10:45:05 -0700 Subject: [PATCH 15/22] fixed unnecessary include --- Projections/hydro_MacProjector.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/Projections/hydro_MacProjector.cpp b/Projections/hydro_MacProjector.cpp index 103fb1c9f..8dd8fafc3 100644 --- a/Projections/hydro_MacProjector.cpp +++ b/Projections/hydro_MacProjector.cpp @@ -4,7 +4,6 @@ #include #include -#include #include From 9ce8942bca133fef4377d6b958d0c94e3cf8f965 Mon Sep 17 00:00:00 2001 From: Mukul Dave Date: Tue, 2 Jul 2024 10:48:47 -0700 Subject: [PATCH 16/22] added new file to make.package for gnu make users --- Utils/Make.package | 1 + 1 file changed, 1 insertion(+) 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 From e0570652cdea1562a7b6ed7e599cc504b37f9c22 Mon Sep 17 00:00:00 2001 From: Mukul Dave Date: Tue, 2 Jul 2024 15:02:17 -0700 Subject: [PATCH 17/22] added docs for in-out solvability --- Docs/source/InOutSolvability.rst | 23 +++++++++++++++++++++++ Docs/source/Utilities.rst | 1 + 2 files changed, 24 insertions(+) create mode 100644 Docs/source/InOutSolvability.rst diff --git a/Docs/source/InOutSolvability.rst b/Docs/source/InOutSolvability.rst new file mode 100644 index 000000000..86f225a91 --- /dev/null +++ b/Docs/source/InOutSolvability.rst @@ -0,0 +1,23 @@ +.. 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}. 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 From 6d567d5a4fdaf5de83c1232a116c648d2e67d1af Mon Sep 17 00:00:00 2001 From: Mukul Dave Date: Tue, 2 Jul 2024 17:02:36 -0700 Subject: [PATCH 18/22] replacing separate inflow and outflow masks with a single mask --- Utils/hydro_enforce_inout_solvability.cpp | 59 ++++++++++------------- 1 file changed, 26 insertions(+), 33 deletions(-) diff --git a/Utils/hydro_enforce_inout_solvability.cpp b/Utils/hydro_enforce_inout_solvability.cpp index 22c629ddc..c56e399e1 100644 --- a/Utils/hydro_enforce_inout_solvability.cpp +++ b/Utils/hydro_enforce_inout_solvability.cpp @@ -13,8 +13,7 @@ namespace { void set_inout_masks( const int lev, const Vector>& vels_vec, - Array& inflow_masks, - Array& outflow_masks, + Array& inout_masks, const BCRec* bc_type, const Box& domain, const bool corners) @@ -29,9 +28,8 @@ void set_inout_masks( // Multifab for normal velocity auto& vel_mf = vels_vec[lev][dir]; - // mask iMFs for the respective velocity direction - auto& inflow_mask = inflow_masks[dir]; - auto& outflow_mask = outflow_masks[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 @@ -83,17 +81,16 @@ void set_inout_masks( Box box2d(box); box2d.setRange(dir, bndry); auto vel_arr = vel_mf->array(mfi); - auto in_mask = inflow_mask.array(mfi); - auto out_mask = outflow_mask.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)) { - in_mask(i,j,k) = 1; + inout_mask_arr(i,j,k) = -1; } else { - out_mask(i,j,k) = 1; + inout_mask_arr(i,j,k) = +1; } }); } @@ -105,8 +102,7 @@ void set_inout_masks( void compute_influx_outflux( const int lev, const Vector>& vels_vec, - const Array& inflow_masks, - const Array& outflow_masks, + const Array& inout_masks, const Real* a_dx, Real& influx, Real& outflux, @@ -133,14 +129,12 @@ void compute_influx_outflux( ngrow[(idim+2)%AMREX_SPACEDIM] = 1; } - // mask iMFs for the respective velocity direction - auto& inflow_mask = inflow_masks[idim]; - auto& outflow_mask = outflow_masks[idim]; + // mask iMF for the respective velocity direction + auto& inout_mask = inout_masks[idim]; - // define "multi-arrays" and perform reduction using masks + // define "multi-arrays" and perform reduction using the mask auto const& vel_ma = vel_mf->const_arrays(); - auto const& inflow_mask_ma = inflow_mask.const_arrays(); - auto const& outflow_mask_ma = outflow_mask.const_arrays(); + auto const& inout_mask_ma = inout_mask.const_arrays(); influx += ds * ParReduce(TypeList{}, @@ -149,7 +143,7 @@ void compute_influx_outflux( [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept -> GpuTuple { - if (inflow_mask_ma[box_no](i,j,k)) { + if (inout_mask_ma[box_no](i,j,k) == -1) { return { std::abs(vel_ma[box_no](i,j,k)) }; } else { return { 0. }; @@ -163,7 +157,7 @@ void compute_influx_outflux( [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept -> GpuTuple { - if (outflow_mask_ma[box_no](i,j,k)) { + if (inout_mask_ma[box_no](i,j,k) == 1) { return { std::abs(vel_ma[box_no](i,j,k)) }; } else { return { 0. }; @@ -177,7 +171,7 @@ void compute_influx_outflux( void correct_outflow( const int lev, const Vector>& vels_vec, - const Array& outflow_masks, + const Array& inout_masks, const BCRec* bc_type, const Box& domain, const Real alpha_fcf, @@ -193,8 +187,8 @@ void correct_outflow( // Multifab for normal velocity auto& vel_mf = vels_vec[lev][dir]; - // mask iMFs for the respective velocity direction - auto& outflow_mask = outflow_masks[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 @@ -237,11 +231,11 @@ void correct_outflow( Box box2d(box); box2d.setRange(dir, bndry); auto vel_arr = vel_mf->array(mfi); - auto out_mask = outflow_mask.array(mfi); + auto inout_mask_arr = inout_mask.array(mfi); ParallelFor(box2d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - if (out_mask(i,j,k)) { + if (inout_mask_arr(i,j,k) == 1) { vel_arr(i,j,k) *= alpha_fcf; } }); @@ -267,8 +261,9 @@ void enforceInOutSolvability ( // masks to tag inflow/outflow at the boundaries // separate iMultifab for each velocity direction - Array inflow_masks; - Array outflow_masks; + // 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++) @@ -286,20 +281,18 @@ void enforceInOutSolvability ( ngrow[(idim+2)%AMREX_SPACEDIM] = 1; } - inflow_masks[idim].define(vel_mf->boxArray(), vel_mf->DistributionMap(), 1, ngrow); - inflow_masks[idim].setVal(0); - outflow_masks[idim].define(vel_mf->boxArray(), vel_mf->DistributionMap(), 1, ngrow); - outflow_masks[idim].setVal(0); + inout_masks[idim].define(vel_mf->boxArray(), vel_mf->DistributionMap(), 1, ngrow); + inout_masks[idim].setVal(0); } - set_inout_masks(lev, vels_vec, inflow_masks, outflow_masks, bc_type, domain, include_bndry_corners); + 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, inflow_masks, outflow_masks, a_dx, influx, outflux, include_bndry_corners); + 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, outflow_masks, bc_type, domain, alpha_fcf, include_bndry_corners); + correct_outflow(lev, vels_vec, inout_masks, bc_type, domain, alpha_fcf, include_bndry_corners); } // levels loop } From 1a529b8ac22b2987ec1f01fe428515480bfd4d09 Mon Sep 17 00:00:00 2001 From: Mukul Dave Date: Tue, 2 Jul 2024 18:36:31 -0700 Subject: [PATCH 19/22] adds clarification to the docs regarding the function not accounting for pure inflow or outflow BCs --- Docs/source/InOutSolvability.rst | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/Docs/source/InOutSolvability.rst b/Docs/source/InOutSolvability.rst index 86f225a91..c74caa515 100644 --- a/Docs/source/InOutSolvability.rst +++ b/Docs/source/InOutSolvability.rst @@ -21,3 +21,8 @@ The new flux-conserving velocities to be used for the MAC/nodal projections, .. 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 From 3bb2f8dbbeb6e483a4d7e624ef940245eff6f8df Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Sun, 21 Jul 2024 16:22:16 -0700 Subject: [PATCH 20/22] Update hydro_utils.H --- Utils/hydro_utils.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Utils/hydro_utils.H b/Utils/hydro_utils.H index d6aaa8a7f..ac3699bf0 100644 --- a/Utils/hydro_utils.H +++ b/Utils/hydro_utils.H @@ -29,7 +29,7 @@ void enforceInOutSolvability ( const amrex::Vector>& vels_vec, amrex::BCRec const* bc_type, const amrex::Vector& geom, - const bool include_bndry_corners = false); + bool include_bndry_corners = false); /** * \brief Compute edge state and flux. Most general version for use with multilevel synchronization. From 71ba07b15cf8a3e1b99e458d7ad6dadc14581e31 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Sun, 21 Jul 2024 16:22:48 -0700 Subject: [PATCH 21/22] Update hydro_enforce_inout_solvability.cpp --- Utils/hydro_enforce_inout_solvability.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Utils/hydro_enforce_inout_solvability.cpp b/Utils/hydro_enforce_inout_solvability.cpp index c56e399e1..093fd827e 100644 --- a/Utils/hydro_enforce_inout_solvability.cpp +++ b/Utils/hydro_enforce_inout_solvability.cpp @@ -251,7 +251,7 @@ void enforceInOutSolvability ( const Vector>& vels_vec, const BCRec* bc_type, const Vector& geom, - const bool include_bndry_corners + bool include_bndry_corners ) { const Box domain = geom[0].Domain(); From b5c73ad30328c4fa4c3ed8b870ad5d67b1a16a3d Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Sun, 21 Jul 2024 16:38:02 -0700 Subject: [PATCH 22/22] Update hydro_enforce_inout_solvability.cpp --- Utils/hydro_enforce_inout_solvability.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/Utils/hydro_enforce_inout_solvability.cpp b/Utils/hydro_enforce_inout_solvability.cpp index 093fd827e..3a73f0e0c 100644 --- a/Utils/hydro_enforce_inout_solvability.cpp +++ b/Utils/hydro_enforce_inout_solvability.cpp @@ -26,7 +26,7 @@ void set_inout_masks( const auto ishigh = ori.isHigh(); // Multifab for normal velocity - auto& vel_mf = vels_vec[lev][dir]; + const auto& vel_mf = vels_vec[lev][dir]; // mask iMF for the respective velocity direction auto& inout_mask = inout_masks[dir]; @@ -117,7 +117,7 @@ void compute_influx_outflux( a_dx[(idim+1) % AMREX_SPACEDIM] * a_dx[(idim+2) % AMREX_SPACEDIM]; // Multifab for normal velocity - auto& vel_mf = vels_vec[lev][idim]; + const auto& vel_mf = vels_vec[lev][idim]; // grow in the respective direction if vel is cell-centered IndexType index_type = vel_mf->ixType(); @@ -130,7 +130,7 @@ void compute_influx_outflux( } // mask iMF for the respective velocity direction - auto& inout_mask = inout_masks[idim]; + const auto& inout_mask = inout_masks[idim]; // define "multi-arrays" and perform reduction using the mask auto const& vel_ma = vel_mf->const_arrays(); @@ -185,10 +185,10 @@ void correct_outflow( const auto ishigh = ori.isHigh(); // Multifab for normal velocity - auto& vel_mf = vels_vec[lev][dir]; + const auto& vel_mf = vels_vec[lev][dir]; // mask iMF for the respective velocity direction - auto& inout_mask = inout_masks[dir]; + const auto& inout_mask = inout_masks[dir]; IndexType::CellIndex dir_index_type = (vel_mf->ixType()).ixType(dir); // domain extent indices for the velocities @@ -268,7 +268,7 @@ void enforceInOutSolvability ( // defining the mask iMultifabs in each direction for (int idim = 0; idim < AMREX_SPACEDIM; idim++) { - auto& vel_mf = vels_vec[lev][idim]; // normal velocity multifab + 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