use cc project and face-centered SF

src/convection/incflo_compute_MAC_projected_velocities.cpp

@@ -178,8 +178,9 @@ incflo::compute_MAC_projected_velocities (
                                       allow_inflow_on_outflow, BC_MF.get());
 
         //add surface tension
-        //if(m_vof_advect_tracer)
-        //  get_volume_of_fluid ()->velocity_face_source(lev,l_dt, AMREX_D_DECL(*u_mac[lev], *v_mac[lev], *w_mac[lev]));
+        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){

src/incflo.H

@@ -103,14 +103,16 @@ public:
                              amrex::Vector<amrex::MultiFab const*> const& density,
                              amrex::Vector<amrex::MultiFab const*> const& tracer_old,
                              amrex::Vector<amrex::MultiFab const*> const& tracer_new,
-                             bool include_pressure_gradient = true);
+                             bool include_pressure_gradient = true,
+                             bool include_SF = false);
     void compute_vel_forces_on_level ( int lev,
                                        amrex::MultiFab& vel_forces,
                                        const amrex::MultiFab& velocity,
                                        const amrex::MultiFab& density,
                                        const amrex::MultiFab& tracer_old,
                                        const amrex::MultiFab& tracer_new,
-                                       bool include_pressure_gradient = true);
+                                       bool include_pressure_gradient = true,
+                                       bool include_SF = false);
 
 
     ///////////////////////////////////////////////////////////////////////////
@@ -225,8 +227,8 @@ public:
                            AMREX_D_DECL(amrex::Vector<amrex::MultiFab const*> const& u_mac,
                                         amrex::Vector<amrex::MultiFab const*> const& v_mac,
                                         amrex::Vector<amrex::MultiFab const*> const& w_mac));
-    void update_vof_density (amrex::Vector<amrex::MultiFab*> const& density,
-                             amrex::Vector<amrex::MultiFab*> const& tracer);
+    void update_vof_density (int lev, amrex::Vector<amrex::MultiFab*> const& density,
+                                      amrex::Vector<amrex::MultiFab*> const& tracer);
 
 
     [[nodiscard]] amrex::Array<amrex::MultiFab,AMREX_SPACEDIM>
@@ -662,6 +664,7 @@ private:
 
         // cell-centered pressure gradient
         amrex::MultiFab gp;
+        amrex::MultiFab gp_mac;
 
         amrex::MultiFab conv_velocity;
         amrex::MultiFab conv_velocity_o;

src/incflo.cpp

@@ -82,7 +82,8 @@ void incflo::InitData ()
         }
 
         InitialIterations();
-
+//get_volume_of_fluid()->WriteTecPlotFile (m_cur_time,m_nstep);
+//amrex::Abort("finish initial projection");
         // Set m_nstep to 0 before entering time loop
         m_nstep = 0;
 

src/incflo_apply_predictor.cpp

@@ -128,7 +128,8 @@ void incflo::ApplyPredictor (bool incremental_projection)
     //when VOF method is used to advect the tracer, density and viscosity of each cell will
     // depend the VOF field value of the cell.
     if (m_vof_advect_tracer)
-        update_vof_density (get_density_old(),get_tracer_old());
+      for (int lev = 0; lev <= finest_level; ++lev)
+        update_vof_density (lev, get_density_old(),get_tracer_old());
 
     // *************************************************************************************
     // Compute explicit viscous term

src/incflo_compute_forces.cpp

