add a string-function temperature IC/dirichlet bc

include/EnthalpyEquationSystem.h

@@ -82,10 +82,14 @@ class EnthalpyEquationSystem : public EquationSystem
   void initialize();
   void reinitialize_linear_system();
 
-  virtual void register_initial_condition_fcn(
+  void register_initial_condition_fcn(
+    stk::mesh::Part* /* part */,
+    const std::map<std::string, std::string>& /* theNames */,
+    const std::map<std::string, std::vector<double>>& /* theParams */) final;
+
+  void register_initial_condition_string_function(
     stk::mesh::Part* part,
-    const std::map<std::string, std::string>& theNames,
-    const std::map<std::string, std::vector<double>>& theParams);
+    const std::map<std::string, std::string>& func) final;
 
   void predict_state();
 

include/Enums.h

@@ -111,6 +111,7 @@ enum UserDataType {
   CONSTANT_UD = 0,
   FUNCTION_UD = 1,
   USER_SUB_UD = 2,
+  STRING_FUNCTION_UD = 3,
   UserDataType_END
 };
 

include/EquationSystem.h

@@ -233,6 +233,12 @@ class EquationSystem
   {
   }
 
+  virtual void register_initial_condition_string_function(
+    stk::mesh::Part* /*part*/,
+    const std::map<std::string, std::string>& /*func*/)
+  {
+  }
+
   // rip through the propertyAlg_
   virtual void evaluate_properties();
 

include/EquationSystems.h

@@ -111,6 +111,9 @@ class EquationSystems
   void register_initial_condition_fcn(
     stk::mesh::Part* part, const UserFunctionInitialConditionData& fcnIC);
 
+  void register_initial_condition_string_function(
+    stk::mesh::Part* part, const std::map<std::string, std::string>& func);
+
   void initialize();
   void reinitialize_linear_system();
   void populate_derived_quantities();

include/NaluParsing.h

@@ -37,6 +37,7 @@ struct UserData
   std::map<std::string, std::string> userFunctionMap_;
   std::map<std::string, std::vector<double>> functionParams_;
   std::map<std::string, std::vector<std::string>> functionStringParams_;
+  std::map<std::string, std::string> functions;
 
   // FIXME: must elevate temperature due to the temperature_bc_setup method
   Temperature temperature_;
@@ -112,6 +113,7 @@ struct InflowUserData : public UserData
   bool mixFracSpec_;
   bool massFractionSpec_;
   bool gammaSpec_;
+
   InflowUserData()
     : UserData(),
       uSpec_(false),
@@ -362,6 +364,11 @@ struct UserFunctionInitialConditionData : public InitialCondition
   std::map<std::string, std::vector<double>> functionParams_;
 };
 
