Vof

src/Make.incflo

@@ -14,6 +14,7 @@ Bdirs += src/projection
 Bdirs += src/rheology
 Bdirs += src/setup
 Bdirs += src/utilities
+Bdirs += src/vof
 
 ifeq ($(USE_PARTICLES), TRUE)
   DEFINES += -DINCFLO_USE_PARTICLES

src/boundary_conditions/boundary_conditions.cpp

@@ -34,13 +34,22 @@ void incflo::init_bcs ()
             m_bc_type[ori] = BC::pressure_inflow;
 
             pp.get("pressure", m_bc_pressure[ori]);
-
+            pp.queryarr("tracer", m_bc_tracer[ori], 0, m_ntrac);
             // Set mathematical BCs here also
             AMREX_D_TERM(m_bcrec_velocity[0].set(ori, BCType::foextrap);,
                          m_bcrec_velocity[1].set(ori, BCType::foextrap);,
                          m_bcrec_velocity[2].set(ori, BCType::foextrap););
             m_bcrec_density[0].set(ori, BCType::foextrap);
-            for (auto& b : m_bcrec_tracer) { b.set(ori, BCType::foextrap); }
+            //when the VOF method is used, the default BC for tracer (i.e., the keyword 'tracer'
+            //is not explicitly included in the bcid) is symmetrical.
+            if ( pp.contains("tracer") ) {
+                for (auto& b : m_bcrec_tracer) { b.set(ori, BCType::ext_dir); }
+            }else if(m_vof_advect_tracer){
+                for (auto& b : m_bcrec_tracer) { b.set(ori, BCType::reflect_even); }
+            }
+            else {
+                for (auto& b : m_bcrec_tracer) { b.set(ori, BCType::foextrap); }
+            }
         }
         else if (bc_type == "pressure_outflow" || bc_type == "po")
         {
@@ -49,13 +58,31 @@ void incflo::init_bcs ()
             m_bc_type[ori] = BC::pressure_outflow;
 
             pp.get("pressure", m_bc_pressure[ori]);
-
+            pp.queryarr("tracer", m_bc_tracer[ori], 0, m_ntrac);
             // Set mathematical BCs here also
             AMREX_D_TERM(m_bcrec_velocity[0].set(ori, BCType::foextrap);,
                          m_bcrec_velocity[1].set(ori, BCType::foextrap);,
                          m_bcrec_velocity[2].set(ori, BCType::foextrap););
+            // Only normal oriection has reflect_even
+            //for (int dim = 0; dim < AMREX_SPACEDIM; dim++){
+                //if (dim !=ori.coordDir())
+                 // m_bcrec_velocity[ori.coordDir()].set(ori, BCType::reflect_even);
+                //else
+                // m_bcrec_velocity[dim].set(ori, BCType::ext_dir);
+            //}
             m_bcrec_density[0].set(ori, BCType::foextrap);
-            for (auto& b : m_bcrec_tracer) { b.set(ori, BCType::foextrap); }
+            //when the VOF method is used, the default BC for tracer (i.e., the keyword 'tracer'
+            //is not explicitly included in the bcid) is symmetrical.
+            if ( pp.contains("tracer") ) {
+                for (auto& b : m_bcrec_tracer) { b.set(ori, BCType::ext_dir); }
+            }else if(m_vof_advect_tracer){
+                for (auto& b : m_bcrec_tracer) { b.set(ori, BCType::reflect_even); }
+            }
+            else {
+                for (auto& b : m_bcrec_tracer) { b.set(ori, BCType::foextrap); }
+            }
+
+
         }
         else if (bc_type == "mass_inflow" || bc_type == "mi")
         {
@@ -132,8 +159,12 @@ void incflo::init_bcs ()
                          m_bcrec_velocity[1].set(ori, BCType::ext_dir);,
                          m_bcrec_velocity[2].set(ori, BCType::ext_dir););
             m_bcrec_density[0].set(ori, BCType::foextrap);
