Consolidate code (AMReX-Fluids#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

@@ -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
    )

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

src/boundary_conditions/incflo_set_velocity_bcs.cpp → 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<Real> 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)

src/embedded_boundaries/eb_box.cpp

@@ -33,6 +33,7 @@ void incflo::make_eb_box()
 
         Vector<Real> 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<Real, 2> point_lox{xlo, 0.0};
-        Array<Real, 2> normal_lox{-1.0, 0.0};
-        Array<Real, 2> point_hix{xhi, 0.0};
-        Array<Real, 2> normal_hix{1.0, 0.0};
-
-        Array<Real, 2> point_loy{0.0, ylo};
-        Array<Real, 2> normal_loy{0.0, -1.0};
-        Array<Real, 2> point_hiy{0.0, yhi};
-        Array<Real, 2> 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<Real, 3> point_lox{xlo, 0.0, 0.0};
-        Array<Real, 3> normal_lox{-1.0, 0.0, 0.0};
-        Array<Real, 3> point_hix{xhi, 0.0, 0.0};
-        Array<Real, 3> normal_hix{1.0, 0.0, 0.0};
-
-        Array<Real, 3> point_loy{0.0, ylo, 0.0};
-        Array<Real, 3> normal_loy{0.0, -1.0, 0.0};
-        Array<Real, 3> point_hiy{0.0, yhi, 0.0};
-        Array<Real, 3> normal_hiy{0.0, 1.0, 0.0};
-
-        Array<Real, 3> point_loz{0.0, 0.0, zlo};
-        Array<Real, 3> normal_loz{0.0, 0.0, -1.0};
-        Array<Real, 3> point_hiz{0.0, 0.0, zhi};
-        Array<Real, 3> 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;

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<amrex::Real> const& v, int lev, amrex::Real time);
 
 #include "incflo_prob_I.H"
 #include "incflo_prob_usr_I.H"

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

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
 

src/projection/incflo_apply_nodal_projection.cpp

@@ -1,5 +1,7 @@
 #include <AMReX_BC_TYPES.H>
+#include <AMReX_PhysBCFunct.H>
 #include <incflo.H>
+#include <prob_bc.H>
 #include <memory>
 
 using namespace amrex;
@@ -107,7 +109,27 @@ void incflo::ApplyNodalProjection (Vector<MultiFab const*> 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<amrex::BCRec> 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<GpuBndryFuncFab<IncfloVelFill> > 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());
         }
     }
 

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: