Multiphase stabilization and mass-momentum consistency updates

include/LowMachEquationSystem.h

@@ -297,6 +297,7 @@ class ContinuityEquationSystem : public EquationSystem
   ScalarFieldType* pressure_;
   VectorFieldType* dpdx_;
   ScalarFieldType* massFlowRate_;
+  ScalarFieldType* massVOFBalancedFlowRate_;
   VectorFieldType* coordinates_;
 
   ScalarFieldType* pTmp_;

include/SolutionOptions.h

@@ -100,6 +100,8 @@ class SolutionOptions
 
   bool has_set_boussinesq_time_scale();
 
+  bool realm_has_vof_{false};
+
   double hybridDefault_;
   double alphaDefault_;
   double alphaUpwDefault_;

include/VolumeOfFluidEquationSystem.h

@@ -94,6 +94,8 @@ class VolumeOfFluidEquationSystem : public EquationSystem
 
   virtual void solve_and_update();
 
+  virtual void predict_state();
+
   const bool managePNG_;
   ScalarFieldType* volumeOfFluid_;
   VectorFieldType* dvolumeOfFluiddx_;

include/edge_kernels/MomentumEdgeSolverAlg.h

@@ -32,7 +32,9 @@ class MomentumEdgeSolverAlg : public AssembleEdgeSolverAlgorithm
   unsigned dudx_{stk::mesh::InvalidOrdinal};
   unsigned edgeAreaVec_{stk::mesh::InvalidOrdinal};
   unsigned massFlowRate_{stk::mesh::InvalidOrdinal};
+  unsigned massVofBalancedFlowRate_{stk::mesh::InvalidOrdinal};
   unsigned viscosity_{stk::mesh::InvalidOrdinal};
+  unsigned density_{stk::mesh::InvalidOrdinal};
   unsigned pecletFactor_{stk::mesh::InvalidOrdinal};
   unsigned maskNodeField_{stk::mesh::InvalidOrdinal};
 };

include/edge_kernels/VOFAdvectionEdgeAlg.h

@@ -39,7 +39,11 @@ class VOFAdvectionEdgeAlg : public AssembleEdgeSolverAlgorithm
   unsigned dqdx_{stk::mesh::InvalidOrdinal};
   unsigned edgeAreaVec_{stk::mesh::InvalidOrdinal};
   unsigned massFlowRate_{stk::mesh::InvalidOrdinal};
+  unsigned massVofBalancedFlowRate_{stk::mesh::InvalidOrdinal};
   unsigned density_{stk::mesh::InvalidOrdinal};
+
+  double density_liquid_;
+  double density_gas_;
 };
 
 } // namespace nalu

include/user_functions/DropletVOFAuxFunction.h

@@ -20,7 +20,7 @@ namespace nalu {
 class DropletVOFAuxFunction : public AuxFunction
 {
 public:
-  DropletVOFAuxFunction();
+  DropletVOFAuxFunction(const std::vector<double>& params);
 
   virtual ~DropletVOFAuxFunction() {}
 
@@ -34,6 +34,15 @@ class DropletVOFAuxFunction : public AuxFunction
     const unsigned fieldSize,
     const unsigned beginPos,
     const unsigned endPos) const;
+
+  int surf_idx_;
+  double droppos_x_;
+  double droppos_y_;
+  double droppos_z_;
+  double surf_pos_;
+  double surf_idx_dbl_;
+  double radius_;
+  double interface_thickness_;
 };
 
 } // namespace nalu

include/user_functions/DropletVelocityAuxFunction.h

@@ -0,0 +1,54 @@
+// Copyright 2017 National Technology & Engineering Solutions of Sandia, LLC
+// (NTESS), National Renewable Energy Laboratory, University of Texas Austin,
+// Northwest Research Associates. Under the terms of Contract DE-NA0003525
+// with NTESS, the U.S. Government retains certain rights in this software.
+//
+// This software is released under the BSD 3-clause license. See LICENSE file
+// for more details.
+//
+
+#ifndef DropletVelocityAuxFunction_h
+#define DropletVelocityAuxFunction_h
+
+#include <AuxFunction.h>
+
+#include <vector>
+
+namespace sierra {
+namespace nalu {
+
+class DropletVelocityAuxFunction : public AuxFunction
+{
+public:
+  DropletVelocityAuxFunction(
+    const unsigned beginPos,
+    const unsigned endPos,
+    const std::vector<double>& params);
+
+  virtual ~DropletVelocityAuxFunction() {}
+
+  using AuxFunction::do_evaluate;
+  virtual void do_evaluate(
+    const double* coords,
+    const double time,
+    const unsigned spatialDimension,
+    const unsigned numPoints,
+    double* fieldPtr,
+    const unsigned fieldSize,
+    const unsigned beginPos,
+    const unsigned endPos) const;
+
+  double droppos_x_;
+  double droppos_y_;
+  double droppos_z_;
+  double dropvel_x_;
+  double dropvel_y_;
+  double dropvel_z_;
+  double radius_;
+  double interface_thickness_;
+};
+
+} // namespace nalu
+} // namespace sierra
+
+#endif

