From 0256838bb5525d413286c798e2247ce708cb2b07 Mon Sep 17 00:00:00 2001 From: Candace Gilet Date: Thu, 14 Mar 2024 21:36:09 -0400 Subject: [PATCH] Consolidate code (#88) * Update make_eb_box() to use amrex::BoxIF. This allows specifying internal or external flow. * Remove incflo::set_velocity_bc functions in favor of using PhysBCFunct, which is needed to set the inflow bcs during FillPatch. Now we only need one function to set the inflow bc. * Rename file to better reflect contents which includes setting density and tracer EB bcs. Renamed: incflo_set_velocity_bcs.cpp -> incflo_set_bcs.cpp --- src/boundary_conditions/CMakeLists.txt | 2 +- src/boundary_conditions/Make.package | 2 +- ...et_velocity_bcs.cpp => incflo_set_bcs.cpp} | 35 ------ src/embedded_boundaries/eb_box.cpp | 55 ++-------- src/incflo.H | 3 - src/prob/CMakeLists.txt | 1 - src/prob/Make.package | 1 - src/prob/prob_bc.cpp | 101 ------------------ .../incflo_apply_nodal_projection.cpp | 24 ++++- src/projection/incflo_projection_bc.cpp | 4 +- 10 files changed, 35 insertions(+), 193 deletions(-) rename src/boundary_conditions/{incflo_set_velocity_bcs.cpp => incflo_set_bcs.cpp} (85%) delete mode 100644 src/prob/prob_bc.cpp diff --git a/src/boundary_conditions/CMakeLists.txt b/src/boundary_conditions/CMakeLists.txt index 85a6c0ea..ce9bd625 100644 --- a/src/boundary_conditions/CMakeLists.txt +++ b/src/boundary_conditions/CMakeLists.txt @@ -3,5 +3,5 @@ target_sources(incflo boundary_conditions.cpp incflo_fillpatch.cpp incflo_fillphysbc.cpp - incflo_set_velocity_bcs.cpp + incflo_set_bcs.cpp ) diff --git a/src/boundary_conditions/Make.package b/src/boundary_conditions/Make.package index f7d0b486..9de2ae89 100644 --- a/src/boundary_conditions/Make.package +++ b/src/boundary_conditions/Make.package @@ -1,4 +1,4 @@ CEXE_sources += boundary_conditions.cpp CEXE_sources += incflo_fillpatch.cpp incflo_fillphysbc.cpp -CEXE_sources += incflo_set_velocity_bcs.cpp +CEXE_sources += incflo_set_bcs.cpp diff --git a/src/boundary_conditions/incflo_set_velocity_bcs.cpp b/src/boundary_conditions/incflo_set_bcs.cpp similarity index 85% rename from src/boundary_conditions/incflo_set_velocity_bcs.cpp rename to src/boundary_conditions/incflo_set_bcs.cpp index d01a854a..a2ef9713 100644 --- a/src/boundary_conditions/incflo_set_velocity_bcs.cpp +++ b/src/boundary_conditions/incflo_set_bcs.cpp @@ -6,41 +6,6 @@ using namespace amrex; -void -incflo::set_inflow_velocity (int lev, amrex::Real time, MultiFab& vel, int nghost) -{ - Geometry const& gm = Geom(lev); - Box const& domain = gm.growPeriodicDomain(nghost); - for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { - Orientation olo(dir,Orientation::low); - Orientation ohi(dir,Orientation::high); - if (m_bc_type[olo] == BC::mass_inflow || m_bc_type[ohi] == BC::mass_inflow) { - Box dlo = (m_bc_type[olo] == BC::mass_inflow) ? amrex::adjCellLo(domain,dir,nghost) : Box(); - Box dhi = (m_bc_type[ohi] == BC::mass_inflow) ? amrex::adjCellHi(domain,dir,nghost) : Box(); -#ifdef _OPENMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - for (MFIter mfi(vel); mfi.isValid(); ++mfi) { - Box const& gbx = amrex::grow(mfi.validbox(),nghost); - Box blo = gbx & dlo; - Box bhi = gbx & dhi; - Array4 const& v = vel[mfi].array(); - int gid = mfi.index(); - if (blo.ok()) { - prob_set_inflow_velocity(gid, olo, blo, v, lev, time); - } - if (bhi.ok()) { - prob_set_inflow_velocity(gid, ohi, bhi, v, lev, time); - } - } - } - } - // We make sure to only fill "nghost" ghost cells so we don't accidentally - // over-write good ghost cell values with unfilled ghost cell values - IntVect ng_vect(AMREX_D_DECL(nghost,nghost,nghost)); - vel.EnforcePeriodicity(0,AMREX_SPACEDIM,ng_vect,gm.periodicity()); -} - #ifdef AMREX_USE_EB void incflo::set_eb_velocity (int lev, amrex::Real /*time*/, MultiFab& eb_vel, int nghost) diff --git a/src/embedded_boundaries/eb_box.cpp b/src/embedded_boundaries/eb_box.cpp index 967f96e6..3c361e1d 100644 --- a/src/embedded_boundaries/eb_box.cpp +++ b/src/embedded_boundaries/eb_box.cpp @@ -33,6 +33,7 @@ void incflo::make_eb_box() Vector boxLo(AMREX_SPACEDIM), boxHi(AMREX_SPACEDIM); Real offset = 1.0e-15; + bool inside = true; for(int i = 0; i < AMREX_SPACEDIM; i++) { @@ -44,6 +45,7 @@ void incflo::make_eb_box() pp.queryarr("Hi", boxHi, 0, AMREX_SPACEDIM); pp.query("offset", offset); + pp.query("internal_flow", inside); Real xlo = boxLo[0] + offset; Real xhi = boxHi[0] - offset; @@ -67,27 +69,7 @@ void incflo::make_eb_box() yhi = 2.0 * geom[0].ProbHi(1) - geom[0].ProbLo(1); } -#if (AMREX_SPACEDIM == 2) - Array point_lox{xlo, 0.0}; - Array normal_lox{-1.0, 0.0}; - Array point_hix{xhi, 0.0}; - Array normal_hix{1.0, 0.0}; - - Array point_loy{0.0, ylo}; - Array normal_loy{0.0, -1.0}; - Array point_hiy{0.0, yhi}; - Array normal_hiy{0.0, 1.0}; - - EB2::PlaneIF plane_lox(point_lox, normal_lox); - EB2::PlaneIF plane_hix(point_hix, normal_hix); - - EB2::PlaneIF plane_loy(point_loy, normal_loy); - EB2::PlaneIF plane_hiy(point_hiy, normal_hiy); - - // Generate GeometryShop - auto gshop = EB2::makeShop(EB2::makeUnion(plane_lox, plane_hix, - plane_loy, plane_hiy)); -#else +#if (AMREX_SPACEDIM > 2) Real zlo = boxLo[2] + offset; Real zhi = boxHi[2] - offset; @@ -98,36 +80,15 @@ void incflo::make_eb_box() zlo = 2.0 * geom[0].ProbLo(2) - geom[0].ProbHi(2); zhi = 2.0 * geom[0].ProbHi(2) - geom[0].ProbLo(2); } +#endif - Array point_lox{xlo, 0.0, 0.0}; - Array normal_lox{-1.0, 0.0, 0.0}; - Array point_hix{xhi, 0.0, 0.0}; - Array normal_hix{1.0, 0.0, 0.0}; - - Array point_loy{0.0, ylo, 0.0}; - Array normal_loy{0.0, -1.0, 0.0}; - Array point_hiy{0.0, yhi, 0.0}; - Array normal_hiy{0.0, 1.0, 0.0}; - - Array point_loz{0.0, 0.0, zlo}; - Array normal_loz{0.0, 0.0, -1.0}; - Array point_hiz{0.0, 0.0, zhi}; - Array normal_hiz{0.0, 0.0, 1.0}; - - EB2::PlaneIF plane_lox(point_lox, normal_lox); - EB2::PlaneIF plane_hix(point_hix, normal_hix); - - EB2::PlaneIF plane_loy(point_loy, normal_loy); - EB2::PlaneIF plane_hiy(point_hiy, normal_hiy); + RealArray lo {AMREX_D_DECL(xlo, ylo, zlo)}; + RealArray hi {AMREX_D_DECL(xhi, yhi, zhi)}; - EB2::PlaneIF plane_loz(point_loz, normal_loz); - EB2::PlaneIF plane_hiz(point_hiz, normal_hiz); + EB2::BoxIF my_box(lo, hi, inside); // Generate GeometryShop - auto gshop = EB2::makeShop(EB2::makeUnion(plane_lox, plane_hix, - plane_loy, plane_hiy, - plane_loz, plane_hiz)); -#endif + auto gshop = EB2::makeShop(my_box); // Build index space int max_level_here = 0; diff --git a/src/incflo.H b/src/incflo.H index b8a81475..b8368f30 100644 --- a/src/incflo.H +++ b/src/incflo.H @@ -122,7 +122,6 @@ public: // /////////////////////////////////////////////////////////////////////////// - void set_inflow_velocity (int lev, amrex::Real time, amrex::MultiFab& vel, int nghost); #ifdef AMREX_USE_EB void set_eb_velocity (int lev, amrex::Real time, amrex::MultiFab& eb_vel, int nghost); void set_eb_density (int lev, amrex::Real time, amrex::MultiFab& eb_density, int nghost); @@ -226,8 +225,6 @@ public: /////////////////////////////////////////////////////////////////////////// void prob_init_fluid (int lev); - void prob_set_inflow_velocity (int grid_id, amrex::Orientation ori, amrex::Box const& bx, - amrex::Array4 const& v, int lev, amrex::Real time); #include "incflo_prob_I.H" #include "incflo_prob_usr_I.H" diff --git a/src/prob/CMakeLists.txt b/src/prob/CMakeLists.txt index 951190f8..0b8791aa 100644 --- a/src/prob/CMakeLists.txt +++ b/src/prob/CMakeLists.txt @@ -4,7 +4,6 @@ target_sources(incflo PRIVATE incflo_prob_I.H incflo_prob_usr_I.H - prob_bc.cpp prob_bc.H prob_init_fluid.cpp prob_init_fluid_usr.cpp diff --git a/src/prob/Make.package b/src/prob/Make.package index fd22168a..a922cb47 100644 --- a/src/prob/Make.package +++ b/src/prob/Make.package @@ -1,4 +1,3 @@ -CEXE_sources += prob_bc.cpp CEXE_sources += prob_init_fluid.cpp CEXE_sources += prob_init_fluid_usr.cpp diff --git a/src/prob/prob_bc.cpp b/src/prob/prob_bc.cpp deleted file mode 100644 index 42ef1936..00000000 --- a/src/prob/prob_bc.cpp +++ /dev/null @@ -1,101 +0,0 @@ -#include - -using namespace amrex; - -void incflo::prob_set_inflow_velocity (int /*grid_id*/, Orientation ori, Box const& bx, - Array4 const& vel, int lev, Real /*time*/) -{ - if (6 == m_probtype) - { - AMREX_D_TERM(Real u = m_ic_u;, - Real v = m_ic_v;, - Real w = m_ic_w;); - - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - AMREX_D_TERM(vel(i,j,k,0) = u;, - vel(i,j,k,1) = v;, - vel(i,j,k,2) = w;); - }); - } - else if (31 == m_probtype) - { - Real dyinv = Real(1.0) / Geom(lev).Domain().length(1); - Real u = m_ic_u; - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real y = Real(j+0.5)*dyinv; - vel(i,j,k,0) = Real(6.0) * u * y * (Real(1.0)-y); - }); - } -#if (AMREX_SPACEDIM == 3) - else if (311 == m_probtype) - { - Real dzinv = Real(1.0) / Geom(lev).Domain().length(2); - Real u = m_ic_u; - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real z = Real(k+0.5)*dzinv; - vel(i,j,k,0) = Real(6.0) * u * z * (Real(1.0)-z); - }); - } - else if (41 == m_probtype) - { - Real dzinv = Real(1.0) / Geom(lev).Domain().length(2); - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real z = Real(k+0.5)*dzinv; - vel(i,j,k,0) = Real(0.5)*z; - }); - } - else if (32 == m_probtype) - { - Real dzinv = Real(1.0) / Geom(lev).Domain().length(2); - Real v = m_ic_v; - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real z = Real(k+0.5)*dzinv; - vel(i,j,k,1) = Real(6.0) * v * z * (Real(1.0)-z); - }); - } -#endif - else if (322 == m_probtype) - { - Real dxinv = Real(1.0) / Geom(lev).Domain().length(0); - Real v = m_ic_v; - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real x = Real(i+0.5)*dxinv; - vel(i,j,k,1) = Real(6.0) * v * x * (Real(1.0)-x); - }); - } - else if (33 == m_probtype) - { - Real dxinv = Real(1.0) / Geom(lev).Domain().length(0); - Real w = m_ic_w; - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real x = Real(i+0.5)*dxinv; - vel(i,j,k,2) = Real(6.0) * w * x * (Real(1.0)-x); - }); - } - else if (333 == m_probtype) - { - Real dyinv = Real(1.0) / Geom(lev).Domain().length(1); - Real w = m_ic_w; - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real y = Real(j+0.5)*dyinv; - vel(i,j,k,2) = Real(6.0) * w * y * (Real(1.0)-y); - }); - } - else - { - const int dir = ori.coordDir(); - const Real bcv = m_bc_velocity[ori][dir]; - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - vel(i,j,k,dir) = bcv; - }); - }; -} diff --git a/src/projection/incflo_apply_nodal_projection.cpp b/src/projection/incflo_apply_nodal_projection.cpp index 81d0f5c4..71044a18 100644 --- a/src/projection/incflo_apply_nodal_projection.cpp +++ b/src/projection/incflo_apply_nodal_projection.cpp @@ -1,5 +1,7 @@ #include +#include #include +#include #include using namespace amrex; @@ -107,7 +109,27 @@ void incflo::ApplyNodalProjection (Vector density, vel.push_back(&(m_leveldata[lev]->velocity)); vel[lev]->setBndry(0.0); if (!proj_for_small_dt && !incremental) { - set_inflow_velocity(lev, time, *vel[lev], 1); + // Only the inflow boundary gets set here + IntVect nghost(1); + amrex::Vector inflow_bcr; + inflow_bcr.resize(AMREX_SPACEDIM); + for (OrientationIter oit; oit; ++oit) { + if (m_bc_type[oit()] == BC::mass_inflow) { + AMREX_D_TERM(inflow_bcr[0].set(oit(), BCType::ext_dir);, + inflow_bcr[1].set(oit(), BCType::ext_dir);, + inflow_bcr[2].set(oit(), BCType::ext_dir)); + } + } + + PhysBCFunct > physbc + (geom[lev], inflow_bcr, IncfloVelFill{m_probtype, m_bc_velocity}); + physbc(*vel[lev], 0, AMREX_SPACEDIM, nghost, time, 0); + + // Note that we wouldn't need this if we trusted users to supply inflow data + // that obeys periodicity + // We make sure to only fill "nghost" ghost cells so we don't accidentally + // over-write good ghost cell values with unfilled ghost cell values + vel[lev]->EnforcePeriodicity(0, AMREX_SPACEDIM, nghost, geom[lev].periodicity()); } } diff --git a/src/projection/incflo_projection_bc.cpp b/src/projection/incflo_projection_bc.cpp index 17b695d5..de4bd2de 100644 --- a/src/projection/incflo_projection_bc.cpp +++ b/src/projection/incflo_projection_bc.cpp @@ -20,9 +20,9 @@ incflo::get_nodal_projection_bc (Orientation::Side side) const noexcept break; } case BC::mass_inflow: - { + { r[dir] = LinOpBCType::inflow; - break; + break; } case BC::slip_wall: case BC::no_slip_wall: