WIP - add function to create mixedBC mask, although maybe it would

be better to use a PhysBCFunct...

src/boundary_conditions/incflo_set_bcs.cpp

@@ -41,6 +41,36 @@ incflo::set_inflow_velocity (int lev, amrex::Real time, MultiFab& vel, int nghos
     vel.EnforcePeriodicity(0,AMREX_SPACEDIM,ng_vect,gm.periodicity());
 }
 
+void
+incflo::set_mixedBC_mask (int lev, amrex::Real time, MultiFab& mask)
+{
+    Geometry const& gm = Geom(lev);
+    Box const& domain = gm.Domain();
+    for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) {
+        Orientation olo(dir,Orientation::low);
+        Orientation ohi(dir,Orientation::high);
+        if (m_bc_type[olo] == BC::mixed || m_bc_type[ohi] == BC::mixed) {
+            Box dlo = (m_bc_type[olo] == BC::mixed) ? surroundingNodes(bdryLo(domain,dir)) : Box();
+            Box dhi = (m_bc_type[ohi] == BC::mixed) ? surroundingNodes(bdryhi(domain,dir)) : Box();
+#ifdef _OPENMP
+#pragma omp parallel if (Gpu::notInLaunchRegion())
+#endif
+            for (MFIter mfi(mask); mfi.isValid(); ++mfi) {
+                Box blo = mfi.validbox() & dlo;
+                Box bhi = mfi.validbox() & dhi;
+                Array4<Real> const& mask_arr = mask[mfi].array();
+                int gid = mfi.index();
+                if (blo.ok()) {
+                    prob_set_mixedBC_mask(olo, blo, mask_arr, lev, time);
+                }
+                if (bhi.ok()) {
+                    prob_set_mixedBC_mask(ohi, bhi, mask_arr, lev, time);
+                }
+            }
+        }
+    }
+}
+
 #ifdef AMREX_USE_EB
 void
 incflo::set_eb_velocity (int lev, amrex::Real /*time*/, MultiFab& eb_vel, int nghost)

src/incflo.H

@@ -111,6 +111,7 @@ public:
     ///////////////////////////////////////////////////////////////////////////
 
     void set_inflow_velocity (int lev, amrex::Real time, amrex::MultiFab& vel, int nghost);
+    void set_mixedBC_mask (int lev, amrex::Real time, amrex::MultiFab& mask);
 #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);
@@ -213,6 +214,8 @@ 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);
+    void prob_set_mixedBC_mask (amrex::Orientation ori, amrex::Box const& bx,
+                                amrex::Array4<amrex::Real> const& mask, int lev, amrex::Real time);
 
 #include "incflo_prob_I.H"
 #include "incflo_prob_usr_I.H"
@@ -558,6 +561,11 @@ private:
         amrex::MultiFab divtau_o;
         amrex::MultiFab laps;
         amrex::MultiFab laps_o;
+
+        // mask for mixed BCs: 1 = inflow (Neumann for NodalProj, Dirichlet for advection),
+        //                     0 = outflow (Homogeneous Dirichlet for NodalProj, FOEXTRAP for advection)
+        amrex::iMultiFab mixedBC_mask;
+
     };
 
     amrex::Vector<std::unique_ptr<LevelData> > m_leveldata;
@@ -569,9 +577,6 @@ private:
         periodic, mixed, undefined
     };
 
-    // mask for mixed BCs: 1 = Neumann, 0 = Dirichlet
-    amrex::Vector<std::unique_ptr<amrex::iMultiFab>> m_BC_mask;
-
     amrex::GpuArray<BC                         , AMREX_SPACEDIM*2> m_bc_type;
     amrex::GpuArray<amrex::Real                , AMREX_SPACEDIM*2> m_bc_pressure;
     amrex::GpuArray<amrex::Real                , AMREX_SPACEDIM*2> m_bc_density;

src/prob/prob_bc.cpp

@@ -97,5 +97,61 @@ void incflo::prob_set_inflow_velocity (int /*grid_id*/, Orientation ori, Box con
         {
             vel(i,j,k,dir) = bcv;
         });
-    };
+    }
+}
+
+void incflo::prob_set_mixedBC_mask (Orientation ori, Box const& bx,
+                                    Array4<Real> const& mask, int lev, Real /*time*/)
+{
+    if (1100 == m_probtype || 1101 == m_probtype || 1102 == m_probtype)
+    {
+        const int direction = ori.coordDir();
+        int half_num_cells  = domain.length(direction) / 2;
+
+        Orientation::Side side = ori.faceDir();
+        if (side == Orientation::low) {
+            amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            {
+                if (direction == 0) {
+                    if (i <= half_num_cells) {
+                        mask(i,j,k,0) = 0; // outflow on bottom
+                    }
+                }
+                else if (direction == 1) {
+                    if (j <= half_num_cells) {
+                        mask(i,j,k,0) = 0;
+                    }
+                }
+                else if (direction == 2) {
+                    if (k <= half_num_cells) {
+                        mask(i,j,k,0) = 0;
+                    }
+                }
+            });
+        } else {
+            amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            {
+                if (direction == 0) {
+                    if (i > half_num_cells) {
+                        mask(i,j,k,0) = 0; // outflow on top
+                    }
+                }
+                else if (direction == 1) {
+                    if (j > half_num_cells) {
+                        mask(i,j,k,0) = 0;
+                    }
+                }
+                else if (direction == 2) {
+                    if (k > half_num_cells) {
+                        mask(i,j,k,0) = 0;
+                    }
+                }
+            });
+        }
+    }
+    else
+    {
+        Abort("incflo::prob_set_mixedBC_mask: No masking function for probtype "
+              +std::to_string(probtype));
+    }
 }

src/setup/incflo_arrays.cpp

@@ -46,6 +46,16 @@ incflo::LevelData::LevelData (amrex::BoxArray const& ba,
             laps_o.define(ba, dm, ntrac, 0, MFInfo(), fact);
         }
     }
+
+    // Only need BC_mask for mixed
+    for (OrientationIter oit; oit; ++oit) {
+        if (m_bc_type[oit()] == BC::mixed )
+        {
+            // MLNodeLap does not require any ghost cells...
+            mixedBC_mask.define(amrex::convert(ba,IntVect::TheNodeVector()), dm, 1, 0);
+            break;
+        }
+    }
 }
 
 // Resize all arrays when instance of incflo class is constructed.