include/user_functions/SloshingTankVOFAuxFunction.h

@@ -0,0 +1,47 @@
+// Copyright 2017 National Technology & Engineering Solutions of Sandia, LLC
+// (NTESS), National Renewable Energy Laboratory, University of Texas Austin,
+// Northwest Research Associates. Under the terms of Contract DE-NA0003525
+// with NTESS, the U.S. Government retains certain rights in this software.
+//
+// This software is released under the BSD 3-clause license. See LICENSE file
+// for more details.
+//
+
+#ifndef SloshingTankVOFAuxFunction_h
+#define SloshingTankVOFAuxFunction_h
+
+#include <AuxFunction.h>
+
+#include <vector>
+
+namespace sierra {
+namespace nalu {
+
+class SloshingTankVOFAuxFunction : public AuxFunction
+{
+public:
+  SloshingTankVOFAuxFunction(const std::vector<double>& params);
+
+  virtual ~SloshingTankVOFAuxFunction() {}
+
+  using AuxFunction::do_evaluate;
+  virtual void do_evaluate(
+    const double* coords,
+    const double time,
+    const unsigned spatialDimension,
+    const unsigned numPoints,
+    double* fieldPtr,
+    const unsigned fieldSize,
+    const unsigned beginPos,
+    const unsigned endPos) const;
+
+  double water_level_;
+  double amplitude_;
+  double kappa_;
+  double interface_thickness_;
+};
+
+} // namespace nalu
+} // namespace sierra
+
+#endif

reg_tests/CTestList.cmake

@@ -193,6 +193,7 @@ if(NOT ENABLE_CUDA AND NOT ENABLE_ROCM)
     add_test_r(taylorGreenVortex_p3 4)
     add_test_r(vortexOpen 4)
     add_test_r(VOFDroplet 4)
+    add_test_r(VOFInertialDroplet 4)
     add_test_r(VOFZalDisk 4)
 
     if (ENABLE_FFTW)

reg_tests/test_files/VOFDroplet/VOFDroplet.norm.gold

@@ -1,3 +1,3 @@
-2.673434136755543e-05 11 0.0001
-6.859619665424398e-06 12 0.0002
-3.313464463487801e-06 13 0.0003
+0.0006181350623480881 11 0.0001
+0.0001148227599440996 12 0.0002
+0.0001728609404847864 13 0.0003

reg_tests/test_files/VOFDroplet/VOFDroplet.yaml

@@ -58,6 +58,8 @@ realms:
         target_name: block_1
         user_function_name:
           volume_of_fluid: droplet
+        user_function_parameters:
+          volume_of_fluid: [0.0, 0.25, 0.0, 0.1, 0.0, 1, 0.0025]
       - constant: ic_2
         target_name: block_1
         value:

reg_tests/test_files/VOFInertialDroplet/VOFInertialDroplet.norm.gold

@@ -0,0 +1,3 @@
+0.0004280828978167022 11 0.005
+0.0002425402745256439 12 0.01
+0.0001207937198159064 13 0.015

reg_tests/test_files/VOFInertialDroplet/VOFInertialDroplet.yaml

