put in basic data structures / operations for cc-projection as altern…

src/boundary_conditions/CMakeLists.txt

@@ -4,4 +4,5 @@ target_sources(incflo
    incflo_fillpatch.cpp
    incflo_fillphysbc.cpp
    incflo_set_bcs.cpp
+   incflo_set_velocity_bcs.cpp
    )

src/boundary_conditions/Make.package

@@ -2,3 +2,4 @@
 CEXE_sources += boundary_conditions.cpp
 CEXE_sources += incflo_fillpatch.cpp incflo_fillphysbc.cpp
 CEXE_sources += incflo_set_bcs.cpp
+CEXE_sources += incflo_set_velocity_bcs.cpp

src/boundary_conditions/incflo_set_velocity_bcs.cpp

@@ -0,0 +1,39 @@
+#ifdef AMREX_USE_EB
+#include <AMReX_EBMultiFabUtil.H>
+#endif
+
+#include <incflo.H>
+
+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);
+                }
+            }
+        }
+    }
+    vel.EnforcePeriodicity(gm.periodicity());
+}

src/incflo.H

@@ -123,6 +123,8 @@ public:
     amrex::iMultiFab make_nodalBC_mask (int lev);
     amrex::Vector<amrex::MultiFab> make_robinBC_MFs(int lev, amrex::MultiFab* state = nullptr);
 
+    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);
@@ -249,6 +251,9 @@ public:
                                       amrex::Array4<amrex::Real const> const& bcval,
                                       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"
 
@@ -258,7 +263,10 @@ public:
     //
     ///////////////////////////////////////////////////////////////////////////
 
-    void ApplyProjection(amrex::Vector<amrex::MultiFab const*> density,
+    void ApplyProjection(             amrex::Vector<amrex::MultiFab const*> const& density,
+                         AMREX_D_DECL(amrex::Vector<amrex::MultiFab*> const& u_mac,
+                                      amrex::Vector<amrex::MultiFab*> const& v_mac,
+                                      amrex::Vector<amrex::MultiFab*> const& w_mac),
                          amrex::Real time, amrex::Real scaling_factor, bool incremental);
     void ApplyNodalProjection(amrex::Vector<amrex::MultiFab const*> density,
                               amrex::Real time, amrex::Real scaling_factor, bool incremental);
@@ -267,6 +275,11 @@ public:
                               amrex::Vector<amrex::MultiFab      *> const& divu_Source,
                               amrex::Real time, amrex::Real scaling_factor, bool incremental,
                               bool set_inflow_bc);
+    void ApplyCCProjection(amrex::Vector<amrex::MultiFab const*> density,
+                           AMREX_D_DECL(amrex::Vector<amrex::MultiFab*> const& u_mac,
+                                        amrex::Vector<amrex::MultiFab*> const& v_mac,
+                                        amrex::Vector<amrex::MultiFab*> const& w_mac),
+                           amrex::Real time, amrex::Real scaling_factor, bool incremental);
 
     ///////////////////////////////////////////////////////////////////////////
     //
@@ -508,6 +521,9 @@ private:
     //    use gradient of mac phi which contains the full pressure
     bool m_use_mac_phi_in_godunov = false;
 
+    // If true use CC projection; if false use nodal projection
+    bool m_use_cc_proj = false;
+
     enum struct DiffusionType {
         Invalid, Explicit, Crank_Nicolson, Implicit
     };
@@ -554,7 +570,7 @@ private:
     int m_plt_gpz         = 1;
     int m_plt_rho         = 1;
     int m_plt_tracer      = 1;
-    int m_plt_p_nd        = 0;
+    int m_plt_p           = 0;
     int m_plt_macphi      = 0;
     int m_plt_eta         = 0;
     int m_plt_magvel      = 1;
@@ -615,7 +631,8 @@ private:
 
         amrex::MultiFab mac_phi; // cell-centered pressure used in MAC projection
 
-        // nodal pressure multifab
+        // Pressure MultiFabs (only one will actually be used)
+        amrex::MultiFab p_cc;
         amrex::MultiFab p_nd;
 
         // cell-centered pressure gradient
@@ -819,7 +836,7 @@ private:
         return m_bcrec_force_d.data(); }
 
     [[nodiscard]] amrex::Array<amrex::LinOpBCType,AMREX_SPACEDIM>
-    get_nodal_projection_bc (amrex::Orientation::Side side) const noexcept;
+    get_projection_bc (amrex::Orientation::Side side) const noexcept;
     [[nodiscard]] amrex::Array<amrex::LinOpBCType,AMREX_SPACEDIM>
     get_mac_projection_bc (amrex::Orientation::Side side) const noexcept;
 

