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
+
+        }
+        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) {
 
@@ -319,7 +319,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 +446,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 +512,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 +663,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 +705,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 +783,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/incflo.H

@@ -20,7 +20,8 @@
 
 #include <DiffusionTensorOp.H>
 #include <DiffusionScalarOp.H>
-
+#include <VolumeOfFluid.H>
+#define VOF_NODATA std::numeric_limits<Real>::max()
 enum struct StepType {
     Predictor, Corrector
 };
@@ -31,6 +32,7 @@ public:
 
     friend DiffusionTensorOp;
     friend DiffusionScalarOp;
+    friend VolumeOfFluid;
 
     enum struct FluidModel {
         Newtonian, powerlaw, Bingham, HerschelBulkley, deSouzaMendesDutra
@@ -101,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);
 
 
     ///////////////////////////////////////////////////////////////////////////
@@ -213,6 +217,19 @@ public:
                            amrex::Vector<amrex::MultiFab const*> const& eta,
                            amrex::Real dt_diff);
 
+    ///////////////////////////////////////////////////////////////////////////
+    //
+    //  tracer advection by VOF method
+    //
+    ////////////////////////////////////////////////////////////////////////////
+
+    void tracer_vof_advection (amrex::Vector<amrex::MultiFab*> const& tracer,
+                           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 (int lev, amrex::MultiFab & density, amrex::MultiFab & tracer);
+
+
     [[nodiscard]] amrex::Array<amrex::MultiFab,AMREX_SPACEDIM>
     average_scalar_eta_to_faces (int lev, int comp, amrex::MultiFab const& cc_eta) const;
 
@@ -298,11 +315,13 @@ public:
     void compute_viscosity (amrex::Vector<amrex::MultiFab*> const& eta,
                             amrex::Vector<amrex::MultiFab*> const& rho,
                             amrex::Vector<amrex::MultiFab*> const& vel,
+                            amrex::Vector<amrex::MultiFab*> const& tracer,
                             amrex::Real time, int nghost);
     void compute_viscosity_at_level (int lev,
                                      amrex::MultiFab* eta,
                                      amrex::MultiFab* rho,
                                      amrex::MultiFab* vel,
+                                     amrex::MultiFab* tracer,
                                      amrex::Geometry& lev_geom,
                                      amrex::Real time, int nghost);
     void compute_tracer_diff_coeff (amrex::Vector<amrex::MultiFab*> const& tra_eta, int nghost);
@@ -545,6 +564,18 @@ private:
     amrex::Real m_papa_reg = 0.0;
     amrex::Real m_eta_0 = 0.0;
 
+    //the number of averaging process for density
+    int m_number_of_averaging = 0;
+
+    // VOF advection parameters
+    bool m_vof_advect_tracer = false;
+    //Use VOF to update the density
+    bool m_update_density_from_vof = false;
+    // density of the phase represented by VOF=1
+    amrex::Vector<amrex::Real> m_ro_s;
+    // surface tension
+    amrex::Vector<amrex::Real> m_sigma;
+
     int m_plot_int = -1;
 
     // Dump plotfiles at as close as possible to the designated period *without* changing dt
@@ -641,6 +672,7 @@ private:
 
         // cell-centered pressure gradient
         amrex::MultiFab gp;
+        amrex::MultiFab gp_mac;
 
         amrex::MultiFab conv_velocity;
         amrex::MultiFab conv_velocity_o;
@@ -697,6 +729,8 @@ private:
     std::unique_ptr<DiffusionTensorOp> m_diffusion_tensor_op;
     std::unique_ptr<DiffusionScalarOp> m_diffusion_scalar_op;
 
+    //vof class pointer
+    std::unique_ptr<VolumeOfFluid> p_volume_of_fluid;
     //
     // end of member variables
     //
@@ -777,6 +811,8 @@ private:
     DiffusionTensorOp* get_diffusion_tensor_op ();
     DiffusionScalarOp* get_diffusion_scalar_op ();
 
+    VolumeOfFluid* get_volume_of_fluid ();
+
     amrex::Vector<amrex::MultiFab*> get_velocity_old () noexcept;
     amrex::Vector<amrex::MultiFab*> get_velocity_new () noexcept;
     amrex::Vector<amrex::MultiFab*> get_velocity_eb () noexcept;

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;
 
@@ -165,7 +166,15 @@ void incflo::Evolve()
                 printGridSummary(amrex::OutStream(), 0, finest_level);
             }
         }
-
+        if(m_vof_advect_tracer){
+        get_volume_of_fluid()->output_droplet(m_cur_time,m_nstep);
+        //if (m_nstep<10)
+        //get_volume_of_fluid()->apply_velocity_field(m_cur_time,m_nstep);
+        }
+        if (writeNow()&& m_vof_advect_tracer){
+            get_volume_of_fluid()->WriteTecPlotFile (m_cur_time,m_nstep);
+            get_volume_of_fluid()->write_tecplot_surface(m_cur_time,m_nstep);
+        }
         // Advance to time t + dt
         Advance();
         m_nstep++;
@@ -174,6 +183,8 @@ void incflo::Evolve()
         if (writeNow())
         {
             WritePlotFile();
+            //get_volume_of_fluid()->WriteTecPlotFile (m_cur_time,m_nstep);
+            //get_volume_of_fluid()->write_tecplot_surface(m_cur_time,m_nstep);
             m_last_plt = m_nstep;
         }
         if (writeNow(m_smallplot_int, m_smallplot_per_approx, -1.))

src/incflo_apply_corrector.cpp

@@ -117,7 +117,7 @@ void incflo::ApplyCorrector()
     bool include_pressure_gradient = !(m_use_mac_phi_in_godunov);
     compute_vel_forces(GetVecOfPtrs(vel_forces), get_velocity_new_const(),
                        get_density_new_const(), get_tracer_new_const(), get_tracer_new_const(),
-                       include_pressure_gradient);
+                       include_pressure_gradient, m_use_cc_proj?false:true);
     compute_MAC_projected_velocities(get_velocity_new_const(), get_density_new_const(),
                                      AMREX_D_DECL(GetVecOfPtrs(u_mac), GetVecOfPtrs(v_mac),
                                      GetVecOfPtrs(w_mac)), GetVecOfPtrs(vel_forces), new_time);
@@ -133,7 +133,7 @@ void incflo::ApplyCorrector()
     // *************************************************************************************
     // Compute viscosity / diffusive coefficients
     // *************************************************************************************
-    compute_viscosity(GetVecOfPtrs(vel_eta), get_density_new(), get_velocity_new(), new_time, 1);
+    compute_viscosity(GetVecOfPtrs(vel_eta), get_density_new(), get_velocity_new(),get_tracer_new(), new_time, 1);
 
     // Here we create divtau of the (n+1,*) state that was computed in the predictor
     if ( (m_diff_type == DiffusionType::Explicit) || use_tensor_correction )