From c9b947677792c16ff4e3b2b51d1961a6998816d8 Mon Sep 17 00:00:00 2001 From: cgilet Date: Fri, 19 Jan 2024 11:35:12 -0500 Subject: [PATCH 1/9] Rename file to better reflect contents. renamed: incflo_set_velocity_bcs.cpp -> incflo_set_bcs.cpp --- .../{incflo_set_velocity_bcs.cpp => incflo_set_bcs.cpp} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename src/boundary_conditions/{incflo_set_velocity_bcs.cpp => incflo_set_bcs.cpp} (100%) diff --git a/src/boundary_conditions/incflo_set_velocity_bcs.cpp b/src/boundary_conditions/incflo_set_bcs.cpp similarity index 100% rename from src/boundary_conditions/incflo_set_velocity_bcs.cpp rename to src/boundary_conditions/incflo_set_bcs.cpp From 4c2082b2deb6dbe399931bdc8b52a649899303e2 Mon Sep 17 00:00:00 2001 From: cgilet Date: Mon, 15 Jan 2024 14:10:28 -0500 Subject: [PATCH 2/9] Update make_eb_box() to use amrex::BoxIF. This allows specifying internal or external flow. --- src/embedded_boundaries/eb_box.cpp | 61 ++++++------------------------ 1 file changed, 11 insertions(+), 50 deletions(-) diff --git a/src/embedded_boundaries/eb_box.cpp b/src/embedded_boundaries/eb_box.cpp index 967f96e6..ec04acc2 100644 --- a/src/embedded_boundaries/eb_box.cpp +++ b/src/embedded_boundaries/eb_box.cpp @@ -33,7 +33,8 @@ 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++) { boxLo[i] = geom[0].ProbLo(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,37 +80,16 @@ 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); - - EB2::PlaneIF plane_loz(point_loz, normal_loz); - EB2::PlaneIF plane_hiz(point_hiz, normal_hiz); + RealArray lo {AMREX_D_DECL(xlo, ylo, zlo)}; + RealArray hi {AMREX_D_DECL(xhi, yhi, zhi)}; + + 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; int max_coarsening_level = 100; From 99f2d28cc97e1281aaceb46d244fec4da722a4e6 Mon Sep 17 00:00:00 2001 From: cgilet Date: Fri, 19 Jan 2024 17:05:30 -0500 Subject: [PATCH 3/9] Remove incflo set_velocity_bc functions in favor of using PhysBCFunct, which we use to set the inflow bcs during FillPatch. --- src/boundary_conditions/incflo_set_bcs.cpp | 35 ------ src/incflo.H | 3 - src/prob/prob_bc.cpp | 101 ------------------ .../incflo_apply_nodal_projection.cpp | 11 +- 4 files changed, 10 insertions(+), 140 deletions(-) delete mode 100644 src/prob/prob_bc.cpp diff --git a/src/boundary_conditions/incflo_set_bcs.cpp b/src/boundary_conditions/incflo_set_bcs.cpp index d01a854a..a2ef9713 100644 --- a/src/boundary_conditions/incflo_set_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/incflo.H b/src/incflo.H index 63af0acd..93264290 100644 --- a/src/incflo.H +++ b/src/incflo.H @@ -110,7 +110,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); @@ -211,8 +210,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/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..166b8956 100644 --- a/src/projection/incflo_apply_nodal_projection.cpp +++ b/src/projection/incflo_apply_nodal_projection.cpp @@ -107,7 +107,16 @@ 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); + //set_inflow_velocity(lev, time, *vel[lev], 1); + PhysBCFunct > physbc + (geom[lev], bcrec, IncfloVelFill{m_probtype, m_bc_velocity}); + physbcf(*vel[lev], 0, AMREX_SPACEDIM, 1, time, 0); + + //FIXME? Not sure we really need this + // 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.EnforcePeriodicity(0,AMREX_SPACEDIM,IntVect::Unit,gm.periodicity()); + } } From 08a16bb855d5d04e3cfe7cfcacec1165cc806155 Mon Sep 17 00:00:00 2001 From: cgilet Date: Fri, 19 Jan 2024 17:06:38 -0500 Subject: [PATCH 4/9] Fix spacing --- src/projection/incflo_projection_bc.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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: From 54c9b94e4490036950e020cd6788dd9d9d594eb3 Mon Sep 17 00:00:00 2001 From: cgilet Date: Fri, 19 Jan 2024 17:08:42 -0500 Subject: [PATCH 5/9] Update Make.package --- src/boundary_conditions/Make.package | 1 - src/prob/Make.package | 1 - 2 files changed, 2 deletions(-) diff --git a/src/boundary_conditions/Make.package b/src/boundary_conditions/Make.package index f7d0b486..a2751032 100644 --- a/src/boundary_conditions/Make.package +++ b/src/boundary_conditions/Make.package @@ -1,4 +1,3 @@ CEXE_sources += boundary_conditions.cpp CEXE_sources += incflo_fillpatch.cpp incflo_fillphysbc.cpp -CEXE_sources += incflo_set_velocity_bcs.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 From f4c34a1d2b62fca3e1f2d546e5ad850035f80d0b Mon Sep 17 00:00:00 2001 From: cgilet Date: Fri, 19 Jan 2024 17:22:14 -0500 Subject: [PATCH 6/9] Fixes to get compiling --- src/boundary_conditions/Make.package | 1 + src/projection/incflo_apply_nodal_projection.cpp | 10 ++++++---- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/src/boundary_conditions/Make.package b/src/boundary_conditions/Make.package index a2751032..9de2ae89 100644 --- a/src/boundary_conditions/Make.package +++ b/src/boundary_conditions/Make.package @@ -1,3 +1,4 @@ CEXE_sources += boundary_conditions.cpp CEXE_sources += incflo_fillpatch.cpp incflo_fillphysbc.cpp +CEXE_sources += incflo_set_bcs.cpp diff --git a/src/projection/incflo_apply_nodal_projection.cpp b/src/projection/incflo_apply_nodal_projection.cpp index 166b8956..685c5f16 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,15 +109,15 @@ 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); + IntVect nghost(1); PhysBCFunct > physbc - (geom[lev], bcrec, IncfloVelFill{m_probtype, m_bc_velocity}); - physbcf(*vel[lev], 0, AMREX_SPACEDIM, 1, time, 0); + (geom[lev], get_velocity_bcrec(), IncfloVelFill{m_probtype, m_bc_velocity}); + physbc(*vel[lev], 0, AMREX_SPACEDIM, nghost, time, 0); //FIXME? Not sure we really need this // 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.EnforcePeriodicity(0,AMREX_SPACEDIM,IntVect::Unit,gm.periodicity()); + vel[lev]->EnforcePeriodicity(0, AMREX_SPACEDIM, nghost, geom[lev].periodicity()); } } From 09d0ea182a00374d42eadfe1e6adbc435eb02aa3 Mon Sep 17 00:00:00 2001 From: cgilet Date: Fri, 19 Jan 2024 17:41:35 -0500 Subject: [PATCH 7/9] Fix trailing whitespace --- src/embedded_boundaries/eb_box.cpp | 6 +++--- src/projection/incflo_apply_nodal_projection.cpp | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/embedded_boundaries/eb_box.cpp b/src/embedded_boundaries/eb_box.cpp index ec04acc2..3c361e1d 100644 --- a/src/embedded_boundaries/eb_box.cpp +++ b/src/embedded_boundaries/eb_box.cpp @@ -34,7 +34,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++) { boxLo[i] = geom[0].ProbLo(i); @@ -84,12 +84,12 @@ void incflo::make_eb_box() RealArray lo {AMREX_D_DECL(xlo, ylo, zlo)}; RealArray hi {AMREX_D_DECL(xhi, yhi, zhi)}; - + EB2::BoxIF my_box(lo, hi, inside); // Generate GeometryShop auto gshop = EB2::makeShop(my_box); - + // Build index space int max_level_here = 0; int max_coarsening_level = 100; diff --git a/src/projection/incflo_apply_nodal_projection.cpp b/src/projection/incflo_apply_nodal_projection.cpp index 685c5f16..ee7e2d93 100644 --- a/src/projection/incflo_apply_nodal_projection.cpp +++ b/src/projection/incflo_apply_nodal_projection.cpp @@ -109,7 +109,7 @@ void incflo::ApplyNodalProjection (Vector density, vel.push_back(&(m_leveldata[lev]->velocity)); vel[lev]->setBndry(0.0); if (!proj_for_small_dt && !incremental) { - IntVect nghost(1); + IntVect nghost(1); PhysBCFunct > physbc (geom[lev], get_velocity_bcrec(), IncfloVelFill{m_probtype, m_bc_velocity}); physbc(*vel[lev], 0, AMREX_SPACEDIM, nghost, time, 0); From a3d786d22d1771e1f45431e9c1b0b3ce2cdc2fd9 Mon Sep 17 00:00:00 2001 From: cgilet Date: Fri, 19 Jan 2024 17:45:54 -0500 Subject: [PATCH 8/9] Fix CMakeLists.txt --- src/boundary_conditions/CMakeLists.txt | 2 +- src/prob/CMakeLists.txt | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) 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/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 From 04a31a629ea000ac6900ba0f8d1bc890a66cdc0c Mon Sep 17 00:00:00 2001 From: cgilet Date: Fri, 19 Jan 2024 21:26:24 -0500 Subject: [PATCH 9/9] Make it so only the inflow BCs get set in apply_nodal_projection The rest should stay zero. --- .../incflo_apply_nodal_projection.cpp | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/src/projection/incflo_apply_nodal_projection.cpp b/src/projection/incflo_apply_nodal_projection.cpp index ee7e2d93..71044a18 100644 --- a/src/projection/incflo_apply_nodal_projection.cpp +++ b/src/projection/incflo_apply_nodal_projection.cpp @@ -109,16 +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) { + // 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], get_velocity_bcrec(), IncfloVelFill{m_probtype, m_bc_velocity}); + (geom[lev], inflow_bcr, IncfloVelFill{m_probtype, m_bc_velocity}); physbc(*vel[lev], 0, AMREX_SPACEDIM, nghost, time, 0); - //FIXME? Not sure we really need this + // 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()); - } }