VOF BC changes to address sticking fluid (#1354)

* New VOF BC

* Formatting

---------

Co-authored-by: whorne <whorne@sandia.gov>

Author: wjhorne

src/LowMachEquationSystem.C

@@ -1958,14 +1958,6 @@ MomentumEquationSystem::register_wall_bc(
     stk::mesh::put_field_on_mesh(
       *wallNormalDistanceBip, *part, numScsBip, nullptr);
 
-    if (realm_.solutionOptions_->use_balanced_buoyancy_force_) {
-      if (!buoyancySrcMask_) {
-        buoyancySrcMask_.reset(new NodalBuoyancyFuncUtil(realm_, part));
-      } else {
-        buoyancySrcMask_->partVec_.push_back(part);
-      }
-    }
-
     // need wall friction velocity for TKE boundary condition
     if (RANSAblBcApproach_) {
       const AlgorithmType wfAlgType = WALL_FCN;
@@ -2105,6 +2097,14 @@ MomentumEquationSystem::register_wall_bc(
     }
   }
 
+  if (realm_.solutionOptions_->use_balanced_buoyancy_force_) {
+    if (!buoyancySrcMask_) {
+      buoyancySrcMask_.reset(new NodalBuoyancyFuncUtil(realm_, part));
+    } else {
+      buoyancySrcMask_->partVec_.push_back(part);
+    }
+  }
+
   // Dirichlet wall boundary condition.
   if (!anyWallFunctionActivated || userData.isNoSlip_) {
     const AlgorithmType algType = WALL;

src/edge_kernels/VOFAdvectionEdgeAlg.C

@@ -118,6 +118,13 @@ VOFAdvectionEdgeAlg::execute()
   const auto velocity_rtm = fieldMgr.get_field<double>(
     get_field_ordinal(realm_.meta_data(), velocity_rtm_name));
 
+  const bool using_balanced_force =
+    realm_.solutionOptions_->use_balanced_buoyancy_force_;
+  const std::string wall_mask_name =
+    using_balanced_force ? "buoyancy_source_mask" : "density";
+  const auto wall_mask = fieldMgr.get_field<double>(
+    get_field_ordinal(realm_.meta_data(), wall_mask_name));
+
   run_algorithm(
     realm_.bulk_data(),
     KOKKOS_LAMBDA(
@@ -192,7 +199,6 @@ VOFAdvectionEdgeAlg::execute()
       smdata.lhs(0, 1) += alhsfac;
 
       // Compression term
-
       DblType dOmegadxMag = 0.0;
       DblType interface_gradient[3] = {0.0, 0.0, 0.0};
 
@@ -207,9 +213,15 @@ VOFAdvectionEdgeAlg::execute()
 
       dOmegadxMag = stk::math::sqrt(dOmegadxMag);
 
+      const DblType left_mask =
+        using_balanced_force ? wall_mask.get(nodeL, 0) : 1.0;
+      const DblType right_mask =
+        using_balanced_force ? wall_mask.get(nodeR, 0) : 1.0;
+
       // No gradient == no interface
-      if (dOmegadxMag < gradient_eps)
+      if (dOmegadxMag < gradient_eps) {
         return;
+      }
 
       for (int d = 0; d < ndim; ++d)
         interface_normal[d] = interface_gradient[d] / dOmegadxMag;
@@ -236,25 +248,28 @@ VOFAdvectionEdgeAlg::execute()
       local_velocity += vdot;
       local_velocity = stk::math::abs(local_velocity) / face_area;
 
-      const DblType velocity_scale = sharpening_scaling * local_velocity;
+      const DblType velocity_scale =
+        sharpening_scaling * local_velocity * left_mask * right_mask;
 
       diffusion_coef = stk::math::sqrt(diffusion_coef) * diffusion_scaling;
 
       const DblType inv_axdx = 1.0 / axdx;
 
       const DblType dlhsfac = -velocity_scale * diffusion_coef * asq * inv_axdx;
 
-      smdata.rhs(0) -= dlhsfac * (qNp1R - qNp1L);
-      smdata.rhs(1) += dlhsfac * (qNp1R - qNp1L);
+      smdata.rhs(0) -= dlhsfac * (qNp1R - qNp1L) +
+                       inv_axdx * (1.0 - left_mask) * (qNp1R - qNp1L);
+      smdata.rhs(1) += dlhsfac * (qNp1R - qNp1L) +
+                       inv_axdx * (1.0 - right_mask) * (qNp1R - qNp1L);
 
       massVofBalancedFlowRate.get(edge, 0) =
         dlhsfac * (qNp1R - qNp1L) * (density_liquid - density_gas);
 
-      smdata.lhs(0, 0) -= dlhsfac;
-      smdata.lhs(0, 1) += dlhsfac;
+      smdata.lhs(0, 0) -= dlhsfac + inv_axdx * (1.0 - left_mask);
+      smdata.lhs(0, 1) += dlhsfac + inv_axdx * (1.0 - left_mask);
 
-      smdata.lhs(1, 0) += dlhsfac;
-      smdata.lhs(1, 1) -= dlhsfac;
+      smdata.lhs(1, 0) += dlhsfac + inv_axdx * (1.0 - right_mask);
+      smdata.lhs(1, 1) -= dlhsfac + inv_axdx * (1.0 - right_mask);
 
       const DblType omegaL =
         diffusion_coef * stk::math::log((qNp1L + eps) / (1.0 - qNp1L + eps));
@@ -296,6 +311,8 @@ VOFAdvectionEdgeAlg::execute()
                    stk::math::tanh(0.5 * omegaIp / diffusion_coef)) *
           interface_normal[d] * av[d];
 
+      compression = compression * left_mask * right_mask;
+
       smdata.rhs(0) -= compression;
       smdata.rhs(1) += compression;