src/incflo.cpp

@@ -196,11 +196,20 @@ void incflo::Evolve()
 }
 
 void
-incflo::ApplyProjection (Vector<MultiFab const*> density,
+incflo::ApplyProjection (Vector<MultiFab const*> const& density,
+                         AMREX_D_DECL(Vector<MultiFab*> const& u_mac,
+                                      Vector<MultiFab*> const& v_mac,
+                                      Vector<MultiFab*> const& w_mac),
                          Real time, Real scaling_factor, bool incremental)
 {
     BL_PROFILE("incflo::ApplyProjection");
-    ApplyNodalProjection(std::move(density),time,scaling_factor,incremental);
+    if (m_use_cc_proj)
+    {
+        ApplyCCProjection(density,AMREX_D_DECL(u_mac,v_mac,w_mac),
+                          time,scaling_factor,incremental);
+    }
+    else
+        ApplyNodalProjection(density,time,scaling_factor,incremental);
 }
 
 // Make a new level from scratch using provided BoxArray and DistributionMapping.

src/incflo_apply_corrector.cpp

@@ -161,7 +161,9 @@ void incflo::ApplyCorrector()
     // Project velocity field, update pressure
     // **********************************************************************************************
     bool incremental_projection = false;
-    ApplyProjection(get_density_nph_const(), new_time,m_dt,incremental_projection);
+    ApplyProjection(get_density_nph_const(),
+                    AMREX_D_DECL(GetVecOfPtrs(u_mac), GetVecOfPtrs(v_mac),
+                    GetVecOfPtrs(w_mac)),new_time,m_dt,incremental_projection);
 
 #ifdef AMREX_USE_EB
     // **********************************************************************************************

src/incflo_apply_predictor.cpp

@@ -194,7 +194,9 @@ void incflo::ApplyPredictor (bool incremental_projection)
     // **********************************************************************************************
     // Project velocity field, update pressure
     // **********************************************************************************************
-    ApplyProjection(get_density_nph_const(),new_time,m_dt,incremental_projection);
+    ApplyProjection(get_density_nph_const(),
+                    AMREX_D_DECL(GetVecOfPtrs(u_mac), GetVecOfPtrs(v_mac),
+                    GetVecOfPtrs(w_mac)),new_time,m_dt,incremental_projection);
 
 #ifdef INCFLO_USE_PARTICLES
     // **************************************************************************************

src/incflo_regrid.cpp

@@ -34,6 +34,8 @@ void incflo::MakeNewLevelFromCoarse (int lev,
         fillcoarsepatch_tracer(lev, time, new_leveldata->tracer, 0);
     }
     fillcoarsepatch_gradp(lev, time, new_leveldata->gp, 0);
+
+    new_leveldata->p_cc.setVal(0.0);
     new_leveldata->p_nd.setVal(0.0);
 
     m_leveldata[lev] = std::move(new_leveldata);
@@ -82,6 +84,8 @@ void incflo::RemakeLevel (int lev, Real time, const BoxArray& ba,
         fillpatch_tracer(lev, time, new_leveldata->tracer, 0);
     }
     fillpatch_gradp(lev, time, new_leveldata->gp, 0);
+
+    new_leveldata->p_cc.setVal(0.0);
     new_leveldata->p_nd.setVal(0.0);
 
     m_leveldata[lev] = std::move(new_leveldata);

src/prob/prob_bc.cpp

@@ -45,13 +45,15 @@ void incflo::prob_set_BC_MF (Orientation const& ori, Box const& bx,
                             mask(i,j,k,n) = inflow_val;
                         }
                     }
+#if (AMREX_SPACEDIM == 3)
                     else if (direction == 2) {
                         if (k <= half_num_cells) {
                             mask(i,j,k,n) = outflow_val;
                         } else {
                             mask(i,j,k,n) = inflow_val;
                         }
                     }
+#endif
                 }
             });
         } else {
@@ -72,13 +74,15 @@ void incflo::prob_set_BC_MF (Orientation const& ori, Box const& bx,
                             mask(i,j,k,n) = inflow_val;
                         }
                     }
+#if (AMREX_SPACEDIM == 3)
                     else if (direction == 2) {
                         if (k > half_num_cells) {
                             mask(i,j,k,n) = outflow_val;
                         } else {
                             mask(i,j,k,n) = inflow_val;
                         }
                     }
