Add first-generation smoothing/sharpening (#531)

* pde-based compression/sharpening,

    d/dt(alpha) + u_j d/dxj(alpha) + d/dxj(s_j) = 0.0

    with s_j = ucj*alpha(1-alpha)

  activate via:

    volume_of_fluid: [vof, sucv_nso, sharpen]

* explicit smoothing post-solve; dial in through input file.

        - VolumeOfFluid:
            name: myV
            max_iterations: 1
            convergence_tolerance: 1.e-2
            fourier_number: 0.25
            compression_constant: 0.05
            activate_smoothing: yes
            smoothing_iterations: 5

    alphaSmoothed := alpha + Fo*dxMin*dxMin*d^2/dx^2(alpha)

* clarify explicit sharpening output

* add proper updates for overset; assembledRHS and interface normal for SCS
  smoothing approach

* modify vof test to include new sharpening terms; diffs due to sharpening..

Notes:

a) sharpening works great for sample driven cavity cases.
b) advection limitation for stability when mesh spacing between
   overset and background are different.
c) Core VOF transport scheme in... Time to modify the PPE

include/NaluParsing.h

@@ -189,6 +189,7 @@ struct WallUserData : public UserData {
   TurbKinEnergy tke_;
   MixtureFraction mixFrac_;
   MassFraction massFraction_;
+  VolumeOfFluid vof_;
   Emissivity emissivity_;
   Irradiation irradiation_;
   Transmissivity transmissivity_;

include/SolutionOptions.h

@@ -137,7 +137,7 @@ class SolutionOptions
   double mdotAlgAccumulation_;
   double mdotAlgInflow_;
   double mdotAlgOpen_;
- 
+
   // turbulence model coeffs
   std::map<TurbulenceModelConstant, double> turbModelConstantMap_;
 

include/VolumeOfFluidEquationSystem.h

@@ -34,7 +34,11 @@ class VolumeOfFluidEquationSystem : public EquationSystem {
   VolumeOfFluidEquationSystem(
     EquationSystems& equationSystems,
     const bool outputClippingDiag,
-    const double deltaZClip);
+    const double deltaVofClip,
+    const double Fo,
+    const double cAlpha,
+    const bool smooth,
+    const int smoothIter);
   virtual ~VolumeOfFluidEquationSystem();
 
   void populate_derived_quantities();
@@ -86,22 +90,33 @@ class VolumeOfFluidEquationSystem : public EquationSystem {
 
   void sharpen_interface_explicit();
   void smooth_vof();
+  void smooth_vof_execute();
   void compute_interface_normal();
 
   const bool managePNG_;
-  const bool outputClippingDiag_;
-  const double deltaVofClip_;
 
   ScalarFieldType *vof_;
   ScalarFieldType *vofSmoothed_;
-  ScalarFieldType *interfaceNormal_;
+  ScalarFieldType *smoothedRhs_;
+  VectorFieldType *interfaceNormal_;
   VectorFieldType *dvofdx_;
   ScalarFieldType *vofTmp_;
 
   AssembleNodalGradAlgorithmDriver *assembleNodalGradAlgDriver_;
 
   ProjectedNodalGradientEquationSystem *projectedNodalGradEqs_;
 
+  // clipping
+  const bool outputClippingDiag_;
+  const double deltaVofClip_;
+
+  // smoothing and sharepening params
+  const double Fo_;
+  const double cAlpha_;
+  double dxMin_;
+  const bool smooth_;
+  const int smoothIter_;
+
   bool isInit_;
 };
 

include/kernel/VolumeOfFluidSharpenElemKernel.h

@@ -0,0 +1,65 @@
+/*------------------------------------------------------------------------*/
+/*  Copyright 2014 Sandia Corporation.                                    */
+/*  This software is released under the license detailed                  */
+/*  in the file, LICENSE, which is located in the top-level Nalu          */
+/*  directory structure                                                   */
+/*------------------------------------------------------------------------*/
+
+#ifndef VolumeOfFluidSharpenElemKernel_H
+#define VolumeOfFluidSharpenElemKernel_H
+
+#include "kernel/Kernel.h"
+#include "FieldTypeDef.h"
+
+#include <stk_mesh/base/BulkData.hpp>
+#include <stk_mesh/base/Entity.hpp>
+
+#include <Kokkos_Core.hpp>
+
+namespace sierra {
+namespace nalu {
+
+class TimeIntegrator;
+class SolutionOptions;
+class MasterElement;
+class ElemDataRequests;
+
+template<typename AlgTraits>
+class VolumeOfFluidSharpenElemKernel: public Kernel
+{
+public:
+  VolumeOfFluidSharpenElemKernel(
+    const stk::mesh::BulkData&,
+    const SolutionOptions&,
+    ScalarFieldType*,
+    const double,
+    ElemDataRequests&);
+
+  virtual ~VolumeOfFluidSharpenElemKernel();
+
+  /** Execute the kernel within a Kokkos loop and populate the LHS and RHS for
+   *  the linear solve
+   */
+  virtual void execute(
+    SharedMemView<DoubleType**>&,
+    SharedMemView<DoubleType*>&,
+    ScratchViews<DoubleType>&);
+
+private:
+  VolumeOfFluidSharpenElemKernel() = delete;
+
+  VectorFieldType *coordinates_{nullptr};
+  VectorFieldType *velocityRTM_{nullptr};
+  VectorFieldType *interfaceNormal_{nullptr};
+  ScalarFieldType *vofNp1_{nullptr};
+
+  /// Integration point to node mapping
+  const double cAlpha_;
+  const int* ipNodeMap_;
+
+};
+
+}  // nalu
+}  // sierra
+
+#endif /* VolumeOfFluidSharpenElemKernel_H */

include/user_functions/RayleighTaylorMixFracAuxFunction.h

@@ -39,6 +39,7 @@ class RayleighTaylorMixFracAuxFunction : public AuxFunction
   double tX_;
   double yTr_;
   double dTr_;
+  double surf_;
   const double pi_;
 };
 