@@ -0,0 +1,136 @@
+Simulations:
+  - name: sim1
+    time_integrator: ti_1
+    optimizer: opt1
+
+linear_solvers:
+
+  - name: solve_scalar
+    type: tpetra
+    method: gmres
+    preconditioner: sgs
+    tolerance: 1e-12
+    max_iterations: 50
+    kspace: 50
+    output_level: 0
+
+  - name: solve_cont
+    type: tpetra
+    method: gmres
+    preconditioner: muelu
+    tolerance: 1e-12
+    max_iterations: 75
+    kspace: 75
+    output_level: 0
+    recompute_preconditioner: yes
+    muelu_xml_file_name: ../../xml/vof_resolved.xml
+
+
+realms:
+
+  - name: realm_1
+    mesh: "generated:75x75x75|bbox:-0.5,-0.5,-0.5,0.5,0.5,0.5|sideset:xXyYzZ|show" 
+    use_edges: yes
+    automatic_decomposition_type: rcb
+
+    equation_systems:
+      name: theEqSys
+      max_iterations: 3
+
+      solver_system_specification:
+        volume_of_fluid: solve_scalar
+        velocity: solve_scalar
+        pressure: solve_cont
+
+      systems:
+        - VolumeOfFluid:
+            name: myVOF
+            max_iterations: 1
+            convergence_tolerance: 1e-8
+        - LowMachEOM:
+            name: myLowMach
+            max_iterations: 1
+            convergence_tolerance: 1e-8
+
+    initial_conditions:
+
+      - user_function: ic_1
+        target_name: block_1
+        user_function_name:
+          volume_of_fluid: droplet
+          velocity: droplet
+        user_function_parameters:
+          volume_of_fluid: [0.0, 0.0, 0.0, 0.15, -100, 0, 0.01]
+          velocity: [0.0, 0.0, 0.0, 0.2, 0.2, 0.2, 0.15, 0.01]
+      - constant: ic_2
+        target_name: block_1
+        value:
+          pressure: 0.0
+
+    material_properties:
+      target_name: block_1
+      specifications:
+        - name: density
+          type: volume_of_fluid
+          primary_value: 1000.0
+          secondary_value: 1.0
+
+        - name: viscosity
+          type: volume_of_fluid
+          primary_value: 1.0e-3
+          secondary_value: 1.0e-5
+
+    boundary_conditions:
+      - periodic_boundary_condition: bc_1
+        target_name: [surface_1, surface_2]
+        periodic_user_data:
+          search_tolerance: 0.001
+      - periodic_boundary_condition: bc_2
+        target_name: [surface_3, surface_4]
+        periodic_user_data:
+          search_tolerance: 0.001
+      - periodic_boundary_condition: bc_3
+        target_name: [surface_5, surface_6]
+        periodic_user_data:
+          search_tolerance: 0.001
+
+    solution_options:
+      name: myOptions
+      projected_timescale_type: momentum_diag_inv
+
+      options:
+        - hybrid_factor:
+            velocity: 0.05
+            volume_of_fluid: 0.0
+
+        - limiter:
+            volume_of_fluid: yes
+            pressure: no
+            velocity: no
+
+        - consistent_mass_matrix_png:
+            pressure: no
+
+    output:
+      output_data_base_name: out/inertialDroplet.e
+      output_frequency: 50
+      output_node_set: yes
+      output_variables:
+       - density
+       - volume_of_fluid
+       - velocity
+       - pressure
+       - dvolume_of_fluiddx
+Time_Integrators:
+  - StandardTimeIntegrator:
+      name: ti_1
+      start_time: 0
+      termination_time: 0.015
+      time_step: 0.005
+      time_stepping_type: fixed
+      time_step_count: 10
+      second_order_accuracy: yes
+
+
+      realms:
+        - realm_1

src/LowMachEquationSystem.C

@@ -174,6 +174,8 @@
 #include <user_functions/PerturbedShearLayerAuxFunctions.h>
 #include <user_functions/GaussJetVelocityAuxFunction.h>
 
+#include <user_functions/DropletVelocityAuxFunction.h>
+
 // deprecated
 
 // stk_util
@@ -663,6 +665,13 @@ LowMachEquationSystem::register_initial_condition_fcn(
       theAuxFunc = new SinProfileChannelFlowVelocityAuxFunction(0, nDim);
     } else if (fcnName == "PerturbedShearLayer") {
       theAuxFunc = new PerturbedShearLayerVelocityAuxFunction(0, nDim);
+    } else if (fcnName == "droplet") {
+      std::map<std::string, std::vector<double>>::const_iterator iterParams =
+        theParams.find(dofName);
+      std::vector<double> fcnParams = (iterParams != theParams.end())
+                                        ? (*iterParams).second
+                                        : std::vector<double>();
+      theAuxFunc = new DropletVelocityAuxFunction(0, nDim, fcnParams);
     } else {
       throw std::runtime_error(
         "InitialCondFunction::non-supported velocity IC");
@@ -2895,6 +2904,12 @@ ContinuityEquationSystem::register_edge_fields(
   massFlowRate_ = &(meta_data.declare_field<double>(
     stk::topology::EDGE_RANK, "mass_flow_rate"));
   stk::mesh::put_field_on_mesh(*massFlowRate_, selector, nullptr);
+
+  if (realm_.solutionOptions_->realm_has_vof_) {
+    auto massVofBalancedFlowRate_ = &(meta_data.declare_field<double>(
+      stk::topology::EDGE_RANK, "mass_vof_balanced_flow_rate"));
+    stk::mesh::put_field_on_mesh(*massVofBalancedFlowRate_, selector, nullptr);
+  }
 }
 
 //--------------------------------------------------------------------------