+#endif
                 }
             });
         }
@@ -136,6 +140,7 @@ void incflo::prob_set_MAC_robinBCs (Orientation const& ori, Box const& bx,
                         robin_b(i,j,k) = 1.;
                     }
                 }
+#if (AMREX_SPACEDIM == 3)
                 else if (direction == 2) {
                     if (k <= half_num_cells) {
                         robin_a(i,j,k) = 1.;
@@ -145,6 +150,7 @@ void incflo::prob_set_MAC_robinBCs (Orientation const& ori, Box const& bx,
                         robin_b(i,j,k) = 1.;
                     }
                 }
+#endif
             });
         } else {
             ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
@@ -169,6 +175,7 @@ void incflo::prob_set_MAC_robinBCs (Orientation const& ori, Box const& bx,
                         robin_b(i,j,k) = 1.;
                     }
                 }
+#if (AMREX_SPACEDIM == 3)
                 else if (direction == 2) {
                     if (k > half_num_cells) {
                         robin_a(i,j,k) = 1.;
@@ -178,6 +185,7 @@ void incflo::prob_set_MAC_robinBCs (Orientation const& ori, Box const& bx,
                         robin_b(i,j,k) = 1.;
                     }
                 }
+#endif
             });
         }
     }
@@ -237,6 +245,7 @@ void incflo::prob_set_diffusion_robinBCs (Orientation const& ori, Box const& bx,
                         robin_f(i,j,k,0) = bcval(i,j,k,0);
                     }
                 }
+#if (AMREX_SPACEDIM == 3)
                 else if (direction == 2) {
                     if (k <= half_num_cells) {
                         robin_a(i,j,k,0) = 0.;
@@ -248,6 +257,7 @@ void incflo::prob_set_diffusion_robinBCs (Orientation const& ori, Box const& bx,
                         robin_f(i,j,k,0) = bcval(i,j,k,0);
                     }
                 }
+#endif
             });
         } else {
             ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
@@ -276,6 +286,7 @@ void incflo::prob_set_diffusion_robinBCs (Orientation const& ori, Box const& bx,
                         robin_f(i,j,k,0) = bcval(i,j,k,0);
                     }
                 }
+#if (AMREX_SPACEDIM == 3)
                 else if (direction == 2) {
                     if (k > half_num_cells) {
                         robin_a(i,j,k,0) = 0.;
@@ -287,6 +298,7 @@ void incflo::prob_set_diffusion_robinBCs (Orientation const& ori, Box const& bx,
                         robin_f(i,j,k,0) = bcval(i,j,k,0);
                     }
                 }
+#endif
             });
         }
     }
@@ -296,3 +308,101 @@ void incflo::prob_set_diffusion_robinBCs (Orientation const& ori, Box const& bx,
               +std::to_string(m_probtype));
     }
 }
+
+void incflo::prob_set_inflow_velocity (int /*grid_id*/, Orientation ori, Box const& bx,
+                                       Array4<Real> 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 = 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 = (j+0.5)*dyinv;
+            vel(i,j,k,0) = 6. * u * y * (1.-y);
+        });
+    }
+#if (AMREX_SPACEDIM == 3)
+    else if (311 == m_probtype)
+    {
+        Real dzinv = 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 = (k+0.5)*dzinv;
+            vel(i,j,k,0) = 6. * u * z * (1.-z);
+        });
+    }
+    else if (41 == m_probtype)
+    {
+        Real dzinv = 1.0 / Geom(lev).Domain().length(2);
+        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+        {
+            Real z = (k+0.5)*dzinv;
+            vel(i,j,k,0) = 0.5*z;
+        });
+    }
+    else if (32 == m_probtype)
+    {
+        Real dzinv = 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 = (k+0.5)*dzinv;
+            vel(i,j,k,1) = 6. * v * z * (1.-z);
+        });
+    }
+#endif
+    else if (322 == m_probtype)
+    {
+        Real dxinv = 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 = (i+0.5)*dxinv;
+            vel(i,j,k,1) = 6. * v * x * (1.-x);
+        });
+    }
+    else if (33 == m_probtype)
+    {
+        Real dxinv = 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 = (i+0.5)*dxinv;
+            vel(i,j,k,2) = 6. * w * x * (1.-x);
+        });
+    }
+    else if (333 == m_probtype)
+    {
+        Real dyinv = 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 = (j+0.5)*dyinv;
+            vel(i,j,k,2) = 6. * w * y * (1.-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;
+        });
+    };
+}