reg_tests/test_files/vofSphereDrop/vofSphereDrop.i

@@ -51,6 +51,9 @@ realms:
             name: myV
             max_iterations: 1
             convergence_tolerance: 1.e-2
+            activate_smoothing: yes
+            fourier_number: 0.25
+            compression_constant: 0.05
 
     initial_conditions:
       - constant: ic_1
@@ -149,7 +152,7 @@ realms:
         - element_source_terms:
             momentum: [lumped_momentum_time_derivative, advection_diffusion, buoyancy]
             continuity: [advection]
-            volume_of_fluid: [vof, sucv_nso]
+            volume_of_fluid: [vof, sucv_nso, sharpen]
 
         - user_constants:
             gravity: [0.0, -10.0]
@@ -169,6 +172,10 @@ realms:
        - intersected_element
        - mesh_velocity
        - density
+       - volume_of_fluid
+       - volume_of_fluid_smoothed
+       - interface_normal
+       - dvofdx
 
 Time_Integrators:
   - StandardTimeIntegrator:

reg_tests/test_files/vofSphereDrop/vofSphereDrop.norm.gold

@@ -1,20 +1,20 @@
-0.02432116329066863 1 0.005
-0.001966663692513982 2 0.01
-0.0003478067418379743 3 0.015
-0.000128807125321143 4 0.02
-0.0001186665578123386 5 0.025
-6.159589527807734e-05 6 0.03
-2.85084248670952e-05 7 0.035
-1.804536444197644e-05 8 0.04
-1.696814237987382e-05 9 0.045
-1.82415952817005e-05 10 0.05
-2.013158697076528e-05 11 0.055
-2.227717639712913e-05 12 0.06
-2.459260313863444e-05 13 0.065
-2.704331805762574e-05 14 0.07
-2.961771509211273e-05 15 0.075
-3.230282985030123e-05 16 0.08
-3.50955110186631e-05 17 0.085
-3.796971046972843e-05 18 0.09
-4.093875506528796e-05 19 0.095
-4.399335630296404e-05 20 0.1
+0.01447384532257732 1 0.005
+0.01216518334441016 2 0.01
+0.01073768042783061 3 0.015
+0.01051559281054295 4 0.02
+0.01048860402882239 5 0.025
+0.01041243777094482 6 0.03
+0.01035905603141019 7 0.035
+0.01032539481913046 8 0.04
+0.01029914072916844 9 0.045
+0.0102748391713822 10 0.05
+0.01025087973701757 11 0.055
+0.01022698076916539 12 0.06
+0.01020317723200227 13 0.065
+0.01017951831366472 14 0.07
+0.01015602635420198 15 0.075
+0.01013270865108439 16 0.08
+0.01010957223961167 17 0.085
+0.01008660623117398 18 0.09
+0.01006382646893597 19 0.095
+0.01004123400529326 20 0.1

src/EquationSystems.C

@@ -104,7 +104,17 @@ void EquationSystems::load(const YAML::Node & y_node)
           get_if_present_no_default(y_eqsys, "output_clipping_diagnostic", outputClipDiag);
           double deltaVofClip = 0.0;
           get_if_present_no_default(y_eqsys, "clipping_delta", deltaVofClip);
-          eqSys = new VolumeOfFluidEquationSystem(*this, outputClipDiag, deltaVofClip);
+          // vof smoothing/sharpening; provide defaults...
+          double fourierNumber = 0.25;
+          double cAlpha = 0.05;
+          bool smooth = true;
+          int smoothIter = 5;
+          get_if_present_no_default(y_eqsys, "fourier_number", fourierNumber);
+          get_if_present_no_default(y_eqsys, "compression_constant", cAlpha);
+          get_if_present_no_default(y_eqsys, "activate_smoothing", smooth);
+          get_if_present_no_default(y_eqsys, "smoothing_iterations", smoothIter);
+          eqSys = new VolumeOfFluidEquationSystem(*this, outputClipDiag, deltaVofClip, 
+                                                  fourierNumber, cAlpha, smooth, smoothIter);
         }
         else if ( expect_map(y_system, "LowMachEOM", true) ) {
 	  y_eqsys =  expect_map(y_system, "LowMachEOM", true);

src/NaluParsing.C

@@ -821,6 +821,13 @@ namespace YAML
       wallData.bcDataSpecifiedMap_["mass_fraction"] = true;
       wallData.bcDataTypeMap_["mass_fraction"] = sierra::nalu::CONSTANT_UD;
     }
+    if (node["volume_of_fluid"])
+    {
+      wallData.vof_ = node["volume_of_fluid"].as<
+        sierra::nalu::VolumeOfFluid>();
+      wallData.bcDataSpecifiedMap_["volume_of_fluid"] = true;
+      wallData.bcDataTypeMap_["volume_of_fluid"] = sierra::nalu::CONSTANT_UD;
+    }
     if (node["emissivity"])
     {
       wallData.emissivity_ = node["emissivity"].as<sierra::nalu::Emissivity>();