Zalesak Sphere case - designed to match AMR-Wind

include/user_functions/ZalesakSphereMassFlowRateKernel.h

@@ -0,0 +1,38 @@
+// 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 ZALESAKSPHEREMASSFLOWRATEALG_H
+#define ZALESAKSPHEREMASSFLOWRATEALG_H
+
+#include "AssembleEdgeSolverAlgorithm.h"
+
+namespace sierra {
+namespace nalu {
+
+class ZalesakSphereMassFlowRateEdgeAlg : public AssembleEdgeSolverAlgorithm
+{
+public:
+  // TODO: refactor to use FieldManager
+  ZalesakSphereMassFlowRateEdgeAlg(
+    Realm&, stk::mesh::Part*, EquationSystem*, const bool = false);
+
+  virtual ~ZalesakSphereMassFlowRateEdgeAlg() = default;
+
+  virtual void execute();
+
+private:
+  unsigned coordinates_{stk::mesh::InvalidOrdinal};
+  unsigned edgeAreaVec_{stk::mesh::InvalidOrdinal};
+  unsigned massFlowRate_{stk::mesh::InvalidOrdinal};
+};
+
+} // namespace nalu
+} // namespace sierra
+
+#endif /* ZALESAKSPHEREMASSFLOWRATEALG_H */

include/user_functions/ZalesakSphereVOFAuxFunction.h

@@ -0,0 +1,42 @@
+// 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 ZalesakSphereVOFAuxFunction_h
+#define ZalesakSphereVOFAuxFunction_h
+
+#include <AuxFunction.h>
+
+#include <vector>
+
+namespace sierra {
+namespace nalu {
+
+class ZalesakSphereVOFAuxFunction : public AuxFunction
+{
+public:
+  ZalesakSphereVOFAuxFunction();
+
+  virtual ~ZalesakSphereVOFAuxFunction() {}
+
+  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;
+};
+
+} // namespace nalu
+} // namespace sierra
+
+#endif

src/VolumeOfFluidEquationSystem.C

@@ -29,6 +29,7 @@
 // edge kernels
 #include "edge_kernels/VOFAdvectionEdgeAlg.h"
 #include "user_functions/ZalesakDiskMassFlowRateKernel.h"
+#include "user_functions/ZalesakSphereMassFlowRateKernel.h"
 
 // node kernels
 #include "node_kernels/NodeKernelUtils.h"
@@ -40,6 +41,7 @@
 #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 "ngp_utils/NgpFieldBLAS.h"
 #include "ngp_utils/NgpLoopUtils.h"
@@ -491,6 +493,27 @@ VolumeOfFluidEquationSystem::register_initial_condition_fcn(
           new ZalesakDiskMassFlowRateEdgeAlg(realm_, part, this, useAvgMdot);
         realm_.initCondAlg_.push_back(VOFSetMassFlowRate);
       }