+            //when the VOF method is used, the default BC for tracer (i.e., the keyword 'tracer'
+            //is not explicitly included in the bcid) is symmetrical.
             if ( pp.contains("tracer") ) {
                 for (auto& b : m_bcrec_tracer) { b.set(ori, BCType::ext_dir); }
+            }else if(m_vof_advect_tracer){
+                for (auto& b : m_bcrec_tracer) { b.set(ori, BCType::reflect_even); }
             } else {
                 for (auto& b : m_bcrec_tracer) { b.set(ori, BCType::foextrap); }
             }

src/boundary_conditions/incflo_fillpatch.cpp

@@ -30,16 +30,34 @@ void incflo::fillpatch_velocity (int lev, Real time, MultiFab& vel, int ng)
 #else
         Interpolater* mapper = &cell_cons_interp;
 #endif
-        FillPatchTwoLevels(vel, IntVect(ng), time,
-                           {&(m_leveldata[lev-1]->velocity_o),
-                            &(m_leveldata[lev-1]->velocity)},
-                           {m_t_old[lev-1], m_t_new[lev-1]},
-                           {&(m_leveldata[lev]->velocity_o),
-                            &(m_leveldata[lev]->velocity)},
-                           {m_t_old[lev], m_t_new[lev]},
-                           0, 0, AMREX_SPACEDIM, geom[lev-1], geom[lev],
-                           cphysbc, 0, fphysbc, 0,
-                           refRatio(lev-1), mapper, bcrec, 0);
+        if (!m_fillpatchnlevels){
+            FillPatchTwoLevels(vel, IntVect(ng), time,
+                               {&(m_leveldata[lev-1]->velocity_o),
+                                &(m_leveldata[lev-1]->velocity)},
+                               {m_t_old[lev-1], m_t_new[lev-1]},
+                               {&(m_leveldata[lev]->velocity_o),
+                                &(m_leveldata[lev]->velocity)},
+                               {m_t_old[lev], m_t_new[lev]},
+                               0, 0, AMREX_SPACEDIM, geom[lev-1], geom[lev],
+                               cphysbc, 0, fphysbc, 0,
+                               refRatio(lev-1), mapper, bcrec, 0);
+        }else{
+            //for quad-/Oct-tree like grids, it is safter to use FillPatchNLevels
+            Vector<PhysBCFunct<GpuBndryFuncFab<IncfloVelFill>>> physbcs;
+            for (int ilev = 0; ilev <= finest_level; ++ilev) {
+                physbcs.emplace_back(geom[ilev],bcrec,IncfloVelFill{m_probtype, m_bc_velocity});
+            }
+            Vector<Vector<MultiFab*>> smf(finest_level+1);
+            Vector<Vector<Real>> st(finest_level+1);
+            for (int ilev = 0; ilev <= finest_level; ++ilev) {
+              smf[ilev] = {&(m_leveldata[ilev]->velocity_o), &(m_leveldata[ilev]->velocity)};
+              st[ilev] = {m_t_old[ilev], m_t_new[ilev]};
+            }
+            FillPatchNLevels(vel, lev, IntVect(ng), time, smf, st, 0, 0, AMREX_SPACEDIM, geom,
+                             physbcs, 0, ref_ratio, mapper, bcrec, 0);
+        }
+
+
     }
 }
 
@@ -65,16 +83,33 @@ void incflo::fillpatch_density (int lev, Real time, MultiFab& density, int ng)
 #else
         Interpolater* mapper = &cell_cons_interp;
 #endif
