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/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);

src/convection/incflo_compute_advection_term.cpp

@@ -66,7 +66,7 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
 
     // Make one flux MF at each level to hold all the fluxes (velocity, density, tracers)
     int n_flux_comp = AMREX_SPACEDIM;
-    if (!m_constant_density) n_flux_comp += 1;
+    if (!m_constant_density&&!m_update_density_from_vof) n_flux_comp += 1;
     if ( m_advect_tracer)    n_flux_comp += m_ntrac;
 
     // This will hold state on faces
@@ -139,7 +139,7 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
     //    and compute the tracer forcing terms for the first time
     if (m_advection_type != "MOL") {
 
-        compute_vel_forces(vel_forces, vel, density, tracer, tracer);
+        compute_vel_forces(vel_forces, vel, density, tracer, tracer, m_use_cc_proj?false:true);
 
         if (m_godunov_include_diff_in_forcing) {
 
@@ -259,15 +259,57 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
                 Array<PhysBCFunct<GpuBndryFuncFab<IncfloVelFill>>,AMREX_SPACEDIM>
                     cbndyFuncArr = {AMREX_D_DECL(crse_bndry_func,crse_bndry_func,crse_bndry_func)};
 
-                // Use piecewise constant interpolation in time
-                Array<int, AMREX_SPACEDIM> idx = {AMREX_D_DECL(0,1,2)};
-                FillPatchTwoLevels(u_fine, IntVect(nghost_mac()), time_nph,
-                                   {u_crse}, {time_nph},
-                                   {u_fine}, {time_nph},
-                                   0, 0, 1,
-                                   geom[lev-1], geom[lev],
-                                   cbndyFuncArr, idx, fbndyFuncArr, idx,
-                                   rr, mapper, bcrecArr, idx);
+               if (/*1*/m_fillpatch_method == 0){
+                  // Use piecewise constant interpolation in time
+                  Array<int, AMREX_SPACEDIM> idx = {AMREX_D_DECL(0,1,2)};
+                  FillPatchTwoLevels(u_fine, IntVect(nghost_mac()), time_nph,
+                                     {u_crse}, {time_nph},
+                                     {u_fine}, {time_nph},
+                                     0, 0, 1,
+                                     geom[lev-1], geom[lev],
+                                     cbndyFuncArr, idx, fbndyFuncArr, idx,
+                                     rr, mapper, bcrecArr, idx);
+               }else{
+                  mapper = &face_cons_linear_interp;
+                  for (int dir=0;dir<AMREX_SPACEDIM;++dir){
+                    if (m_fillpatch_method == 1){
+                      // This FillPatch operation interpolates using the ghost cells of the coarser level
+                      // via `PhysBCFunctUseCoarseGhost`, which is defined in `AMReX_PhysBCFunct.h`.
+                      // For implementation details, see `AMReX_FillPatchUtil_I.h`.
+                      //
+                      // When the `blocking_factor` is small (e.g., 1, 2, or 4), specifically used for generating
+                      // quad-/octree-like grids, this FillPatch method is necessary instead of the previous one.
+                      FillPatchTwoLevels (*u_fine[dir], IntVect(nghost_mac()), IntVect (0), time_nph,
+                                          {u_crse[dir]}, {time_nph}, {u_fine[dir]}, {time_nph},
+                                          0, 0, 1, geom[lev-1], geom[lev],
+                                          refRatio(lev-1), mapper, bcrecArr[dir], 0);
+                      //The physical boundary condition is not enforced in the above fillpatch, so we have to do it here.
+                      fbndyFuncArr[dir].FillBoundary(*u_fine[dir], 0, 1, IntVect(nghost_mac()), time_nph, 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],m_bcrec_velocity,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] = dir < 1 ? Vector<MultiFab*>{u_mac[ilev]} :
+#if (AMREX_SPACEDIM == 3)
+                                    dir > 1 ? Vector<MultiFab*>{w_mac[ilev]} :
+#endif
+                                            Vector<MultiFab*>{v_mac[ilev]};
+                        st[ilev] = {time_nph};
+                      }
+                      FillPatchNLevels(*u_fine[dir], lev, IntVect(nghost_mac()), time_nph,
+                                       smf, st, 0, 0, 1, geom,
+                                       physbcs, 0, ref_ratio, mapper, m_bcrec_velocity, 0);
+
+                    }
+                  }
+               }
+
             }
         } // end umac fill
 
@@ -319,7 +361,7 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
                     Multiply(vel_nph, rho_nph, 0, n, 1, 1);
                 }
             }
-
+            //vel_nph.setVal(0.);
             if (m_advect_tracer && (m_ntrac>0)) {
                 trac_nph.setVal(0.);
                 fillphysbc_tracer(lev, time_nph, trac_nph, 1);
@@ -446,7 +488,7 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
             // ************************************************************************
             // Density
             // ************************************************************************
-            if (!m_constant_density)
+            if (!m_constant_density&&!m_update_density_from_vof)
             {
                 face_comp = AMREX_SPACEDIM;
                 ncomp = 1;
@@ -512,7 +554,7 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
                     });
                 }
 
