Communication for VOF conservation (Exawind#602)

* implementing communication of the VOF field between stages of the directionally split scheme. This is crucial for conservation across processor and periodic boundaries.

* changes to output momentum and volume fraction conservation when MultiPhase is verbose.

* a fix to compare VOF to tiny - avoids instability that can occur with geometric operations on very small VOF values.

* removing an unnecessary fillpatch

* Allow user to specify whether VOF debris should be removed. Default is on, off is intended for conservation tests.

* adding unit test to check vof conservation in multiple directions, as well as accuracy in 1D flow.

* fixing a conditional

* fixing initialization loop for umac and eliminating an unneeded line

* addressing some of the CI feedback

* addressing more of the CI feedback - namely, putting the google test checks in ordinary loops instead of parallelfor.

* small edits to unit test

* making velocity initialization array a GpuArray

* eliminating unnecessary scalar equation loop and moving initialization outside of time iteration loop.

* another initialization for icheck

* setting up variable nx so that lambda function doesn't reference private variable

* removing lambda functions from protected class variables in unit tests.

amr-wind/equation_systems/vof/SplitAdvection.H

@@ -6,11 +6,12 @@
 namespace amr_wind {
 namespace multiphase {
 
-void split_advection(
+void split_advection_step(
     int lev,
     amrex::Box const& bx,
     int isweep,
     amrex::Array4<amrex::Real> const& volfrac,
+    amrex::Array4<amrex::Real> const& fluxC,
     amrex::Array4<amrex::Real const> const& umac,
     amrex::Array4<amrex::Real const> const& vmac,
     amrex::Array4<amrex::Real const> const& wmac,
@@ -20,6 +21,15 @@ void split_advection(
     amrex::Real dt,
     bool is_lagrangian);
 
+void cmask_loop(
+    amrex::Box const& bx,
+    amrex::Array4<amrex::Real> const& volfrac,
+    amrex::Array4<amrex::Real> const& fluxC,
+    bool is_lagrangian);
+
+void debris_loop(
+    amrex::Box const& bx, amrex::Array4<amrex::Real> const& volfrac);
+
 void sweep(
     const int dir,
     amrex::Box const& bx,

amr-wind/equation_systems/vof/SplitAdvection.cpp

@@ -7,11 +7,12 @@ using namespace amrex;
 
 namespace amr_wind {
 
-void multiphase::split_advection(
+void multiphase::split_advection_step(
     int lev,
     amrex::Box const& bx,
     int isweep,
     amrex::Array4<amrex::Real> const& volfrac,
+    amrex::Array4<amrex::Real> const& fluxC,
     amrex::Array4<amrex::Real const> const& umac,
     amrex::Array4<amrex::Real const> const& vmac,
     amrex::Array4<amrex::Real const> const& wmac,
@@ -37,49 +38,40 @@ void multiphase::split_advection(
 
     Array4<Real> fluxL = makeArray4(p, bxg1, 1);
     p += fluxL.size();
-    Array4<Real> fluxC = makeArray4(p, bxg1, 1);
-    p += fluxC.size();
     Array4<Real> fluxR = makeArray4(p, bxg1, 1);
 
-    if (!is_lagrangian) {
-        amrex::ParallelFor(
-            bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
-                multiphase::c_mask(i, j, k, volfrac, fluxC);
-            });
-    }
-
     if (isweep % 3 == 0) {
         sweep(
             2, bx, dtdz, wmac, volfrac, fluxL, fluxC, fluxR, pbc, domlo.z,
             domhi.z, is_lagrangian);
-        sweep(
-            0, bx, dtdx, umac, volfrac, fluxL, fluxC, fluxR, pbc, domlo.x,
-            domhi.x, is_lagrangian);
-        sweep(
-            1, bx, dtdy, vmac, volfrac, fluxL, fluxC, fluxR, pbc, domlo.y,
-            domhi.y, is_lagrangian);
     } else if (isweep % 3 == 1) {
         sweep(
             1, bx, dtdy, vmac, volfrac, fluxL, fluxC, fluxR, pbc, domlo.y,
             domhi.y, is_lagrangian);
-        sweep(
-            2, bx, dtdz, wmac, volfrac, fluxL, fluxC, fluxR, pbc, domlo.z,
-            domhi.z, is_lagrangian);
-        sweep(
-            0, bx, dtdx, umac, volfrac, fluxL, fluxC, fluxR, pbc, domlo.x,
-            domhi.x, is_lagrangian);
     } else {
         sweep(
             0, bx, dtdx, umac, volfrac, fluxL, fluxC, fluxR, pbc, domlo.x,
             domhi.x, is_lagrangian);
-        sweep(
-            1, bx, dtdy, vmac, volfrac, fluxL, fluxC, fluxR, pbc, domlo.y,
-            domhi.y, is_lagrangian);
-        sweep(
-            2, bx, dtdz, wmac, volfrac, fluxL, fluxC, fluxR, pbc, domlo.z,
-            domhi.z, is_lagrangian);
     }
+}
+
+void multiphase::cmask_loop(
+    amrex::Box const& bx,
+    amrex::Array4<amrex::Real> const& volfrac,
+    amrex::Array4<amrex::Real> const& fluxC,
+    bool is_lagrangian)
+{
+    if (!is_lagrangian) {
+        amrex::ParallelFor(
+            bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+                multiphase::c_mask(i, j, k, volfrac, fluxC);
+            });
+    }
+}
 
+void multiphase::debris_loop(
+    amrex::Box const& bx, amrex::Array4<amrex::Real> const& volfrac)
+{
     // Remove VOF debris
     amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
         multiphase::remove_vof_debris(i, j, k, volfrac);

amr-wind/equation_systems/vof/split_advection.H

@@ -38,7 +38,7 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void lagrangian_explicit(
         fluxL(i, j, k) = amrex::max(-aL, 0.0);
         fluxC(i, j, k) = 1.0 - amrex::max(aL, 0.0) + amrex::min(aR, 0.0);
         fluxR(i, j, k) = amrex::max(aR, 0.0);
-    } else if (volfrac(i, j, k) > 0.0) {
+    } else if (volfrac(i, j, k) > tiny) {
         fit_plane(i, j, k, volfrac, mx, my, mz, alpha);
         mmL = amrex::max(aL, 0.0);
         mmR = 1.0 - mmL + amrex::min(0.0, aR);
@@ -102,7 +102,7 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void eulerian_implicit(
     if (std::abs(volfrac(i, j, k) - 1.0) <= tiny) {
         fluxL(i, j, k) = amrex::max(-aL, 0.0);
         fluxR(i, j, k) = amrex::max(aR, 0.0);
-    } else if (volfrac(i, j, k) > 0.0) {
+    } else if (volfrac(i, j, k) > tiny) {
         fit_plane(i, j, k, volfrac, mx, my, mz, alpha);
         // Eulerian advection
         x0 = 0.0;

amr-wind/equation_systems/vof/vof_advection.H

@@ -23,6 +23,7 @@ struct AdvectionOp<VOF, fvm::Godunov>
     {
         amrex::ParmParse pp_multiphase("VOF");
         pp_multiphase.query("use_lagrangian", m_use_lagrangian);
+        pp_multiphase.query("remove_debris", m_rm_debris);
     }
 
     void preadvect(const FieldState /*unused*/, const amrex::Real /*unused*/) {}
@@ -41,14 +42,15 @@ struct AdvectionOp<VOF, fvm::Godunov>
         // Implicit Eulerian Sweeping method with PLIC reconstruction
         //
 
+        // Scratch field for fluxC
+        auto fluxC = repo.create_scratch_field(1, 0, amr_wind::FieldLoc::CELL);
+
         // Define the sweep time
         isweep += 1;
         if (isweep > 3) {
             isweep = 1;
         }
 
-        dof_field.fillpatch(0.0);
-
         for (int lev = 0; lev < repo.num_active_levels(); ++lev) {
             amrex::MFItInfo mfi_info;
             if (amrex::Gpu::notInLaunchRegion()) {
@@ -61,15 +63,57 @@ struct AdvectionOp<VOF, fvm::Godunov>
             for (amrex::MFIter mfi(dof_field(lev), mfi_info); mfi.isValid();
                  ++mfi) {
                 const auto& bx = mfi.tilebox();
-                amrex::FArrayBox tmpfab(amrex::grow(bx, 1), 3 * VOF::ndim);
+                amrex::FArrayBox tmpfab(amrex::grow(bx, 1), 2 * VOF::ndim);
+                tmpfab.setVal<amrex::RunOn::Device>(0.0);
+
+                multiphase::cmask_loop(
+                    bx, dof_field(lev).array(mfi), (*fluxC)(lev).array(mfi),
+                    m_use_lagrangian);
+
+                multiphase::split_advection_step(
+                    lev, bx, isweep + 0, dof_field(lev).array(mfi),
+                    (*fluxC)(lev).array(mfi), u_mac(lev).const_array(mfi),
+                    v_mac(lev).const_array(mfi), w_mac(lev).const_array(mfi),
+                    dof_field.bcrec_device().data(), tmpfab.dataPtr(), geom, dt,
+                    m_use_lagrangian);
+
+                amrex::Gpu::streamSynchronize();
+            }
+
+            dof_field(lev).FillBoundary(geom[lev].periodicity());
+
+            for (amrex::MFIter mfi(dof_field(lev), mfi_info); mfi.isValid();
+                 ++mfi) {
+                const auto& bx = mfi.tilebox();
+                amrex::FArrayBox tmpfab(amrex::grow(bx, 1), 2 * VOF::ndim);
+                tmpfab.setVal<amrex::RunOn::Device>(0.0);
+                multiphase::split_advection_step(
+                    lev, bx, isweep + 1, dof_field(lev).array(mfi),
+                    (*fluxC)(lev).array(mfi), u_mac(lev).const_array(mfi),
+                    v_mac(lev).const_array(mfi), w_mac(lev).const_array(mfi),
+                    dof_field.bcrec_device().data(), tmpfab.dataPtr(), geom, dt,
+                    m_use_lagrangian);
+
+                amrex::Gpu::streamSynchronize();
+            }
+
+            dof_field(lev).FillBoundary(geom[lev].periodicity());
+
+            for (amrex::MFIter mfi(dof_field(lev), mfi_info); mfi.isValid();
+                 ++mfi) {
+                const auto& bx = mfi.tilebox();
+                amrex::FArrayBox tmpfab(amrex::grow(bx, 1), 2 * VOF::ndim);
                 tmpfab.setVal<amrex::RunOn::Device>(0.0);
-                multiphase::split_advection(
-                    lev, bx, isweep, dof_field(lev).array(mfi),
-                    u_mac(lev).const_array(mfi), v_mac(lev).const_array(mfi),
-                    w_mac(lev).const_array(mfi),
+                multiphase::split_advection_step(
+                    lev, bx, isweep + 2, dof_field(lev).array(mfi),
+                    (*fluxC)(lev).array(mfi), u_mac(lev).const_array(mfi),
+                    v_mac(lev).const_array(mfi), w_mac(lev).const_array(mfi),
                     dof_field.bcrec_device().data(), tmpfab.dataPtr(), geom, dt,
                     m_use_lagrangian);
 
+                if (m_rm_debris) {
+                    multiphase::debris_loop(bx, dof_field(lev).array(mfi));
+                }
                 amrex::Gpu::streamSynchronize();
             }
         }
@@ -81,6 +125,7 @@ struct AdvectionOp<VOF, fvm::Godunov>
     Field& w_mac;
     int isweep = 0;
     bool m_use_lagrangian{false};
+    bool m_rm_debris{true};
 };
 
 } // namespace pde

amr-wind/physics/multiphase/MultiPhase.H

@@ -50,6 +50,8 @@ public:
 
     amrex::Real volume_fraction_sum();
 
+    amrex::Real momentum_sum(int n);
+
     InterfaceCapturingMethod interface_capturing_method();
 
     amrex::Real rho1() const { return m_rho1; }
@@ -88,6 +90,8 @@ private:
 
     // sum of volume fractions (for vof only)
     amrex::Real m_total_volfrac{0.0};
+
+    amrex::Real q0, q1, q2, sumvof0;
 };
 
 } // namespace amr_wind

amr-wind/physics/multiphase/MultiPhase.cpp

@@ -74,6 +74,11 @@ void MultiPhase::post_init_actions()
             break;
         };
     }
+
+    q0 = momentum_sum(0);
+    q1 = momentum_sum(1);
+    q2 = momentum_sum(2);
+    sumvof0 = volume_fraction_sum();
 }
 
 void MultiPhase::pre_advance_work()
@@ -96,16 +101,23 @@ void MultiPhase::post_advance_work()
 {
     switch (m_interface_capturing_method) {
     case InterfaceCapturingMethod::VOF:
-        // Compute the print the total volume fraction
+        // Compute and print the total volume fraction, momenta, and differences
         if (m_verbose > 0) {
             m_total_volfrac = volume_fraction_sum();
+            amrex::Real mom_x = momentum_sum(0) - q0;
+            amrex::Real mom_y = momentum_sum(1) - q1;
+            amrex::Real mom_z = momentum_sum(2) - q2;
             const auto& geom = m_sim.mesh().Geom();
             const amrex::Real total_vol = geom[0].ProbDomain().volume();
             amrex::Print() << "Volume of Fluid diagnostics:" << std::endl;
-            amrex::Print() << "   Water Volume Fractions Sum : "
-                           << m_total_volfrac << std::endl;
+            amrex::Print() << "   Water Volume Fractions Sum, Difference : "
+                           << m_total_volfrac << " "
+                           << m_total_volfrac - sumvof0 << std::endl;
             amrex::Print() << "   Air Volume Fractions Sum : "
                            << total_vol - m_total_volfrac << std::endl;
+            amrex::Print() << "   Total Momentum Difference (x, y, z) : "
+                           << mom_x << " " << mom_y << " " << mom_z
+                           << std::endl;
             amrex::Print() << " " << std::endl;
         }
         break;
@@ -162,6 +174,56 @@ amrex::Real MultiPhase::volume_fraction_sum()
     return total_volume_frac;
 }
 
+amrex::Real MultiPhase::momentum_sum(int n)
+{
+    using namespace amrex;
+    BL_PROFILE("amr-wind::multiphase::ComputeVolumeFractionSum");
+    const int nlevels = m_sim.repo().num_active_levels();
+    const auto& geom = m_sim.mesh().Geom();
+    const auto& mesh = m_sim.mesh();
+
+    amrex::Real total_momentum = 0.0;
+
+    for (int lev = 0; lev < nlevels; ++lev) {
+
+        amrex::iMultiFab level_mask;
+        if (lev < nlevels - 1) {
+            level_mask = makeFineMask(
+                mesh.boxArray(lev), mesh.DistributionMap(lev),
+                mesh.boxArray(lev + 1), amrex::IntVect(2), 1, 0);
+        } else {
+            level_mask.define(
+                mesh.boxArray(lev), mesh.DistributionMap(lev), 1, 0,
+                amrex::MFInfo());
+            level_mask.setVal(1);
+        }
+
+        auto& velocity = m_sim.repo().get_field("velocity")(lev);
+        auto& density = m_sim.repo().get_field("density")(lev);
+        const amrex::Real cell_vol = geom[lev].CellSize()[0] *
+                                     geom[lev].CellSize()[1] *
+                                     geom[lev].CellSize()[2];
+
+        total_momentum += amrex::ReduceSum(
+            velocity, density, level_mask, 0,
+            [=] AMREX_GPU_HOST_DEVICE(
+                amrex::Box const& bx,
+                amrex::Array4<amrex::Real const> const& vel,
+                amrex::Array4<amrex::Real const> const& dens,
+                amrex::Array4<int const> const& mask_arr) -> amrex::Real {
+                amrex::Real vol_fab = 0.0;
+                amrex::Loop(bx, [=, &vol_fab](int i, int j, int k) noexcept {
+                    vol_fab += vel(i, j, k, n) * dens(i, j, k) *
+                               mask_arr(i, j, k) * cell_vol;
+                });
+                return vol_fab;
+            });
+    }
+    amrex::ParallelDescriptor::ReduceRealSum(total_momentum);
+
+    return total_momentum;
+}
+
 void MultiPhase::set_density_via_levelset()
 {
     const int nlevels = m_sim.repo().num_active_levels();
@@ -278,6 +340,7 @@ void MultiPhase::favre_filtering()
                 });
         }
     }
+    m_velocity.fillpatch(0.0);
 }
 
 // Reconstructing the volume fraction from a levelset function
@@ -328,6 +391,8 @@ void MultiPhase::levelset2vof()
                 });
         }
     }
+    // Fill ghost and boundary cells before simulation begins
+    (*m_vof).fillpatch(0.0);
 }
 
 } // namespace amr_wind

unit_tests/multiphase/CMakeLists.txt

@@ -1,4 +1,5 @@
 target_sources(
   ${amr_wind_unit_test_exe_name} PRIVATE
   test_vof_plic.cpp
+  test_vof_cons.cpp
   )