-        FillPatchTwoLevels(density, IntVect(ng), time,
-                           {&(m_leveldata[lev-1]->density_o),
-                            &(m_leveldata[lev-1]->density)},
-                           {m_t_old[lev-1], m_t_new[lev-1]},
-                           {&(m_leveldata[lev]->density_o),
-                            &(m_leveldata[lev]->density)},
-                           {m_t_old[lev], m_t_new[lev]},
-                           0, 0, 1, geom[lev-1], geom[lev],
-                           cphysbc, 0, fphysbc, 0,
-                           refRatio(lev-1), mapper, bcrec, 0);
+
+        if (!m_fillpatchnlevels){
+          FillPatchTwoLevels(density, IntVect(ng), time,
+                             {&(m_leveldata[lev-1]->density_o),
+                              &(m_leveldata[lev-1]->density)},
+                             {m_t_old[lev-1], m_t_new[lev-1]},
+                             {&(m_leveldata[lev]->density_o),
+                              &(m_leveldata[lev]->density)},
+                             {m_t_old[lev], m_t_new[lev]},
+                             0, 0, 1, geom[lev-1], geom[lev],
+                             cphysbc, 0, fphysbc, 0,
+                             refRatio(lev-1), mapper, bcrec, 0);
+        }else{
+            //for quad-/Oct-tree like grids, it is safter to use FillPatchNLevels
+            Vector<PhysBCFunct<GpuBndryFuncFab<IncfloDenFill>>> physbcs;
+            for (int ilev = 0; ilev <= finest_level; ++ilev) {
+                physbcs.emplace_back(geom[ilev],bcrec,IncfloDenFill{m_probtype, m_bc_density, m_bc_velocity});
+            }
+            Vector<Vector<MultiFab*>> smf(finest_level+1);
+            Vector<Vector<Real>> st(finest_level+1);
+            for (int ilev = 0; ilev <= finest_level; ++ilev) {
+              smf[ilev] = {&(m_leveldata[ilev]->density_o), &(m_leveldata[ilev]->density)};
+              st[ilev] = {m_t_old[ilev], m_t_new[ilev]};
+            }
+            FillPatchNLevels(density, lev, IntVect(ng), time, smf, st, 0, 0, 1, geom,
+                             physbcs, 0, ref_ratio, mapper, bcrec, 0);
+        }
     }
 }
 
@@ -101,16 +136,34 @@ void incflo::fillpatch_tracer (int lev, Real time, MultiFab& tracer, int ng)
 #else
         Interpolater* mapper = &cell_cons_interp;
 #endif
-        FillPatchTwoLevels(tracer, IntVect(ng), time,
-                           {&(m_leveldata[lev-1]->tracer_o),
-                            &(m_leveldata[lev-1]->tracer)},
-                           {m_t_old[lev-1], m_t_new[lev-1]},
-                           {&(m_leveldata[lev]->tracer_o),
-                            &(m_leveldata[lev]->tracer)},
-                           {m_t_old[lev], m_t_new[lev]},
-                           0, 0, m_ntrac, geom[lev-1], geom[lev],
-                           cphysbc, 0, fphysbc, 0,
-                           refRatio(lev-1), mapper, bcrec, 0);
+        if (!m_fillpatchnlevels){
+          FillPatchTwoLevels(tracer, IntVect(ng), time,
+                             {&(m_leveldata[lev-1]->tracer_o),
+                              &(m_leveldata[lev-1]->tracer)},
+                             {m_t_old[lev-1], m_t_new[lev-1]},
+                             {&(m_leveldata[lev]->tracer_o),
+                              &(m_leveldata[lev]->tracer)},
+                             {m_t_old[lev], m_t_new[lev]},
+                             0, 0, m_ntrac, geom[lev-1], geom[lev],
+                             cphysbc, 0, fphysbc, 0,
+                             refRatio(lev-1), mapper, bcrec, 0);
+        }else{
+            //for quad-/Oct-tree like grids, it is safter to use FillPatchNLevels
+            Vector<PhysBCFunct<GpuBndryFuncFab<IncfloTracFill>>> physbcs;
+            for (int ilev = 0; ilev <= finest_level; ++ilev) {
+                physbcs.emplace_back(geom[ilev],bcrec,IncfloTracFill{m_probtype, m_ntrac, m_bc_tracer_d, m_bc_velocity});
+            }
+            Vector<Vector<MultiFab*>> smf(finest_level+1);
+            Vector<Vector<Real>> st(finest_level+1);
+            for (int ilev = 0; ilev <= finest_level; ++ilev) {
+              smf[ilev] = {&(m_leveldata[ilev]->tracer_o), &(m_leveldata[ilev]->tracer)};
+              st[ilev] = {m_t_old[ilev], m_t_new[ilev]};
+              //smf[ilev].push_back(&(m_leveldata[ilev]->tracer_o));
+              //smf[ilev].push_back(&(m_leveldata[ilev]->tracer));
+            }
+            FillPatchNLevels(tracer, lev, IntVect(ng), time, smf, st, 0, 0, m_ntrac, geom,
+                             physbcs, 0, ref_ratio, mapper, bcrec, 0);
+        }
     }
 }
 
