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

…ative to nodal projection

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"
 
@@ -259,6 +264,9 @@ public:
     ///////////////////////////////////////////////////////////////////////////
 
     void ApplyProjection(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);
     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

@@ -197,10 +197,19 @@ void incflo::Evolve()
 
 void
 incflo::ApplyProjection (Vector<MultiFab 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

@@ -296,3 +296,99 @@ 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);
+        });
+    }
+    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);
+        });
+    }
+    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;
+        });
+    };
+}

src/prob/prob_init_fluid.cpp

@@ -10,7 +10,9 @@ void incflo::prob_init_fluid (int lev)
     auto const& problo = geom[lev].ProbLoArray();
     auto const& probhi = geom[lev].ProbHiArray();
 
+    ld.p_cc.setVal(0.0);
     ld.p_nd.setVal(0.0);
+
     ld.gp.setVal(0.0);
 
     ld.density.setVal(m_ro_0);

src/projection/CMakeLists.txt

@@ -3,5 +3,6 @@ target_include_directories(incflo PRIVATE ${CMAKE_CURRENT_LIST_DIR})
 target_sources(incflo
    PRIVATE
    incflo_projection_bc.cpp
+   incflo_apply_cc_projection.cpp
    incflo_apply_nodal_projection.cpp
    )

src/projection/Make.package

@@ -1,2 +1,3 @@
 CEXE_sources += incflo_projection_bc.cpp
+CEXE_sources += incflo_apply_cc_projection.cpp
 CEXE_sources += incflo_apply_nodal_projection.cpp