@@ -40,11 +40,12 @@ void incflo::compute_vel_forces (Vector<MultiFab*> const& vel_forces,
                                  Vector<MultiFab const*> const& density,
                                  Vector<MultiFab const*> const& tracer_old,
                                  Vector<MultiFab const*> const& tracer_new,
-                                 bool include_pressure_gradient)
+                                 bool include_pressure_gradient,
+                                 bool include_SF)
 {
     for (int lev = 0; lev <= finest_level; ++lev)
        compute_vel_forces_on_level (lev, *vel_forces[lev], *velocity[lev], *density[lev],
-                                         *tracer_old[lev], *tracer_new[lev], include_pressure_gradient);
+                                         *tracer_old[lev], *tracer_new[lev], include_pressure_gradient, include_SF);
 }
 
 void incflo::compute_vel_forces_on_level (int lev,
@@ -53,7 +54,8 @@ void incflo::compute_vel_forces_on_level (int lev,
                                           const MultiFab& density,
                                           const MultiFab& tracer_old,
                                           const MultiFab& tracer_new,
-                                          bool include_pressure_gradient)
+                                          bool include_pressure_gradient,
+                                          bool include_SF)
 {
     GpuArray<Real,3> l_gravity{m_gravity[0],m_gravity[1],m_gravity[2]};
     GpuArray<Real,3> l_gp0{m_gp0[0], m_gp0[1], m_gp0[2]};
@@ -157,7 +159,7 @@ void incflo::compute_vel_forces_on_level (int lev,
     // rho:   density
 
     //fixme: we just consider the surface tension for first tracer
-    if (m_vof_advect_tracer && m_sigma[0]!=0.){
+    if (m_vof_advect_tracer && m_sigma[0]!=0.&&!m_use_cc_proj&&include_SF){
 
       VolumeOfFluid*  vof_p = get_volume_of_fluid ();
 
@@ -528,8 +530,8 @@ static int oct[3][2] = { { 1, 2 }, { 0, 2 }, { 0, 1 } };
 
     }
    }
-   else if(choice==4) {
-//cell-centered surface tension force
+      else if (choice ==4) {
+      //cell-centered surface tension force
 #ifdef _OPENMP
 #pragma omp parallel if (Gpu::notInLaunchRegion())
 #endif

src/prob/prob_init_fluid.cpp

@@ -207,6 +207,13 @@ void incflo::prob_init_fluid (int lev)
     if (1109 == m_probtype) {
         get_volume_of_fluid ()->tracer_vof_init_fraction(lev, ld.tracer);
         MultiFab::Copy(ld.tracer_o, ld.tracer, 0, 0, 1, ld.tracer.nGrow());
+        ld.tracer_o.FillBoundary(geom[lev].periodicity());
+        if (m_vof_advect_tracer){
+          update_vof_density (lev, get_density_new(),get_tracer_new());
+          MultiFab::Copy(ld.density_o, ld.density, 0, 0, 1, ld.density.nGrow());
+          fillpatch_density(lev, m_t_new[lev], ld.density_o, 3);
+        }
+
     }
 }
 

src/projection/incflo_apply_cc_projection.cpp

@@ -297,7 +297,7 @@ void incflo::ApplyCCProjection (Vector<MultiFab const*> density,
         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) {
@@ -340,13 +340,14 @@ void incflo::ApplyCCProjection (Vector<MultiFab const*> density,
                      mac_vec[lev][2] = w_mac[lev];);
     }
 
+    Vector<MultiFab> sfu_mac(finest_level+1), sfv_mac(finest_level+1), sfw_mac(finest_level+1);
     // Compute velocity on faces
     for (int lev = 0; lev <= finest_level; ++lev)
     {
         // Predict normal velocity to faces -- note that the {u_mac, v_mac, w_mac}
         //    returned from this call are on face CENTROIDS
         vel[lev]->FillBoundary(geom[lev].periodicity());
-#if 1
+#if 0
         MOL::ExtrapVelToFaces(*vel[lev],
                               AMREX_D_DECL(*u_mac[lev], *v_mac[lev], *w_mac[lev]),
                               geom[lev],
@@ -355,10 +356,27 @@ void incflo::ApplyCCProjection (Vector<MultiFab const*> density,
 #else
         average_ccvel_to_mac(             mac_vec[lev],      *vel[lev]);
 #endif
+
+
+       //add surface tension
+        AMREX_D_TERM(sfu_mac[lev].define(u_mac[lev]->boxArray(), dmap[lev], 1, u_mac[lev]->nGrow(), MFInfo(), Factory(lev));,
+                     sfv_mac[lev].define(v_mac[lev]->boxArray(), dmap[lev], 1, v_mac[lev]->nGrow(), MFInfo(), Factory(lev));,
+                     sfw_mac[lev].define(w_mac[lev]->boxArray(), dmap[lev], 1, w_mac[lev]->nGrow(), MFInfo(), Factory(lev)););
+
+        AMREX_D_TERM(sfu_mac[lev].setVal(0.0);,
+                     sfv_mac[lev].setVal(0.0);,
+                     sfw_mac[lev].setVal(0.0););
+
+
+       if(m_vof_advect_tracer)
+         get_volume_of_fluid ()->velocity_face_source(lev, scaling_factor, AMREX_D_DECL(*u_mac[lev], *v_mac[lev], *w_mac[lev]),
+                                                      AMREX_D_DECL(&sfu_mac[lev], &sfv_mac[lev], &sfw_mac[lev]));
+
     }
 
     macproj->setUMAC(mac_vec);
 
+
     if (m_verbose > 2) amrex::Print() << "CC Projection:\n";
     //
     // Perform MAC projection:  - del dot (dt/rho) grad phi = div(U)
@@ -391,12 +409,22 @@ void incflo::ApplyCCProjection (Vector<MultiFab const*> density,
 
     for (int lev=0; lev <= finest_level; ++lev)
     {
-#ifdef AMREX_USE_EB
-        amrex::Abort("Haven't written mac_to_ccvel for EB");
-#else
+//#ifdef AMREX_USE_EB
+//        amrex::Abort("Haven't written mac_to_ccvel for EB");
+//#else
         average_mac_to_ccvel(GetArrOfPtrs(m_fluxes[lev]),*cc_gphi[lev]);
-#endif
+//#endif
     }
+    // computer the cell-centered surface tension term (see note in VolumeOfFluid:: velocity_face_source)
+    VolumeOfFluid*  vof_p = get_volume_of_fluid ();
+    if(m_vof_advect_tracer)
+      for (int lev=0; lev <= finest_level; ++lev)
+      {
+        AMREX_D_TERM(Copy(m_fluxes[lev][0],sfu_mac[lev], 0, 0, 1, 0);,
+                     Copy(m_fluxes[lev][1],sfv_mac[lev], 0, 0, 1, 0);,
+                     Copy(m_fluxes[lev][2],sfw_mac[lev], 0, 0, 1, 0););
+        average_mac_to_ccvel(GetArrOfPtrs(m_fluxes[lev]),vof_p->force[lev]);
+      }
 
     for(int lev = 0; lev <= finest_level; lev++)
     {
@@ -413,55 +441,71 @@ void incflo::ApplyCCProjection (Vector<MultiFab const*> density,
 
             Array4<Real> const& u = ld.velocity.array(mfi);
             Array4<Real const> const& rho = density[lev]->const_array(mfi);
+            Array4<Real const> const& gsf = vof_p->force[lev].const_array(mfi);
 
             Real r0 = m_ro_0;
 
-            amrex::ParallelFor(tbx, [u,gphi,p_cc,phi,incremental] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
             {
                 AMREX_D_TERM(u(i,j,k,0) += gphi(i,j,k,0);,
                              u(i,j,k,1) += gphi(i,j,k,1);,
                              u(i,j,k,2) += gphi(i,j,k,2););
+                //we need to add the surface-tension effect to the cell-centered velocity
+                if(m_vof_advect_tracer){
+                  AMREX_D_TERM(u(i,j,k,0) -= gsf(i,j,k,0)*scaling_factor;,
+                               u(i,j,k,1) -= gsf(i,j,k,1)*scaling_factor;,
+                               u(i,j,k,2) -= gsf(i,j,k,2)*scaling_factor;);
+                }
                 if (incremental)
                     p_cc (i,j,k) += phi(i,j,k);
                 else
                     p_cc (i,j,k)  = phi(i,j,k);
             });
 
-            if (incremental && m_constant_density) {
+            if (incremental && m_constant_density&&!m_vof_advect_tracer) {
                 amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                 {
                     AMREX_D_TERM(gp_cc(i,j,k,0) -= gphi(i,j,k,0) * r0 / scaling_factor;,
                                  gp_cc(i,j,k,1) -= gphi(i,j,k,1) * r0 / scaling_factor;,
                                  gp_cc(i,j,k,2) -= gphi(i,j,k,2) * r0 / scaling_factor;);
 
                 });
-            } else if (incremental && !m_constant_density) {
+            } else if (incremental && (!m_constant_density||m_vof_advect_tracer)) {
                 amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                 {
                     AMREX_D_TERM(gp_cc(i,j,k,0) -= gphi(i,j,k,0) * rho(i,j,k) / scaling_factor;,
                                  gp_cc(i,j,k,1) -= gphi(i,j,k,1) * rho(i,j,k) / scaling_factor;,
                                  gp_cc(i,j,k,2) -= gphi(i,j,k,2) * rho(i,j,k) / scaling_factor;);
                 });
-            } else if (!incremental && m_constant_density) {
+            } else if (!incremental && m_constant_density&&!m_vof_advect_tracer) {
                 amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                 {
                     AMREX_D_TERM(gp_cc(i,j,k,0) = -gphi(i,j,k,0) * r0 / scaling_factor;,
                                  gp_cc(i,j,k,1) = -gphi(i,j,k,1) * r0 / scaling_factor;,
                                  gp_cc(i,j,k,2) = -gphi(i,j,k,2) * r0 / scaling_factor;);
                 });
-            } else if (!incremental && !m_constant_density) {
+            } else if (!incremental && (!m_constant_density||m_vof_advect_tracer)) {
                 amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                 {
                     AMREX_D_TERM(gp_cc(i,j,k,0) = -gphi(i,j,k,0) * rho(i,j,k) / scaling_factor;,
                                  gp_cc(i,j,k,1) = -gphi(i,j,k,1) * rho(i,j,k) / scaling_factor;,
                                  gp_cc(i,j,k,2) = -gphi(i,j,k,2) * rho(i,j,k) / scaling_factor;);
+                    if(m_vof_advect_tracer){
+                      AMREX_D_TERM(gp_cc(i,j,k,0) += gsf(i,j,k,0) * rho(i,j,k);,
+                                   gp_cc(i,j,k,1) += gsf(i,j,k,1) * rho(i,j,k);,
+                                   gp_cc(i,j,k,2) += gsf(i,j,k,2) * rho(i,j,k););
+                    }
+
                 });
             }
         }
         ld.gp.FillBoundary(geom[lev].periodicity());
         ld.p_cc.FillBoundary(geom[lev].periodicity());
     }
 
+//get_volume_of_fluid()->WriteTecPlotFile (m_cur_time,m_nstep);
+//amrex::Abort("finish initial projection");
+
     // ***************************************************************************************
     // END OF MAC STUFF
     // ***************************************************************************************

src/vof/VolumeOfFluid.H

@@ -20,9 +20,12 @@ public:
     void write_tecplot_surface(amrex::Real time, int nstep);
     void WriteTecPlotFile     (amrex::Real time, int nstep);
     void output_droplet (amrex::Real time, int nstep);
+     int domain_tag_droplets (int finest_level, amrex::Vector<amrex::BoxArray > const &grids, amrex::Vector<amrex::Geometry > const& geom,
+                              amrex::Vector<amrex::MultiFab*> const& vof,amrex::Vector<amrex::MultiFab*> const& tag);
     void apply_velocity_field(amrex::Real time, int nstep);
     void velocity_face_source(int lev,amrex::Real dt, AMREX_D_DECL(amrex::MultiFab& u_mac, amrex::MultiFab& v_mac,
-                                                   amrex::MultiFab& w_mac));
+                                                                   amrex::MultiFab& w_mac),
+                              AMREX_D_DECL(amrex::MultiFab* sfu_mac, amrex::MultiFab* sfv_mac, amrex::MultiFab* sfw_mac));
 
 // normal vector of interface
    amrex::Vector<amrex::MultiFab> normal;