@@ -134,12 +187,30 @@ void incflo::fillpatch_gradp (int lev, Real time, MultiFab& gp, int ng)
 #else
         Interpolater* mapper = &cell_cons_interp;
 #endif
-        FillPatchTwoLevels(gp, IntVect(ng), time,
-                           {&(m_leveldata[lev-1]->gp)}, {time},
-                           {&(m_leveldata[lev]->gp)}, {time},
-                           0, 0, AMREX_SPACEDIM, geom[lev-1], geom[lev],
-                           cphysbc, 0, fphysbc, 0,
-                           refRatio(lev-1), mapper, bcrec, 0);
+
+        if (!m_fillpatchnlevels){
+          FillPatchTwoLevels(gp, IntVect(ng), time,
+                             {&(m_leveldata[lev-1]->gp)}, {time},
+                             {&(m_leveldata[lev]->gp)}, {time},
+                             0, 0, AMREX_SPACEDIM, geom[lev-1], geom[lev],
+                             cphysbc, 0, fphysbc, 0,
+                             refRatio(lev-1), mapper, bcrec, 0);
+        }else{
+            //for quad-/Oct-tree like grids, it is safter to use FillPatchNLevels
+            Vector<PhysBCFunct<GpuBndryFuncFab<IncfloForFill>>> physbcs;
+            for (int ilev = 0; ilev <= finest_level; ++ilev) {
+                physbcs.emplace_back(geom[ilev],bcrec,IncfloForFill{m_probtype});
+            }
+            Vector<Vector<MultiFab*>> smf(finest_level+1);
+            Vector<Vector<Real>> st(finest_level+1);
+            for (int ilev = 0; ilev <= finest_level; ++ilev) {
+              smf[ilev] = {&(m_leveldata[ilev]->gp)};
+              st[ilev] = {time};
+            }
+            FillPatchNLevels(gp, lev, IntVect(ng), time, smf, st, 0, 0, AMREX_SPACEDIM, geom,
+                             physbcs, 0, ref_ratio, mapper, bcrec, 0);
+        }
+
     }
 }
 
@@ -163,12 +234,29 @@ void incflo::fillpatch_force (Real time, Vector<MultiFab*> const& force, int ng)
         PhysBCFunct<GpuBndryFuncFab<IncfloForFill> > fphysbc
             (geom[lev  ], bcrec, IncfloForFill{m_probtype});
         Interpolater* mapper = &pc_interp;
-        FillPatchTwoLevels(*force[lev], IntVect(ng), time,
-                           {force[lev-1]}, {time},
-                           {force[lev  ]}, {time},
-                           0, 0, ncomp, geom[lev-1], geom[lev],
-                           cphysbc, 0, fphysbc, 0,
-                           refRatio(lev-1), mapper, bcrec, 0);
+
+        if (!m_fillpatchnlevels){
+          FillPatchTwoLevels(*force[lev], IntVect(ng), time,
+                             {force[lev-1]}, {time},
+                             {force[lev  ]}, {time},
+                             0, 0, ncomp, geom[lev-1], geom[lev],
+                             cphysbc, 0, fphysbc, 0,
+                             refRatio(lev-1), mapper, bcrec, 0);
+        }else{
+            //for quad-/Oct-tree like grids, it is safter to use FillPatchNLevels
+            Vector<PhysBCFunct<GpuBndryFuncFab<IncfloForFill>>> physbcs;
+            for (int ilev = 0; ilev <= finest_level; ++ilev) {
+                physbcs.emplace_back(geom[ilev],bcrec,IncfloForFill{m_probtype});
+            }
+            Vector<Vector<MultiFab*>> smf(finest_level+1);
+            Vector<Vector<Real>> st(finest_level+1);
+            for (int ilev = 0; ilev <= finest_level; ++ilev) {
+              smf[ilev] = {force[lev-1]};
+              st[ilev] = {time};
+            }
+            FillPatchNLevels(*force[lev], lev, IntVect(ng), time, smf, st, 0, 0, ncomp, geom,
+                             physbcs, 0, ref_ratio, mapper, bcrec, 0);
+        }
     }
 }
 