+struct StringFunctionInitialConditionData : public InitialCondition
+{
+  std::map<std::string, std::string> functions_;
+};
+
 inline bool
 string_represents_positive_integer(std::string v)
 {
@@ -477,6 +484,8 @@ void operator>>(const YAML::Node& node, NonConformalBoundaryConditionData& rhs);
 void operator>>(const YAML::Node& node, ConstantInitialConditionData& rhs);
 
 void operator>>(const YAML::Node& node, UserFunctionInitialConditionData& rhs);
+void
+operator>>(const YAML::Node& node, StringFunctionInitialConditionData& rhs);
 
 void operator>>(const YAML::Node& node, std::map<std::string, bool>& mapName);
 void operator>>(const YAML::Node& node, std::map<std::string, double>& mapName);

include/master_element/MasterElementFunctions.h

@@ -77,10 +77,10 @@ generic_grad_op(
   static_assert(
     std::is_same<ftype, typename OutputViewType::value_type>::value,
     "Incompatiable value type for views");
-  static_assert(GradViewType::Rank == 3, "grad view assumed to be rank 3");
+  static_assert(GradViewType::rank == 3, "grad view assumed to be rank 3");
   static_assert(
-    CoordViewType::Rank == 2, "Coordinate view assumed to be rank 2");
-  static_assert(OutputViewType::Rank == 3, "Weight view assumed to be rank 3");
+    CoordViewType::rank == 2, "Coordinate view assumed to be rank 2");
+  static_assert(OutputViewType::rank == 3, "Weight view assumed to be rank 3");
 
   ThrowAssert(AlgTraits::nodesPerElement_ == referenceGradWeights.extent(1));
   ThrowAssert(AlgTraits::nDim_ == referenceGradWeights.extent(2));
@@ -148,9 +148,9 @@ generic_gij_3d(
   static_assert(
     std::is_same<ftype, typename OutputViewType::value_type>::value,
     "Incompatiable value type for views");
-  static_assert(GradViewType::Rank == 3, "grad view assumed to be 3D");
-  static_assert(CoordViewType::Rank == 2, "Coordinate view assumed to be 2D");
-  static_assert(OutputViewType::Rank == 3, "gij view assumed to be 3D");
+  static_assert(GradViewType::rank == 3, "grad view assumed to be 3D");
+  static_assert(CoordViewType::rank == 2, "Coordinate view assumed to be 2D");
+  static_assert(OutputViewType::rank == 3, "gij view assumed to be 3D");
   static_assert(AlgTraits::nDim_ == 3, "3D method");
 
   for (unsigned ip = 0; ip < referenceGradWeights.extent(0); ++ip) {
@@ -331,9 +331,9 @@ generic_Mij_2d(
   static_assert(
     std::is_same<ftype, typename OutputViewType::value_type>::value,
     "Incompatiable value type for views");
-  static_assert(GradViewType::Rank == 3, "grad view assumed to be 3D");
-  static_assert(CoordViewType::Rank == 2, "Coordinate view assumed to be 2D");
-  static_assert(OutputViewType::Rank == 3, "Mij view assumed to be 3D");
+  static_assert(GradViewType::rank == 3, "grad view assumed to be 3D");
+  static_assert(CoordViewType::rank == 2, "Coordinate view assumed to be 2D");
+  static_assert(OutputViewType::rank == 3, "Mij view assumed to be 3D");
   static_assert(AlgTraits::nDim_ == 2, "2D method");
 
   const int npe = AlgTraits::nodesPerElement_;
@@ -514,9 +514,9 @@ generic_Mij_3d(
   static_assert(
     std::is_same<ftype, typename OutputViewType::value_type>::value,
     "Incompatiable value type for views");
-  static_assert(GradViewType::Rank == 3, "grad view assumed to be 3D");
-  static_assert(CoordViewType::Rank == 2, "Coordinate view assumed to be 2D");
-  static_assert(OutputViewType::Rank == 3, "Mij view assumed to be 3D");
+  static_assert(GradViewType::rank == 3, "grad view assumed to be 3D");
+  static_assert(CoordViewType::rank == 2, "Coordinate view assumed to be 2D");
+  static_assert(OutputViewType::rank == 3, "Mij view assumed to be 3D");
   static_assert(AlgTraits::nDim_ == 3, "3D method");
 
   for (unsigned ip = 0; ip < referenceGradWeights.extent(0); ++ip) {
@@ -593,9 +593,9 @@ generic_determinant_3d(
   static_assert(
     std::is_same<ftype, typename OutputViewType::value_type>::value,
     "Incompatiable value type for views");
-  static_assert(GradViewType::Rank == 3, "grad view assumed to be 3D");
-  static_assert(CoordViewType::Rank == 2, "Coordinate view assumed to be 2D");
-  static_assert(OutputViewType::Rank == 1, "Weight view assumed to be 1D");
+  static_assert(GradViewType::rank == 3, "grad view assumed to be 3D");
+  static_assert(CoordViewType::rank == 2, "Coordinate view assumed to be 2D");
+  static_assert(OutputViewType::rank == 1, "Weight view assumed to be 1D");
   static_assert(AlgTraits::nDim_ == 3, "3D method");
 
   ThrowAssert(AlgTraits::nodesPerElement_ == referenceGradWeights.extent(1));

include/user_functions/StringTimeCoordFunction.h

@@ -0,0 +1,61 @@
+// 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 StringTimeCoordFunction_h
+#define StringTimeCoordFunction_h
+
+#include "stk_expreval/Evaluator.hpp"
+
+#include "Kokkos_Core.hpp"
+
+#include <string>
+
+namespace sierra::nalu {
+
+namespace fcn {
+inline constexpr int UNMAPPED_INDEX = -1;
+}
+
+class StringTimeCoordFunction
+{
+public:
+  StringTimeCoordFunction(std::string fcn);
+  KOKKOS_FUNCTION double
+  operator()(double t, double x, double y, double z) const;
+  [[nodiscard]] KOKKOS_FUNCTION bool is_constant() const { return constant; }
+
+  KOKKOS_FUNCTION int spatial_dim() const
+  {
+    if (z_index != fcn::UNMAPPED_INDEX) {
+      return 3;
+    } else if (y_index != fcn::UNMAPPED_INDEX) {
+      return 2;
+    } else if (x_index != fcn::UNMAPPED_INDEX) {
+      return 1;
+    }
+    return 0;
+  }
+
+  [[nodiscard]] KOKKOS_FUNCTION bool is_spatial() const
+  {
+    return spatial_dim() != 0;
+  }
+
+private:
+  stk::expreval::ParsedEval<> parsed_eval;
+  bool constant = false;
+  int t_index = fcn::UNMAPPED_INDEX;
+  int x_index = fcn::UNMAPPED_INDEX;
+  int y_index = fcn::UNMAPPED_INDEX;
+  int z_index = fcn::UNMAPPED_INDEX;
+};
+
+} // namespace sierra::nalu
+
+#endif

include/user_functions/StringTimeCoordTemperatureAuxFunction.h

@@ -0,0 +1,41 @@
+// 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 StringTimeCoordTemperatureAuxFunction_h
+#define StringTimeCoordTemperatureAuxFunction_h
+
+#include "AuxFunction.h"
+#include "StringTimeCoordFunction.h"
+
+#include <vector>
+
+namespace sierra::nalu {
+
+class StringTimeCoordTemperatureAuxFunction final : public AuxFunction
+{
+public:
+  StringTimeCoordTemperatureAuxFunction(std::string fcn);
+
+  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 final;
+
+private:
+  const StringTimeCoordFunction f_;
+};
+
+} // namespace sierra::nalu
+
+#endif

include/user_functions/CappingInversionTemperatureAuxFunction.h → include/user_functions/TabulatedTemperatureAuxFunction.h

@@ -14,29 +14,29 @@
 
 #include <vector>
 
-namespace sierra {
-namespace nalu {
+namespace sierra::nalu {
 
-class CappingInversionTemperatureAuxFunction : public AuxFunction
+class TabulatedTemperatureAuxFunction final : public AuxFunction
 {
 public:
-  CappingInversionTemperatureAuxFunction();
+  TabulatedTemperatureAuxFunction(
+    std::vector<double> heights, std::vector<double> temperatures);
 
-  virtual ~CappingInversionTemperatureAuxFunction() {}
-
-  using AuxFunction::do_evaluate;
-  virtual void do_evaluate(
+  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;
+    const unsigned endPos) const final;
+
+private:
+  std::vector<double> heights_;
+  std::vector<double> temperatures_;
 };
 
-} // namespace nalu
-} // namespace sierra
+} // namespace sierra::nalu
 
 #endif

reg_tests/test_files/ablHill3dSymPenalty/ablHill3dSymPenalty.norm.gold

@@ -1,5 +1,5 @@
-1.994475350611083e-05 1 0.02
-1.245775182965687e-05 2 0.04
-4.162417946007383e-06 3 0.06
-1.630724817616641e-07 4 0.08
-1.300056783839082e-08 5 0.1
+0.0001340341715347378 1 0.02
+7.668636308766615e-05 2 0.04
+3.862864631832822e-05 3 0.06
+2.872666841740503e-05 4 0.08
+2.453914784048959e-05 5 0.1

reg_tests/test_files/ablHill3dSymPenalty/ablHill3dSymPenalty.yaml

@@ -50,6 +50,7 @@ realms:
       solver_system_specification:
         velocity: solve_scalar
         pressure: solve_cont
+        enthalpy: solve_scalar
 
       # This defines the equation systems, maximum number of inner iterations,
       # and scaled nonlinear residual tolerance.
@@ -60,6 +61,11 @@ realms:
             max_iterations: 1
             convergence_tolerance: 1.0e-5
 
+        - Enthalpy:
+            name: myEnth
+            max_iterations: 1
+            convergence_tolerance: 1.0e-5
+
     # Specify the properties of the fluid, in this case air.
     material_properties:
 
@@ -96,9 +102,14 @@ realms:
         target_name: [fluid]
         value:
           pressure: 0.0
-#          velocity: [0.707, 0.707, 0.0]
           velocity: [0.500, 0.866, 0.0]
 
+      - string_function: ic_2
+        target_name: [fluid]
+        function:
+          temperature: &tfun '(z <= .650) * 300.0 + (z >= .650 && z < .750) * (300 + .08 * (z*1000-650)) + (z >= .750) * (308 + 0.003 * (z*1000-750))'
+
+
 
     # Boundary conditions are inflow on the west, open (outflow) on the east,
     # inflow on the south, open (outflow) on the north.  The lower boundary
@@ -110,23 +121,27 @@ realms:
       target_name: west
       inflow_user_data:
         velocity: [0.500, 0.866, 0.0]
+        temperature: *tfun
 
     - open_boundary_condition: bc_east
       target_name: east
       open_user_data:
         velocity: [0.500, 0.866, 0.0]
         pressure: 0.0
+        temperature: 300
 
     - inflow_boundary_condition: bc_south
       target_name: south
       inflow_user_data:
         velocity: [0.500, 0.866, 0.0]
+        temperature: *tfun
 
     - open_boundary_condition: bc_north
       target_name: north
       open_user_data:
         velocity: [0.500, 0.866, 0.0]
         pressure: 0.0
+        temperature: 300.
 
     - symmetry_boundary_condition: bc_upper
       target_name: top
@@ -147,6 +162,9 @@ realms:
             velocity: no
             enthalpy: yes 
 
+        - laminar_prandtl:
+            enthalpy: 0.7
+
         - peclet_function_form:
             velocity: classic
             enthalpy: tanh
@@ -169,6 +187,7 @@ realms:
       output_variables:
        - velocity
        - pressure
+       - temperature
 
     restart:
       restart_data_base_name: ablHill3dSymPenalty.rst