[pull] master from Exawind:master

include/SolutionOptions.h

@@ -102,6 +102,8 @@ class SolutionOptions
 
   bool use_balanced_buoyancy_force_{false};
   bool realm_has_vof_{false};
+  double vof_sharpening_scaling_factor_{5.0};
+  double vof_diffusion_scaling_factor_{0.6};
 
   double hybridDefault_;
   double alphaDefault_;

src/SolutionOptions.C

@@ -136,6 +136,14 @@ SolutionOptions::load(const YAML::Node& y_node)
       y_solution_options, "use_balanced_buoyancy_force",
       use_balanced_buoyancy_force_, use_balanced_buoyancy_force_);
 
+    get_if_present(
+      y_solution_options, "vof_sharpening_scaling_factor",
+      vof_sharpening_scaling_factor_, vof_sharpening_scaling_factor_);
+
+    get_if_present(
+      y_solution_options, "vof_diffusion_scaling_factor",
+      vof_diffusion_scaling_factor_, vof_diffusion_scaling_factor_);
+
     // Solve for incompressible continuity
     get_if_present(
       y_solution_options, "solve_incompressible_continuity",

src/edge_kernels/VOFAdvectionEdgeAlg.C

@@ -87,6 +87,11 @@ VOFAdvectionEdgeAlg::execute()
   const int ndim = realm_.meta_data().spatial_dimension();
   const auto& meta = realm_.meta_data();
 
+  const DblType sharpening_scaling =
+    realm_.solutionOptions_->vof_sharpening_scaling_factor_;
+  const DblType diffusion_scaling =
+    realm_.solutionOptions_->vof_diffusion_scaling_factor_;
+
   const DblType alphaUpw = realm_.get_alpha_upw_factor("volume_of_fluid");
   const DblType hoUpwind = realm_.get_upw_factor("volume_of_fluid");
   const DblType relaxFac =
@@ -214,16 +219,13 @@ VOFAdvectionEdgeAlg::execute()
         axdx += av[d] * dxj;
       }
 
-      // Hard-coded values comes from Jain, 2022 to enforce
-      // VOF function bounds of [0,1] while maintaining interface
-      // thickness that is ~2 cells
-
       const DblType velocity_scale =
-        5.0 * stk::math::abs(
-                vdot /
-                stk::math::sqrt(av[0] * av[0] + av[1] * av[1] + av[2] * av[2]));
+        sharpening_scaling *
+        stk::math::abs(
+          vdot /
+          stk::math::sqrt(av[0] * av[0] + av[1] * av[1] + av[2] * av[2]));
 
-      diffusion_coef = stk::math::sqrt(diffusion_coef) * 0.3;
+      diffusion_coef = stk::math::sqrt(diffusion_coef) * diffusion_scaling;
 
       const DblType inv_axdx = 1.0 / axdx;