src/convection/incflo_compute_MAC_projected_velocities.cpp

@@ -98,7 +98,7 @@ incflo::compute_MAC_projected_velocities (
         LPInfo lp_info;
         lp_info.setMaxCoarseningLevel(m_mac_mg_max_coarsening_level);
 #ifndef AMREX_USE_EB
-        if (m_constant_density) {
+        if (m_constant_density&&!m_vof_advect_tracer) {
             Vector<BoxArray> ba;
             Vector<DistributionMapping> dm;
             for (auto const& ir : inv_rho) {
@@ -123,7 +123,7 @@ incflo::compute_MAC_projected_velocities (
         }
     } else {
 #ifndef AMREX_USE_EB
-        if (m_constant_density) {
+        if (m_constant_density&&!m_vof_advect_tracer) {
             macproj->updateBeta(l_dt/m_ro_0);  // unnecessary unless m_ro_0 changes.
         } else
 #endif
@@ -176,6 +176,52 @@ incflo::compute_MAC_projected_velocities (
                                       m_godunov_ppm, m_godunov_use_forces_in_trans,
                                       l_advection_type, PPM::default_limiter,
                                       allow_inflow_on_outflow, BC_MF.get());
+
+        //add surface tension
+        if(m_vof_advect_tracer && m_use_cc_proj)
+          get_volume_of_fluid ()->velocity_face_source(lev,0.5*l_dt, AMREX_D_DECL(*u_mac[lev], *v_mac[lev], *w_mac[lev]),
+                                                       AMREX_D_DECL(nullptr, nullptr, nullptr));
+
+
+  if(0){
+     //The following is only used for testing the pure advection of VOF algorithm
+     //Average the cell-centered velocity to face center as MAC velocity
+     #ifdef _OPENMP
+     #pragma omp parallel if (Gpu::notInLaunchRegion())
+     #endif
+       for (MFIter mfi(*vel[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
+       {
+         // Note nodaltilebox will not include the nodal index beyond boundaries between neighboring
+         // titles. Therefore,if we want to use face values (i.e., face_val) immediately below (commented
+         // out), we must create index space for the face-centered values of the tiled region
+         // (i.e., surroundingNodes()).
+         Box const& bx  = mfi.tilebox();
+         AMREX_D_TERM(Box const& xbx = mfi.nodaltilebox(0);,
+                      Box const& ybx = mfi.nodaltilebox(1);,
+                      Box const& zbx = mfi.nodaltilebox(2););
+         Array4<Real const> const& velocity   = vel[lev]->const_array(mfi);
+         AMREX_D_TERM(Array4<Real > const& xfv = u_mac[lev]->array(mfi);,
+                      Array4<Real > const& yfv = v_mac[lev]->array(mfi);,
+                      Array4<Real > const& zfv = w_mac[lev]->array(mfi););
+         AMREX_D_TERM(
+            ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            {
+              xfv(i,j,k) = .5*(velocity(i,j,k,0)+velocity(i-1,j,k,0));
+            });,
+            ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            {
+               yfv(i,j,k) = .5*(velocity(i,j,k,1)+velocity(i,j-1,k,1));
+            });,
+
+            ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            {
+               zfv(i,j,k) = .5*(velocity(i,j,k,2)+velocity(i,j,k-1,2));
+            });
+          )    // end AMREX_D_TERM
+
+        }
+        if (lev == finest_level) return;
+     }//test
     }
 
     Vector<Array<MultiFab*,AMREX_SPACEDIM> > mac_vec(finest_level+1);