Fix user functions and general cleaning

include/matrix_free/PolynomialOrders.h

@@ -56,27 +56,27 @@ enum {
 #define P_INVOKEABLE(func)                                                     \
   template <int p, typename... Args>                                           \
   auto func(Args&&... args)                                                    \
-    -> decltype(IMPLNAME(func) < p > ::invoke(std::forward<Args>(args)...))    \
+    ->decltype(IMPLNAME(func) < p > ::invoke(std::forward<Args>(args)...))     \
   {                                                                            \
     return IMPLNAME(func)<p>::invoke(std::forward<Args>(args)...);             \
   }
 
 // can't return a value dependent on template parameter
-#define SWITCH_INVOKEABLE(func)                                                    \
-  template <typename... Args>                                                      \
-  auto func(int p, Args&&... args)                                                 \
-    -> decltype(IMPLNAME(func) < inst::P1 > ::invoke(std::forward<Args>(args)...)) \
-  {                                                                                \
-    switch (p) {                                                                   \
-    case inst::P2:                                                                 \
-      return IMPLNAME(func)<inst::P2>::invoke(std::forward<Args>(args)...);        \
-    case inst::P3:                                                                 \
-      return IMPLNAME(func)<inst::P3>::invoke(std::forward<Args>(args)...);        \
-    case inst::P4:                                                                 \
-      return IMPLNAME(func)<inst::P4>::invoke(std::forward<Args>(args)...);        \
-    default:                                                                       \
-      return IMPLNAME(func)<inst::P1>::invoke(std::forward<Args>(args)...);        \
-    }                                                                              \
+#define SWITCH_INVOKEABLE(func)                                                   \
+  template <typename... Args>                                                     \
+  auto func(int p, Args&&... args)                                                \
+    ->decltype(IMPLNAME(func) < inst::P1 > ::invoke(std::forward<Args>(args)...)) \
+  {                                                                               \
+    switch (p) {                                                                  \
+    case inst::P2:                                                                \
+      return IMPLNAME(func)<inst::P2>::invoke(std::forward<Args>(args)...);       \
+    case inst::P3:                                                                \
+      return IMPLNAME(func)<inst::P3>::invoke(std::forward<Args>(args)...);       \
+    case inst::P4:                                                                \
+      return IMPLNAME(func)<inst::P4>::invoke(std::forward<Args>(args)...);       \
+    default:                                                                      \
+      return IMPLNAME(func)<inst::P1>::invoke(std::forward<Args>(args)...);       \
+    }                                                                             \
   }
 
 } // namespace matrix_free

src/LowMachEquationSystem.C

@@ -165,9 +165,6 @@
 #include <user_functions/KovasznayVelocityAuxFunction.h>
 #include <user_functions/KovasznayPressureAuxFunction.h>
 
-#include <user_functions/DropletVelocityAuxFunction.h>
-#include <user_functions/FlatDensityAuxFunction.h>
-
 #include <overset/UpdateOversetFringeAlgorithmDriver.h>
 #include <overset/AssembleOversetPressureAlgorithm.h>
 
@@ -177,8 +174,8 @@
 #include <user_functions/GaussJetVelocityAuxFunction.h>
 
 #include <user_functions/DropletVelocityAuxFunction.h>
-
-#include <user_functions/SloshingTankPressureAuxFunction.h>
+#include <user_functions/DropletVelocityAuxFunction.h>
+#include <user_functions/FlatDensityAuxFunction.h>
 #include <user_functions/WaterLevelDensityAuxFunction.h>
 
 // deprecated
@@ -3519,13 +3516,6 @@ ContinuityEquationSystem::register_initial_condition_fcn(
       theAuxFunc = new TaylorGreenPressureAuxFunction();
     } else if (fcnName == "kovasznay") {
       theAuxFunc = new KovasznayPressureAuxFunction();
-    } else if (fcnName == "sloshing_tank") {
-      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 SloshingTankPressureAuxFunction(fcnParams);
     } else {
       throw std::runtime_error("ContinuityEquationSystem::register_initial_"
                                "condition_fcn: limited functions supported");

src/VolumeOfFluidEquationSystem.C

@@ -41,15 +41,16 @@
 #include "ngp_algorithms/NodalGradEdgeAlg.h"
 #include "ngp_algorithms/NodalGradBndryElemAlg.h"
 #include "stk_topology/topology.hpp"
-#include "user_functions/ZalesakDiskVOFAuxFunction.h"
-#include "user_functions/ZalesakSphereVOFAuxFunction.h"
-#include "user_functions/DropletVOFAuxFunction.h"
-#include "user_functions/SloshingTankVOFAuxFunction.h"
 #include "ngp_utils/NgpFieldBLAS.h"
 #include "ngp_utils/NgpLoopUtils.h"
 #include "ngp_utils/NgpFieldUtils.h"
 #include "stk_io/IossBridge.hpp"
 
+#include "user_functions/ZalesakDiskVOFAuxFunction.h"
+#include "user_functions/ZalesakSphereVOFAuxFunction.h"
+#include "user_functions/DropletVOFAuxFunction.h"
+#include "user_functions/SloshingTankVOFAuxFunction.h"
+
 namespace sierra {
 namespace nalu {
 
@@ -118,15 +119,13 @@ VolumeOfFluidEquationSystem::register_nodal_fields(
   const int nDim = meta_data.spatial_dimension();
   stk::mesh::Selector selector = stk::mesh::selectUnion(part_vec);
 
-  // register dof; set it as a restart variable
   const int numStates = realm_.number_of_states();
 
   auto density_ = &(meta_data.declare_field<double>(
     stk::topology::NODE_RANK, "density", numStates));
   stk::mesh::put_field_on_mesh(*density_, selector, nullptr);
   realm_.augment_restart_variable_list("density");
 
-  // push to property list
   realm_.augment_property_map(DENSITY_ID, density_);
 
   volumeOfFluid_ = &(meta_data.declare_field<double>(
@@ -145,7 +144,6 @@ VolumeOfFluidEquationSystem::register_nodal_fields(
   stk::io::set_field_output_type(
     *dvolumeOfFluiddx_, stk::io::FieldOutputType::VECTOR_3D);
 
-  // delta solution for linear solver; share delta since this is a split system
   vofTmp_ =
     &(meta_data.declare_field<double>(stk::topology::NODE_RANK, "vofTmp"));
   stk::mesh::put_field_on_mesh(*vofTmp_, selector, nullptr);
@@ -273,20 +271,17 @@ VolumeOfFluidEquationSystem::register_inflow_bc(
   const stk::topology& /* partTopo */,
   const InflowBoundaryConditionData& inflowBCData)
 {
-  // algorithm type
   const AlgorithmType algType = INFLOW;
   ScalarFieldType& vofNp1 = volumeOfFluid_->field_of_state(stk::mesh::StateNP1);
   VectorFieldType& dvofdxNone =
     dvolumeOfFluiddx_->field_of_state(stk::mesh::StateNone);
 
   stk::mesh::MetaData& meta_data = realm_.meta_data();
 
-  // register boundary data; vof_bc
   ScalarFieldType* theBcField =
     &(meta_data.declare_field<double>(stk::topology::NODE_RANK, "vof_bc"));
   stk::mesh::put_field_on_mesh(*theBcField, *part, nullptr);
 
-  // extract the value for user specified vof and save off the AuxFunction
   InflowUserData userData = inflowBCData.userData_;
   std::string vofName = "volume_of_fluid";
   UserDataType theDataType = get_bc_data_type(userData, vofName);
@@ -307,29 +302,22 @@ VolumeOfFluidEquationSystem::register_inflow_bc(
                              "only constant functions supported");
   }
 
-  // bc data alg
   AuxFunctionAlgorithm* auxAlg = new AuxFunctionAlgorithm(
     realm_, part, theBcField, theAuxFunc, stk::topology::NODE_RANK);
 
-  // how to populate the field?
   if (userData.externalData_) {
-    // xfer will handle population; only need to populate the initial value
     realm_.initCondAlg_.push_back(auxAlg);
   } else {
-    // put it on bcData
     bcDataAlg_.push_back(auxAlg);
   }
 
-  // copy vof_bc to vof_transition np1...
   CopyFieldAlgorithm* theCopyAlg = new CopyFieldAlgorithm(
     realm_, part, theBcField, &vofNp1, 0, 1, stk::topology::NODE_RANK);
   bcDataMapAlg_.push_back(theCopyAlg);
 
-  // non-solver; dvofdx; allow for element-based shifted
   nodalGradAlgDriver_.register_face_algorithm<ScalarNodalGradBndryElemAlg>(
     algType, part, "vof_nodal_grad", &vofNp1, &dvofdxNone, edgeNodalGradient_);
 
-  // Dirichlet bc
   std::map<AlgorithmType, SolverAlgorithm*>::iterator itd =
     solverAlgDriver_->solverDirichAlgMap_.find(algType);
   if (itd == solverAlgDriver_->solverDirichAlgMap_.end()) {
@@ -370,7 +358,6 @@ VolumeOfFluidEquationSystem::register_wall_bc(
   const stk::topology& /*theTopo*/,
   const WallBoundaryConditionData& /*wallBCData*/)
 {
-  // algorithm type
   const AlgorithmType algType = SYMMETRY;
 
   ScalarFieldType& vofNp1 = volumeOfFluid_->field_of_state(stk::mesh::StateNP1);
@@ -391,7 +378,6 @@ VolumeOfFluidEquationSystem::register_symmetry_bc(
   const stk::topology& /*theTopo*/,
   const SymmetryBoundaryConditionData& /* symmetryBCData */)
 {
-  // algorithm type
   const AlgorithmType algType = SYMMETRY;
 
   ScalarFieldType& vofNp1 = volumeOfFluid_->field_of_state(stk::mesh::StateNP1);
@@ -499,10 +485,8 @@ VolumeOfFluidEquationSystem::register_initial_condition_fcn(
           stk::topology::NODE_RANK, "velocity");
         ScalarFieldType* density_ = realm_.meta_data().get_field<double>(
           stk::topology::NODE_RANK, "density");
-        std::vector<double> userSpec(3);
-        userSpec[0] = 1.0;
-        userSpec[1] = 0.0;
-        userSpec[2] = 0.0;
+
+        std::vector<double> userSpec = {1.0, 0.0, 0.0};
         AuxFunction* constantAuxFunc = new ConstantAuxFunction(0, 1, userSpec);
         AuxFunctionAlgorithm* constantAuxAlg = new AuxFunctionAlgorithm(
           realm_, part, density_, constantAuxFunc, stk::topology::NODE_RANK);

src/aero/fsi/FSIturbine.C

@@ -653,7 +653,7 @@ fsiTurbine::write_nc_def_loads(const size_t tStep, const double curTime)
         (0.5 * (brFSIdata_.bld_rloc[i + iStart + 1] -
                 brFSIdata_.bld_rloc[i + iStart - 1]));
     }
-    brFSIdata_.bld_ld[(iStart) * 6 + 4] =
+    brFSIdata_.bld_ld[(iStart)*6 + 4] =
       (0.5 * (brFSIdata_.bld_rloc[iStart + 1] - brFSIdata_.bld_rloc[iStart]));
     brFSIdata_.bld_ld[(iStart + nBldPts - 1) * 6 + 4] =
       (0.5 * (brFSIdata_.bld_rloc[iStart + nBldPts - 1] -

src/aero/fsi/MapLoad.C

@@ -110,11 +110,11 @@ mapTowerLoad(
 
         // Find the interpolated reference position first
         linInterpVec(
-          &twrRefPos[(loadMap_bip) * 6], &twrRefPos[(loadMap_bip + 1) * 6],
+          &twrRefPos[(loadMap_bip)*6], &twrRefPos[(loadMap_bip + 1) * 6],
           loadMapInterp_bip, tmpNodePos.data());
         // Find the interpolated linear displacement
         linInterpVec(
-          &twrDef[(loadMap_bip) * 6], &twrDef[(loadMap_bip + 1) * 6],
+          &twrDef[(loadMap_bip)*6], &twrDef[(loadMap_bip + 1) * 6],
           loadMapInterp_bip, tmpNodeDisp.data());
         // Add displacement to find actual position
         for (auto idim = 0; idim < 3; idim++)
@@ -130,8 +130,8 @@ mapTowerLoad(
         // Split the force and moment into the two surrounding nodes in a
         // variationally consistent manner using the interpolation factor
         fsi::splitForceMoment(
-          tmpForceMoment.data(), loadMapInterp_bip,
-          &(twrLoad[(loadMap_bip) * 6]), &(twrLoad[(loadMap_bip + 1) * 6]));
+          tmpForceMoment.data(), loadMapInterp_bip, &(twrLoad[(loadMap_bip)*6]),
+          &(twrLoad[(loadMap_bip + 1) * 6]));
       }
     }
   }

src/aero/fsi/OpenfastFSI.C

@@ -450,8 +450,8 @@ OpenfastFSI::send_loads(const double curTime)
           (fsiTurbineData_[i]->params_.nBRfsiPtsTwr) * 6, MPI_DOUBLE, MPI_SUM,
           turbProc, bulk_->parallel());
         iError = MPI_Reduce(
-          fsiTurbineData_[i]->brFSIdata_.bld_ld.data(), NULL,
-          (nTotBldNodes) * 6, MPI_DOUBLE, MPI_SUM, turbProc, bulk_->parallel());
+          fsiTurbineData_[i]->brFSIdata_.bld_ld.data(), NULL, (nTotBldNodes)*6,
+          MPI_DOUBLE, MPI_SUM, turbProc, bulk_->parallel());
       }
     }
   }
@@ -664,8 +664,8 @@ OpenfastFSI::map_loads(const int tStep, const double curTime)
           (fsiTurbineData_[i]->params_.nBRfsiPtsTwr) * 6, MPI_DOUBLE, MPI_SUM,
           turbProc, bulk_->parallel());
         iError = MPI_Reduce(
-          fsiTurbineData_[i]->brFSIdata_.bld_ld.data(), NULL,
-          (nTotBldNodes) * 6, MPI_DOUBLE, MPI_SUM, turbProc, bulk_->parallel());
+          fsiTurbineData_[i]->brFSIdata_.bld_ld.data(), NULL, (nTotBldNodes)*6,
+          MPI_DOUBLE, MPI_SUM, turbProc, bulk_->parallel());
       }
 
       fsiTurbineData_[i]->write_nc_def_loads(tStep, curTime);

src/edge_kernels/MomentumEdgeSolverAlg.C

@@ -68,14 +68,14 @@ MomentumEdgeSolverAlg::execute()
   const std::string dofName = "velocity";
   const DblType includeDivU = realm_.get_divU();
   const DblType alpha = realm_.get_alpha_factor(dofName);
-  const DblType alphaUpw = realm_.get_alpha_upw_factor(dofName);
+  const DblType alphaUpw_input = realm_.get_alpha_upw_factor(dofName);
   const DblType hoUpwind = realm_.get_upw_factor(dofName);
   const DblType relaxFacU =
     realm_.solutionOptions_->get_relaxation_factor(dofName);
   const bool useLimiter = realm_.primitive_uses_limiter(dofName);
 
   const DblType om_alpha = 1.0 - alpha;
-  const DblType om_alphaUpw = 1.0 - alphaUpw;
+  const DblType om_alphaUpw_input = 1.0 - alphaUpw_input;
 
   const DblType has_vof = (double)realm_.solutionOptions_->realm_has_vof_;
 
@@ -100,6 +100,12 @@ MomentumEdgeSolverAlg::execute()
       ShmemDataType & smdata, const stk::mesh::FastMeshIndex& edge,
       const stk::mesh::FastMeshIndex& nodeL,
       const stk::mesh::FastMeshIndex& nodeR) {
+      // Variables are read-only when C++ captures by value into a lambda
+      // This avoids this issue, allowing modification of upwinding parameters
+      // for VOF
+      DblType alphaUpw = alphaUpw_input;
+      DblType om_alphaUpw = om_alphaUpw_input;
+
       // Scratch work array for edgeAreaVector
       NALU_ALIGNED DblType av[NDimMax_];
       // Populate area vector work array
@@ -165,18 +171,16 @@ MomentumEdgeSolverAlg::execute()
       // at interfaces with interface
       // widths that are ~2 cells thick
       DblType density_upwinding_factor = 1.0;
-      DblType alphaUpw_w_vof = alphaUpw;
-      DblType om_alphaUpw_w_vof = 1.0 - alphaUpw_w_vof;
       if (has_vof > 0.5) {
         const DblType min_density = stk::math::min(densityL, densityR);
         const DblType density_differential =
           stk::math::abs(densityL - densityR) / min_density;
         density_upwinding_factor =
           1.0 - stk::math::erf(6.0 * density_differential);
 
-        alphaUpw_w_vof = density_upwinding_factor * alphaUpw_w_vof +
-                         (1.0 - density_upwinding_factor);
-        om_alphaUpw_w_vof = 1.0 - alphaUpw_w_vof;
+        alphaUpw = density_upwinding_factor * alphaUpw +
+                   (1.0 - density_upwinding_factor);
+        om_alphaUpw = 1.0 - alphaUpw;
         pecfac =
           1.0 - density_upwinding_factor + density_upwinding_factor * pecfac;
         om_pecfac = 1.0 - pecfac;
@@ -234,9 +238,9 @@ MomentumEdgeSolverAlg::execute()
         const DblType uiIp = 0.5 * (vel.get(nodeR, i) + vel.get(nodeL, i));
 
         // Upwind contribution
-        const DblType uiUpw =
-          (mdot > 0.0) ? (alphaUpw_w_vof * uIpL[i] + om_alphaUpw_w_vof * uiIp)
-                       : (alphaUpw_w_vof * uIpR[i] + om_alphaUpw_w_vof * uiIp);
+        const DblType uiUpw = (mdot > 0.0)
+                                ? (alphaUpw * uIpL[i] + om_alphaUpw * uiIp)
+                                : (alphaUpw * uIpR[i] + om_alphaUpw * uiIp);
 
         const DblType uiHatL = (alpha * uIpL[i] + om_alpha * uiIp);
         const DblType uiHatR = (alpha * uIpR[i] + om_alpha * uiIp);

src/user_functions/CMakeLists.txt

@@ -45,6 +45,5 @@ target_sources(nalu PRIVATE
    ${CMAKE_CURRENT_SOURCE_DIR}/DropletVelocityAuxFunction.C
    ${CMAKE_CURRENT_SOURCE_DIR}/FlatDensityAuxFunction.C
    ${CMAKE_CURRENT_SOURCE_DIR}/SloshingTankVOFAuxFunction.C
-   ${CMAKE_CURRENT_SOURCE_DIR}/SloshingTankPressureAuxFunction.C
    ${CMAKE_CURRENT_SOURCE_DIR}/WaterLevelDensityAuxFunction.C
 )

src/user_functions/FlatDensityAuxFunction.C

@@ -42,11 +42,11 @@ FlatDensityAuxFunction::do_evaluate(
     const double z = coords[2];
     const double interface_thickness = 0.015;
 
-    fieldPtr[0] = 0.0;
-    fieldPtr[0] += -0.5 * (std::erf(y / interface_thickness) + 1.0) + 1.0;
+    const double vof_function =
+      -0.5 * (std::erf(y / interface_thickness) + 1.0) + 1.0;
 
     // air-water
-    fieldPtr[0] = 1000.0 * fieldPtr[0] + (1.0 - fieldPtr[0]);
+    fieldPtr[0] = 1000.0 * vof_function + (1.0 - vof_function);
 
     fieldPtr += fieldSize;
     coords += spatialDimension;