+    } else if (fcnName == "zalesak_sphere") {
+      theAuxFunc = new ZalesakSphereVOFAuxFunction();
+      // Initialize mass flow rate until momentum connection implemented
+      {
+        const bool useAvgMdot = (realm_.solutionOptions_->turbulenceModel_ ==
+                                 TurbulenceModel::SST_AMS)
+                                  ? true
+                                  : false;
+        ScalarFieldType* density_ =
+          realm_.meta_data().get_field<ScalarFieldType>(
+            stk::topology::NODE_RANK, "density");
+        std::vector<double> userSpec(1);
+        userSpec[0] = 1.0;
+        AuxFunction* constantAuxFunc = new ConstantAuxFunction(0, 1, userSpec);
+        AuxFunctionAlgorithm* constantAuxAlg = new AuxFunctionAlgorithm(
+          realm_, part, density_, constantAuxFunc, stk::topology::NODE_RANK);
+        realm_.initCondAlg_.push_back(constantAuxAlg);
+        auto VOFSetMassFlowRate =
+          new ZalesakSphereMassFlowRateEdgeAlg(realm_, part, this, useAvgMdot);
+        realm_.initCondAlg_.push_back(VOFSetMassFlowRate);
+      }
     } else if (fcnName == "droplet") {
       theAuxFunc = new DropletVOFAuxFunction();
     } else {

src/user_functions/CMakeLists.txt

@@ -37,5 +37,7 @@ target_sources(nalu PRIVATE
    ${CMAKE_CURRENT_SOURCE_DIR}/WindEnergyTaylorVortexPressureAuxFunction.C
    ${CMAKE_CURRENT_SOURCE_DIR}/ZalesakDiskVOFAuxFunction.C
    ${CMAKE_CURRENT_SOURCE_DIR}/ZalesakDiskMassFlowRateKernel.C
+   ${CMAKE_CURRENT_SOURCE_DIR}/ZalesakSphereVOFAuxFunction.C
+   ${CMAKE_CURRENT_SOURCE_DIR}/ZalesakSphereMassFlowRateKernel.C
    ${CMAKE_CURRENT_SOURCE_DIR}/DropletVOFAuxFunction.C
 )

src/user_functions/ZalesakSphereMassFlowRateKernel.C

@@ -0,0 +1,79 @@
+// 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.
+//
+
+#include "user_functions/ZalesakSphereMassFlowRateKernel.h"
+#include "EquationSystem.h"
+#include "PecletFunction.h"
+#include "SolutionOptions.h"
+#include "utils/StkHelpers.h"
+#include "edge_kernels/EdgeKernelUtils.h"
+#include "stk_mesh/base/NgpField.hpp"
+#include "stk_mesh/base/Types.hpp"
+
+namespace sierra {
+namespace nalu {
+
+ZalesakSphereMassFlowRateEdgeAlg::ZalesakSphereMassFlowRateEdgeAlg(
+  Realm& realm,
+  stk::mesh::Part* part,
+  EquationSystem* eqSystem,
+  const bool useAverages)
+  : AssembleEdgeSolverAlgorithm(realm, part, eqSystem)
+{
+  const auto& meta = realm.meta_data();
+
+  coordinates_ = get_field_ordinal(meta, realm.get_coordinates_name());
+  edgeAreaVec_ =
+    get_field_ordinal(meta, "edge_area_vector", stk::topology::EDGE_RANK);
+  massFlowRate_ = get_field_ordinal(
+    meta, (useAverages) ? "average_mass_flow_rate" : "mass_flow_rate",
+    stk::topology::EDGE_RANK);
+}
+
+void
+ZalesakSphereMassFlowRateEdgeAlg::execute()
+{
+  const double eps = 1.0e-16;
+  const int ndim = realm_.meta_data().spatial_dimension();
+
+  // Defaults from AMR-Wind
+  const double period = 6.0;
+  const double xrot = 0.5;
+  const double yrot = 0.5;
+
+  // STK stk::mesh::NgpField instances for capture by lambda
+  const auto& fieldMgr = realm_.ngp_field_manager();
+  const auto coordinates = fieldMgr.get_field<double>(coordinates_);
+  const auto edgeAreaVec = fieldMgr.get_field<double>(edgeAreaVec_);
+  const auto massFlowRate = fieldMgr.get_field<double>(massFlowRate_);
+
+  run_algorithm(
+    realm_.bulk_data(),
+    KOKKOS_LAMBDA(
+      ShmemDataType & smdata, const stk::mesh::FastMeshIndex& edge,
+      const stk::mesh::FastMeshIndex& nodeL,
+      const stk::mesh::FastMeshIndex& nodeR) {
+      // Scratch work array for edgeAreaVector
+      NALU_ALIGNED DblType av[NDimMax_];
+      // Populate area vector work array
+      for (int d = 0; d < ndim; ++d)
+        av[d] = edgeAreaVec.get(edge, d);
+
+      NALU_ALIGNED DblType edge_centroid[NDimMax_];
+      for (int d = 0; d < ndim; ++d)
+        edge_centroid[d] =
+          0.5 * coordinates.get(nodeL, d) + 0.5 * coordinates.get(nodeR, d);
+
+      massFlowRate.get(edge, 0) = 2.0 * M_PI / period *
+        ((yrot - edge_centroid[1]) * av[0] + (edge_centroid[0] - xrot) * av[1]);
+    });
+}
+
+} // namespace nalu
+} // namespace sierra

src/user_functions/ZalesakSphereVOFAuxFunction.C

@@ -0,0 +1,66 @@
+// 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.
+//
+
+#include <user_functions/ZalesakSphereVOFAuxFunction.h>
+#include <algorithm>
+
+// basic c++
+#include <cmath>
+#include <vector>
+#include <stdexcept>
+
+namespace sierra {
+namespace nalu {
+
+ZalesakSphereVOFAuxFunction::ZalesakSphereVOFAuxFunction() : AuxFunction(0, 1)
+{
+  // does nothing
+}
+
+void
+ZalesakSphereVOFAuxFunction::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
+{
+  for (unsigned p = 0; p < numPoints; ++p) {
+
+    const double x = coords[0];
+    const double y = coords[1];
+    const double z = coords[2];
+
+    // These are the default arguments from the corresponding case in AMR-Wind
+    const double xc = 0.5;
+    const double yc = 0.72;
+    const double zc = 0.24;
+    const double radius = 0.16;
+    const double depth = 0.2;
+    const double width = 0.04;
+
+    fieldPtr[0] = 0.0;
+    // Put VOF in sphere
+    if ((x - xc) * (x - xc) + (y - yc) * (y - yc) + (z - zc) * (z - zc) < radius * radius)
+      fieldPtr[0] = 1.0;
+
+    // Remove slot
+    if (x - xc > -0.5 * width && x - xc < 0.5 * width && y - yc > radius - depth)
+      fieldPtr[0] = 0.0;
+
+    fieldPtr += fieldSize;
+    coords += spatialDimension;
+  }
+}
+
+} // namespace nalu
+} // namespace sierra