From 74ef02f645a4798f3570595563b3593889d5c168 Mon Sep 17 00:00:00 2001
From: psakievich <psakiev@sandia.gov>
Date: Mon, 28 Aug 2023 22:04:46 -0400
Subject: [PATCH] Revert "Change VOF density flux form in mdot to cancel out
 buoyancy (#1192)"

This reverts commit 923bc265039eaf7a744dbfe146de6a08a9f78523.
---
 .../VOFDroplet/VOFDroplet.norm.gold           |  6 ++---
 .../test_files/VOFDroplet/VOFDroplet.yaml     |  2 +-
 src/edge_kernels/ContinuityEdgeSolverAlg.C    | 23 +++++++++++--------
 src/user_functions/DropletVOFAuxFunction.C    |  5 ++--
 4 files changed, 20 insertions(+), 16 deletions(-)

diff --git a/reg_tests/test_files/VOFDroplet/VOFDroplet.norm.gold b/reg_tests/test_files/VOFDroplet/VOFDroplet.norm.gold
index c69bff960..d58c7fdd0 100644
--- a/reg_tests/test_files/VOFDroplet/VOFDroplet.norm.gold
+++ b/reg_tests/test_files/VOFDroplet/VOFDroplet.norm.gold
@@ -1,3 +1,3 @@
-2.67343413675555e-05 11 0.0001
-6.859619665424383e-06 12 0.0002
-3.313464463487743e-06 13 0.0003
+0.0001417205047386008 11 0.0001
+4.686504088256447e-05 12 0.0002
+1.667681063887448e-05 13 0.0003
diff --git a/reg_tests/test_files/VOFDroplet/VOFDroplet.yaml b/reg_tests/test_files/VOFDroplet/VOFDroplet.yaml
index b3e88225c..9c6823b22 100755
--- a/reg_tests/test_files/VOFDroplet/VOFDroplet.yaml
+++ b/reg_tests/test_files/VOFDroplet/VOFDroplet.yaml
@@ -29,7 +29,7 @@ linear_solvers:
 realms:
 
   - name: realm_1
-    mesh: "generated:75x75x75|bbox:-0.5,-0.5,-0.5,0.5,0.5,0.5|sideset:xXyYzZ|show" 
+    mesh: "generated:150x150x4|bbox:-0.5,-0.5,-0.5,0.5,0.5,0.5|sideset:xXyYzZ|show" 
     use_edges: yes
     automatic_decomposition_type: rcb
 
diff --git a/src/edge_kernels/ContinuityEdgeSolverAlg.C b/src/edge_kernels/ContinuityEdgeSolverAlg.C
index 5507ec8c9..62188924c 100644
--- a/src/edge_kernels/ContinuityEdgeSolverAlg.C
+++ b/src/edge_kernels/ContinuityEdgeSolverAlg.C
@@ -88,15 +88,22 @@ ContinuityEdgeSolverAlg::execute()
       const DblType pressureL = pressure.get(nodeL, 0);
       const DblType pressureR = pressure.get(nodeR, 0);
 
-      const DblType densityL = density.get(nodeL, 0);
-      const DblType densityR = density.get(nodeR, 0);
+      const DblType densityL =
+        om_solveIncompressibleEqn * density.get(nodeL, 0) +
+        solveIncompressibleEqn;
+      const DblType densityR =
+        om_solveIncompressibleEqn * density.get(nodeR, 0) +
+        solveIncompressibleEqn;
+
+      const DblType udiagL =
+        udiag.get(nodeL, 0) * (om_solveIncompressibleEqn +
+                               solveIncompressibleEqn * density.get(nodeL, 0));
+      const DblType udiagR =
+        udiag.get(nodeR, 0) * (om_solveIncompressibleEqn +
+                               solveIncompressibleEqn * density.get(nodeR, 0));
 
-      const DblType udiagL = udiag.get(nodeL, 0);
-      const DblType udiagR = udiag.get(nodeR, 0);
       const DblType projTimeScale = 0.5 * (1.0 / udiagL + 1.0 / udiagR);
       const DblType rhoIp = 0.5 * (densityL + densityR);
-      const DblType denScale =
-        (1.0 / rhoIp) * solveIncompressibleEqn + om_solveIncompressibleEqn;
 
       DblType axdx = 0.0;
       DblType asq = 0.0;
@@ -130,9 +137,7 @@ ContinuityEdgeSolverAlg::execute()
           kxj * GjIp * nocFac;
       }
       tmdot /= tauScale;
-      tmdot *= denScale;
-      const DblType lhsfac =
-        -asq * inv_axdx * projTimeScale * denScale / tauScale;
+      const DblType lhsfac = -asq * inv_axdx * projTimeScale / tauScale;
 
       // Left node entries
       smdata.lhs(0, 0) = -lhsfac;
diff --git a/src/user_functions/DropletVOFAuxFunction.C b/src/user_functions/DropletVOFAuxFunction.C
index dd075daec..2e66ca446 100644
--- a/src/user_functions/DropletVOFAuxFunction.C
+++ b/src/user_functions/DropletVOFAuxFunction.C
@@ -38,13 +38,12 @@ DropletVOFAuxFunction::do_evaluate(
 
     const double x = coords[0];
     const double y = coords[1];
-    const double z = coords[2];
-    const double interface_thickness = 0.025;
+    const double interface_thickness = 0.005;
 
     fieldPtr[0] = 0.0;
     fieldPtr[0] += -0.5 * (std::erf(y / interface_thickness) + 1.0) + 1.0;
 
-    auto radius = std::sqrt(x * x + (y - 0.25) * (y - 0.25) + z * z) - 0.075;
+    auto radius = std::sqrt(x * x + (y - 0.15) * (y - 0.15)) - 0.075;
     fieldPtr[0] += -0.5 * (std::erf(radius / interface_thickness) + 1.0) + 1.0;
 
     fieldPtr += fieldSize;