-                if (m_constant_density)
+                if (m_constant_density||m_update_density_from_vof)
                    face_comp = AMREX_SPACEDIM;
                 else
                    face_comp = AMREX_SPACEDIM+1;
@@ -663,7 +705,7 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
 
         // Note: density is always updated conservatively -- we do not provide an option for
         //       updating density convectively
-        if (!m_constant_density)
+        if (!m_constant_density&&!m_update_density_from_vof)
         {
           int flux_comp = AMREX_SPACEDIM;
 //Was this OMP intentionally left off?
@@ -705,7 +747,7 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
 
         if (m_advect_tracer && m_ntrac > 0)
         {
-          int flux_comp = (m_constant_density) ? AMREX_SPACEDIM : AMREX_SPACEDIM+1;
+          int flux_comp = (m_constant_density||m_update_density_from_vof) ? AMREX_SPACEDIM : AMREX_SPACEDIM+1;
 
 #ifdef _OPENMP
 #pragma omp parallel if (Gpu::notInLaunchRegion())
@@ -783,7 +825,7 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
                               bc_vel, lev);
 
             // density
-            if (!m_constant_density) {
+            if (!m_constant_density&&!m_update_density_from_vof) {
                 auto const& bc_den = get_density_bcrec_device_ptr();
                 redistribute_term(mfi, *conv_r[lev], drdt_tmp,
                                   *density[lev], bc_den, lev);

src/diffusion/DiffusionTensorOp.cpp

@@ -74,6 +74,7 @@ DiffusionTensorOp::DiffusionTensorOp (incflo* a_incflo)
             m_reg_apply_op->setMaxOrder(m_mg_maxorder);
             m_reg_apply_op->setDomainBC(m_incflo->get_diffuse_tensor_bc(Orientation::low),
                                         m_incflo->get_diffuse_tensor_bc(Orientation::high));
+            //m_reg_apply_op->setInterpBndryHalfWidth(1);
         }
     }
 }
@@ -228,9 +229,9 @@ void DiffusionTensorOp::compute_divtau (Vector<MultiFab*> const& a_divtau,
     for (int lev = 0; lev <= finest_level; ++lev) {
         velocity[lev].define(a_velocity[lev]->boxArray(),
                              a_velocity[lev]->DistributionMap(),
-                             AMREX_SPACEDIM, 1, MFInfo(),
+                             AMREX_SPACEDIM, 1/*a_velocity[lev]->nGrow()*/, MFInfo(),
                              a_velocity[lev]->Factory());
-        MultiFab::Copy(velocity[lev], *a_velocity[lev], 0, 0, AMREX_SPACEDIM, 1);
+        MultiFab::Copy(velocity[lev], *a_velocity[lev], 0, 0, AMREX_SPACEDIM, 1/*velocity[lev].nGrow()*/);
     }
 
 #ifdef AMREX_USE_EB
@@ -293,7 +294,32 @@ void DiffusionTensorOp::compute_divtau (Vector<MultiFab*> const& a_divtau,
         MLMG mlmg(*m_reg_apply_op);
         mlmg.apply(a_divtau, GetVecOfPtrs(velocity));
     }
-
+//for (int lev = finest_level; lev >= 0; --lev) {
+//
+//    for (MFIter mfi(*a_divtau[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi) {
+//        const Box& bx = mfi./*tilebox()*/growntilebox(velocity[lev].nGrow());
+//        const Box& bxv = mfi.tilebox();
+//        const auto lo = amrex::lbound(bxv);
+//        const auto hi = amrex::ubound(bxv);
+//        const auto& arr = (*a_divtau[lev])[mfi].array();
+//        const auto& arr0 = (velocity[lev])[mfi].array();
+//        const auto& arr1 = (*a_velocity[lev])[mfi].array();
+//
+//        // Loop over the valid region
+//        if (lo.x==10 && lo.y==16)
+//        amrex::ParallelFor(bx, velocity[lev].nComp(), [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) {
+//            //if (/*std::isnan(arr0(i, j, k, n))&&*/) {
+//                amrex::Print() << "NaN found at index (" << i << ", " << j << ", " << k
+//                               << ") in component " << n<<"  "<<arr0(i,j,k,n)<<" ,"<<arr1(i,j,k,n) << "\n";
+//            //}
+//        });
+//    }
+//
+//
+//if (a_divtau[lev]->contains_nan(0, AMREX_SPACEDIM,0)) {
+//    amrex::Abort("NaN detected in MultiFab!");
+//}
+//}
     bool advect_momentum = m_incflo->AdvectMomentum();
     if (!advect_momentum) {
 #ifdef _OPENMP