diff --git a/include/EnthalpyEquationSystem.h b/include/EnthalpyEquationSystem.h index 3a1fefc6e..49d47b439 100644 --- a/include/EnthalpyEquationSystem.h +++ b/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& /* theNames */, + const std::map>& /* theParams */) final; + + void register_initial_condition_string_function( stk::mesh::Part* part, - const std::map& theNames, - const std::map>& theParams); + const std::map& func) final; void predict_state(); diff --git a/include/Enums.h b/include/Enums.h index 8da40bb3c..e54bc29cc 100644 --- a/include/Enums.h +++ b/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 }; diff --git a/include/EquationSystem.h b/include/EquationSystem.h index cbcb83fe9..a7390afa5 100644 --- a/include/EquationSystem.h +++ b/include/EquationSystem.h @@ -233,6 +233,12 @@ class EquationSystem { } + virtual void register_initial_condition_string_function( + stk::mesh::Part* /*part*/, + const std::map& /*func*/) + { + } + // rip through the propertyAlg_ virtual void evaluate_properties(); diff --git a/include/EquationSystems.h b/include/EquationSystems.h index e8fe1805a..e6c2357c3 100644 --- a/include/EquationSystems.h +++ b/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& func); + void initialize(); void reinitialize_linear_system(); void populate_derived_quantities(); diff --git a/include/NaluParsing.h b/include/NaluParsing.h index 81bca4daf..a29c35e88 100644 --- a/include/NaluParsing.h +++ b/include/NaluParsing.h @@ -37,6 +37,7 @@ struct UserData std::map userFunctionMap_; std::map> functionParams_; std::map> functionStringParams_; + std::map 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> functionParams_; }; +struct StringFunctionInitialConditionData : public InitialCondition +{ + std::map 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& mapName); void operator>>(const YAML::Node& node, std::map& mapName); diff --git a/include/master_element/MasterElementFunctions.h b/include/master_element/MasterElementFunctions.h index aad4bdfc7..cdb0c2bc9 100644 --- a/include/master_element/MasterElementFunctions.h +++ b/include/master_element/MasterElementFunctions.h @@ -77,10 +77,10 @@ generic_grad_op( static_assert( std::is_same::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::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::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::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::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)); diff --git a/include/user_functions/StringTimeCoordFunction.h b/include/user_functions/StringTimeCoordFunction.h new file mode 100644 index 000000000..b5db13792 --- /dev/null +++ b/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 + +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 diff --git a/include/user_functions/StringTimeCoordTemperatureAuxFunction.h b/include/user_functions/StringTimeCoordTemperatureAuxFunction.h new file mode 100644 index 000000000..908fadbcc --- /dev/null +++ b/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 + +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 diff --git a/include/user_functions/CappingInversionTemperatureAuxFunction.h b/include/user_functions/TabulatedTemperatureAuxFunction.h similarity index 68% rename from include/user_functions/CappingInversionTemperatureAuxFunction.h rename to include/user_functions/TabulatedTemperatureAuxFunction.h index f39ec2a07..c33006862 100644 --- a/include/user_functions/CappingInversionTemperatureAuxFunction.h +++ b/include/user_functions/TabulatedTemperatureAuxFunction.h @@ -14,18 +14,15 @@ #include -namespace sierra { -namespace nalu { +namespace sierra::nalu { -class CappingInversionTemperatureAuxFunction : public AuxFunction +class TabulatedTemperatureAuxFunction final : public AuxFunction { public: - CappingInversionTemperatureAuxFunction(); + TabulatedTemperatureAuxFunction( + std::vector heights, std::vector temperatures); - virtual ~CappingInversionTemperatureAuxFunction() {} - - using AuxFunction::do_evaluate; - virtual void do_evaluate( + void do_evaluate( const double* coords, const double time, const unsigned spatialDimension, @@ -33,10 +30,13 @@ class CappingInversionTemperatureAuxFunction : public AuxFunction double* fieldPtr, const unsigned fieldSize, const unsigned beginPos, - const unsigned endPos) const; + const unsigned endPos) const final; + +private: + std::vector heights_; + std::vector temperatures_; }; -} // namespace nalu -} // namespace sierra +} // namespace sierra::nalu #endif diff --git a/reg_tests/test_files/ablHill3dSymPenalty/ablHill3dSymPenalty.norm.gold b/reg_tests/test_files/ablHill3dSymPenalty/ablHill3dSymPenalty.norm.gold index 5c52f03eb..ac354e17a 100644 --- a/reg_tests/test_files/ablHill3dSymPenalty/ablHill3dSymPenalty.norm.gold +++ b/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 diff --git a/reg_tests/test_files/ablHill3dSymPenalty/ablHill3dSymPenalty.yaml b/reg_tests/test_files/ablHill3dSymPenalty/ablHill3dSymPenalty.yaml index da13d8663..c774ff21b 100755 --- a/reg_tests/test_files/ablHill3dSymPenalty/ablHill3dSymPenalty.yaml +++ b/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 diff --git a/src/EnthalpyEquationSystem.C b/src/EnthalpyEquationSystem.C index 14e3a7863..359b243d6 100644 --- a/src/EnthalpyEquationSystem.C +++ b/src/EnthalpyEquationSystem.C @@ -83,7 +83,8 @@ #include #include -#include +#include +#include // overset #include @@ -635,7 +636,7 @@ EnthalpyEquationSystem::register_open_bc( auto& activeKernels = elemSolverAlg->activeKernels_; build_face_topo_kernel_automatic( - partTopo, *this, activeKernels, "turbulent_ke_open", realm_.meta_data(), + partTopo, *this, activeKernels, "temperature_open", realm_.meta_data(), *realm_.solutionOptions_, enthalpy_, enthalpyBc, dataPreReqs); } else { throw std::runtime_error( @@ -1042,11 +1043,12 @@ EnthalpyEquationSystem::register_initial_condition_fcn( } else if (fcnName == "BoussinesqNonIso") { theAuxFunc = new BoussinesqNonIsoTemperatureAuxFunction(); } else if (fcnName == "capping_inversion") { - theAuxFunc = new CappingInversionTemperatureAuxFunction(); + theAuxFunc = new TabulatedTemperatureAuxFunction( + {0, 650.0, 750.0, 1000.0}, {300.0, 300.0, 308.0, 308.75}); } else { - throw std::runtime_error( - "EnthalpyEquationSystem::register_initial_condition_fcn: limited user " - "functions supported"); + throw std::runtime_error("EnthalpyEquationSystem::register_initial_" + "condition_fcn: limited user " + "functions supported"); } // create the algorithm @@ -1058,6 +1060,17 @@ EnthalpyEquationSystem::register_initial_condition_fcn( } } +void +EnthalpyEquationSystem::register_initial_condition_string_function( + stk::mesh::Part* part, const std::map& fcns) +{ + auto func = std::make_unique( + fcns.at("temperature")); + auto auxAlg = std::make_unique( + realm_, part, temperature_, func.release(), stk::topology::NODE_RANK); + realm_.initCondAlg_.push_back(auxAlg.release()); +} + //-------------------------------------------------------------------------- //-------- solve_and_update ------------------------------------------------ //-------------------------------------------------------------------------- @@ -1406,15 +1419,19 @@ EnthalpyEquationSystem::temperature_bc_setup( } else if (fcnName == "BoussinesqNonIso") { theAuxFunc = new BoussinesqNonIsoTemperatureAuxFunction(); } else if (fcnName == "capping_inversion") { - theAuxFunc = new CappingInversionTemperatureAuxFunction(); + theAuxFunc = new TabulatedTemperatureAuxFunction( + {0, 650.0, 750.0, 1000.0}, {300.0, 300.0, 308.0, 308.75}); } else { throw std::runtime_error("EnthalpyEquationSystem::temperature_bc_setup; " "limited user functions supported"); } + } else if (STRING_FUNCTION_UD == theDataType) { + theAuxFunc = new StringTimeCoordTemperatureAuxFunction( + userData.functions.at("temperature")); } else { throw std::runtime_error( - "EnthalpyEquationSystem::temperature_bc_setup: only function and " - "constants supported (and none specified)"); + "EnthalpyEquationSystem::temperature_bc_setup: only function, " + "constants, and string functions supported (and none specified)"); } AuxFunctionAlgorithm* auxTempAlg = new AuxFunctionAlgorithm( diff --git a/src/EquationSystems.C b/src/EquationSystems.C index 231eaffe9..dddf47078 100644 --- a/src/EquationSystems.C +++ b/src/EquationSystems.C @@ -715,6 +715,15 @@ EquationSystems::register_initial_condition_fcn( part, fcnIC.functionNames_, fcnIC.functionParams_); } +void +EquationSystems::register_initial_condition_string_function( + stk::mesh::Part* part, const std::map& func) +{ + // call through to equation systems + for (EquationSystem* eqSys : equationSystemVector_) + eqSys->register_initial_condition_string_function(part, func); +} + //-------------------------------------------------------------------------- //-------- initialize() ---------------------------------------------------- //-------------------------------------------------------------------------- diff --git a/src/InitialConditions.C b/src/InitialConditions.C index b29563479..c1db77ebd 100644 --- a/src/InitialConditions.C +++ b/src/InitialConditions.C @@ -36,6 +36,12 @@ InitialConditionCreator::load_single(const YAML::Node& node) auto* fcnIC = dynamic_cast(ic.get()); node >> *fcnIC; return ic; + } else if (node["string_function"]) { + NaluEnv::self().naluOutputP0() + << "Initial Is Type string-function " << std::endl; + auto string_func = std::make_unique(); + node >> *string_func; + return string_func; } else throw std::runtime_error( "parser error InitialConditions::load; unsupported IC type"); diff --git a/src/NaluParsing.C b/src/NaluParsing.C index c9fdf7ff8..a5b76a136 100644 --- a/src/NaluParsing.C +++ b/src/NaluParsing.C @@ -316,9 +316,9 @@ operator>>( void operator>>(const YAML::Node& node, UserFunctionInitialConditionData& fcnIC) { - fcnIC.theIcType_ = sierra::nalu::FUNCTION_UD; fcnIC.icName_ = node["user_function"].as(); + const YAML::Node& targets = node["target_name"]; if (targets.Type() == YAML::NodeType::Scalar) { fcnIC.targetNames_.resize(1); @@ -343,9 +343,41 @@ operator>>(const YAML::Node& node, UserFunctionInitialConditionData& fcnIC) } void -operator>>(const YAML::Node& node, ConstantInitialConditionData& constIC) +operator>>(const YAML::Node& node, StringFunctionInitialConditionData& fcnIC) { + fcnIC.theIcType_ = sierra::nalu::STRING_FUNCTION_UD; + fcnIC.icName_ = node["string_function"].as(); + const YAML::Node& targets = node["target_name"]; + + if (targets.Type() == YAML::NodeType::Scalar) { + fcnIC.targetNames_.resize(1); + fcnIC.targetNames_[0] = targets.as(); + NaluEnv::self().naluOutputP0() + << "constant IC: name: " << fcnIC.icName_ << " , target[" << 0 + << "] = " << fcnIC.targetNames_[0] << std::endl; + if (fcnIC.targetNames_[0].find(',') != std::string::npos) { + throw std::runtime_error( + "In " + fcnIC.icName_ + + " found ',' in target name - you must enclose in '[...]' for multiple " + "targets"); + } + + } else { + fcnIC.targetNames_.resize(targets.size()); + for (size_t i = 0; i < targets.size(); ++i) { + fcnIC.targetNames_[i] = targets[i].as(); + } + } + + if (expect_map(node, "function", true)) { + fcnIC.functions_ = + node["function"].as>(); + } +} +void +operator>>(const YAML::Node& node, ConstantInitialConditionData& constIC) +{ constIC.theIcType_ = sierra::nalu::CONSTANT_UD; constIC.icName_ = node["constant"].as(); const YAML::Node& targets = node["target_name"]; @@ -750,9 +782,16 @@ convert::decode( sierra::nalu::CONSTANT_UD; } if (node["temperature"]) { - wallData.temperature_ = node["temperature"].as(); + try { + wallData.temperature_ = + node["temperature"].as(); + wallData.bcDataTypeMap_["temperature"] = sierra::nalu::CONSTANT_UD; + } catch (const YAML::BadConversion&) { + wallData.functions.emplace( + "temperature", node["temperature"].as()); + wallData.bcDataTypeMap_["temperature"] = sierra::nalu::STRING_FUNCTION_UD; + } wallData.bcDataSpecifiedMap_["temperature"] = true; - wallData.bcDataTypeMap_["temperature"] = sierra::nalu::CONSTANT_UD; wallData.tempSpec_ = true; } @@ -829,6 +868,7 @@ convert::decode( wallData.bcDataSpecifiedMap_["pressure"] = true; wallData.bcDataTypeMap_["pressure"] = sierra::nalu::CONSTANT_UD; } + if (node["fsi_interface"]) { wallData.isFsiInterface_ = node["fsi_interface"].as(); } @@ -934,8 +974,16 @@ convert::decode( inflowData.massFractionSpec_ = true; } if (node["temperature"]) { - inflowData.temperature_ = - node["temperature"].as(); + try { + inflowData.temperature_ = + node["temperature"].as(); + inflowData.bcDataTypeMap_["temperature"] = sierra::nalu::CONSTANT_UD; + } catch (const YAML::BadConversion&) { + inflowData.functions.emplace( + "temperature", node["temperature"].as()); + inflowData.bcDataTypeMap_["temperature"] = + sierra::nalu::STRING_FUNCTION_UD; + } inflowData.tempSpec_ = true; } @@ -1016,7 +1064,15 @@ convert::decode( openData.massFractionSpec_ = true; } if (node["temperature"]) { - openData.temperature_ = node["temperature"].as(); + try { + openData.temperature_ = + node["temperature"].as(); + openData.bcDataTypeMap_["temperature"] = sierra::nalu::CONSTANT_UD; + } catch (const YAML::BadConversion&) { + openData.functions.emplace( + "temperature", node["temperature"].as()); + openData.bcDataTypeMap_["temperature"] = sierra::nalu::STRING_FUNCTION_UD; + } openData.tempSpec_ = true; } diff --git a/src/Realm.C b/src/Realm.C index f4aaca5f5..f96f62fbd 100644 --- a/src/Realm.C +++ b/src/Realm.C @@ -1327,6 +1327,15 @@ Realm::setup_initial_conditions() equationSystems_.register_initial_condition_fcn(targetPart, fcnIC); } break; + case STRING_FUNCTION_UD: { + const StringFunctionInitialConditionData* fcnIC = + dynamic_cast( + initCond.get()); + assert(fcnIC); + equationSystems_.register_initial_condition_string_function( + targetPart, fcnIC->functions_); + } break; + case USER_SUB_UD: throw std::runtime_error( "Realm::setup_initial_conditions: USER_SUB not supported: "); @@ -5159,7 +5168,7 @@ std::vector Realm::handle_all_element_part_alias( const std::vector& names) const { - if (names.size() == 1u && names.front() == allElementPartAlias) { + if (names.size() == 1u && names.at(0) == allElementPartAlias) { std::vector new_names; for (const auto* part : meta_data().get_mesh_parts()) { ThrowRequire(part); diff --git a/src/master_element/Hex8FEM.C b/src/master_element/Hex8FEM.C index e5fff56b5..090a33d76 100644 --- a/src/master_element/Hex8FEM.C +++ b/src/master_element/Hex8FEM.C @@ -55,9 +55,9 @@ generic_grad_op_3d( static_assert( std::is_same::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, "Weight 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, "Weight view assumed to be 3D"); static_assert(AlgTraits::nDim_ == 3, "3D method"); ThrowAssert(AlgTraits::nodesPerElement_ == referenceGradWeights.extent(1)); diff --git a/src/user_functions/CMakeLists.txt b/src/user_functions/CMakeLists.txt index fa3efa46a..331624d18 100644 --- a/src/user_functions/CMakeLists.txt +++ b/src/user_functions/CMakeLists.txt @@ -4,7 +4,6 @@ target_sources(nalu PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/BoussinesqNonIsoMomentumSrcNodeSuppAlg.C ${CMAKE_CURRENT_SOURCE_DIR}/BoussinesqNonIsoTemperatureAuxFunction.C ${CMAKE_CURRENT_SOURCE_DIR}/BoussinesqNonIsoVelocityAuxFunction.C - ${CMAKE_CURRENT_SOURCE_DIR}/CappingInversionTemperatureAuxFunction.C ${CMAKE_CURRENT_SOURCE_DIR}/ConvectingTaylorVortexPressureAuxFunction.C ${CMAKE_CURRENT_SOURCE_DIR}/ConvectingTaylorVortexVelocityAuxFunction.C ${CMAKE_CURRENT_SOURCE_DIR}/FlowPastCylinderTempAuxFunction.C @@ -19,6 +18,9 @@ target_sources(nalu PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/SteadyTaylorVortexMomentumSrcNodeSuppAlg.C ${CMAKE_CURRENT_SOURCE_DIR}/SteadyTaylorVortexPressureAuxFunction.C ${CMAKE_CURRENT_SOURCE_DIR}/SteadyTaylorVortexVelocityAuxFunction.C + ${CMAKE_CURRENT_SOURCE_DIR}/StringTimeCoordFunction.C + ${CMAKE_CURRENT_SOURCE_DIR}/StringTimeCoordTemperatureAuxFunction.C + ${CMAKE_CURRENT_SOURCE_DIR}/TabulatedTemperatureAuxFunction.C ${CMAKE_CURRENT_SOURCE_DIR}/TaylorGreenPressureAuxFunction.C ${CMAKE_CURRENT_SOURCE_DIR}/TaylorGreenVelocityAuxFunction.C ${CMAKE_CURRENT_SOURCE_DIR}/TornadoAuxFunction.C diff --git a/src/user_functions/CappingInversionTemperatureAuxFunction.C b/src/user_functions/CappingInversionTemperatureAuxFunction.C deleted file mode 100644 index e5a2942aa..000000000 --- a/src/user_functions/CappingInversionTemperatureAuxFunction.C +++ /dev/null @@ -1,66 +0,0 @@ -// 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 -#include -#include - -// basic c++ -#include -#include -#include - -namespace sierra { -namespace nalu { - -CappingInversionTemperatureAuxFunction::CappingInversionTemperatureAuxFunction() - : AuxFunction(0, 1) -{ - // does nothing -} - -void -CappingInversionTemperatureAuxFunction::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 z = coords[2]; - - // heights: [ 0, 650.0, 750.0, 1000.0 ] - // values: [300.0, 300.0, 308.0, 308.75] - - const double slope_1 = (308.0 - 300.0) / (750.0 - 650.0); - const double slope_2 = (308.75 - 308.0) / (1000.0 - 750.0); - - double temp = 300.0; - if (z > 650.0 && z <= 750.0) { - temp = 300.0 + slope_1 * (z - 650.0); - } else if (z > 750.0 && z <= 1000.0) { - temp = 308.0 + slope_2 * (z - 750.0); - } else if (z > 1000.0) { - temp = 308.75; - } - - fieldPtr[0] = temp; - - fieldPtr += fieldSize; - coords += spatialDimension; - } -} - -} // namespace nalu -} // namespace sierra diff --git a/src/user_functions/StringTimeCoordFunction.C b/src/user_functions/StringTimeCoordFunction.C new file mode 100644 index 000000000..c743c8656 --- /dev/null +++ b/src/user_functions/StringTimeCoordFunction.C @@ -0,0 +1,102 @@ +// 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/StringTimeCoordFunction.h" +#include "stk_expreval/Evaluator.hpp" + +#include + +namespace sierra::nalu { + +StringTimeCoordFunction::StringTimeCoordFunction(std::string fcn) +{ + stk::expreval::Eval eval; + + if (fcn.empty()) { + throw std::runtime_error("Empty string function"); + } + + fcn.erase(std::remove_if(fcn.begin(), fcn.end(), ::isspace), fcn.end()); + std::transform(fcn.begin(), fcn.end(), fcn.begin(), ::toupper); + + try { + eval.parse(fcn); + } catch (std::exception& e) { + std::ostringstream os; + os << "unable to parse input string '" << fcn << "'" << std::endl; + throw std::runtime_error(os.str()); + } + + if (eval.undefinedFunction()) { + std::ostringstream os; + os << "found an undefined function in '" << fcn << "'" << std::endl; + throw std::runtime_error(os.str()); + } + + auto bind_index = [&eval](const char* name) { + const auto& names = eval.get_independent_variable_names(); + if (std::find(names.begin(), names.end(), name) != names.end()) { + return eval.get_variable_index(name); + } + return fcn::UNMAPPED_INDEX; + }; + + t_index = bind_index("T"); + x_index = bind_index("X"); + y_index = bind_index("Y"); + z_index = bind_index("Z"); + + std::vector valid_entries{"T", "X", "Y", "Z"}; + for (const auto& name : eval.get_independent_variable_names()) { + if ( + std::find(valid_entries.begin(), valid_entries.end(), name) == + valid_entries.end()) { + throw std::runtime_error( + "Invalid input name in DeviceStringFunction " + name); + } + } + + constant = eval.is_constant_expression(); + parsed_eval = eval.get_parsed_eval(); + try { + eval.evaluate(); + (*this)(0, 0, 0, 0); + } catch (std::exception& e) { + std::ostringstream sout; + sout << "String function was unable to be evaluated with input string '" + << fcn << "'.\n" + << e.what(); + throw std::runtime_error(sout.str()); + } +} + +template +KOKKOS_FUNCTION void +bind_if_valid( + stk::expreval::DeviceVariableMap& var_map, int index, double val) +{ + if (index >= 0 && index < N) { + var_map[index] = val; + } +} + +KOKKOS_FUNCTION double +StringTimeCoordFunction::operator()( + double t, double x, double y, double z) const +{ + constexpr int num_vars = 4; + auto var_map = stk::expreval::DeviceVariableMap(parsed_eval); + bind_if_valid(var_map, t_index, t); + bind_if_valid(var_map, x_index, x); + bind_if_valid(var_map, y_index, y); + bind_if_valid(var_map, z_index, z); + return parsed_eval.evaluate(var_map); +} + +} // namespace sierra::nalu diff --git a/src/user_functions/StringTimeCoordTemperatureAuxFunction.C b/src/user_functions/StringTimeCoordTemperatureAuxFunction.C new file mode 100644 index 000000000..4dc6bf6b1 --- /dev/null +++ b/src/user_functions/StringTimeCoordTemperatureAuxFunction.C @@ -0,0 +1,50 @@ +// 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 +#include + +namespace sierra::nalu { + +StringTimeCoordTemperatureAuxFunction::StringTimeCoordTemperatureAuxFunction( + std::string fcn) + : AuxFunction(0, 1), f_(fcn) +{ +} + +void +StringTimeCoordTemperatureAuxFunction::do_evaluate( + const double* coords, + const double time, + const unsigned dim, + const unsigned numPoints, + double* temperatures, + const unsigned /*fieldSize*/, + const unsigned /*beginPos*/, + const unsigned /*endPos*/) const +{ + if (f_.spatial_dim() > int(dim)) { + // e.g. if someone has "z" coord in 2D + throw std::runtime_error("Dimensional arguments to string function greater " + "than simulation dimension"); + } + + if (dim == 3) { + for (unsigned p = 0; p < numPoints; ++p) { + temperatures[p] = + f_(time, coords[3 * p + 0], coords[3 * p + 1], coords[3 * p + 2]); + } + } else if (dim == 2) { + for (unsigned p = 0; p < numPoints; ++p) { + temperatures[p] = f_(time, coords[2 * p + 0], coords[2 * p + 1], 0); + } + } +} + +} // namespace sierra::nalu diff --git a/src/user_functions/TabulatedTemperatureAuxFunction.C b/src/user_functions/TabulatedTemperatureAuxFunction.C new file mode 100644 index 000000000..7d0b1ebc9 --- /dev/null +++ b/src/user_functions/TabulatedTemperatureAuxFunction.C @@ -0,0 +1,46 @@ +// 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 +#include +#include +#include "utils/LinearInterpolation.h" + +#include +// basic c++ +#include +#include +#include + +namespace sierra::nalu { + +TabulatedTemperatureAuxFunction::TabulatedTemperatureAuxFunction( + std::vector heights, std::vector temperatures) + : AuxFunction(0, 1), heights_(heights), temperatures_(temperatures) +{ +} + +void +TabulatedTemperatureAuxFunction::do_evaluate( + const double* coords, + const double /*time*/, + const unsigned dim, + const unsigned numPoints, + double* temperatures, + const unsigned /*fieldSize*/, + const unsigned /*beginPos*/, + const unsigned /*endPos*/) const +{ + for (unsigned p = 0; p < numPoints; ++p) { + const double height = coords[dim * p + 2]; + utils::linear_interp(heights_, temperatures_, height, temperatures[p]); + } +} + +} // namespace sierra::nalu diff --git a/src/user_functions/ZalesakSphereMassFlowRateKernel.C b/src/user_functions/ZalesakSphereMassFlowRateKernel.C index 268cb9152..2aadf9281 100644 --- a/src/user_functions/ZalesakSphereMassFlowRateKernel.C +++ b/src/user_functions/ZalesakSphereMassFlowRateKernel.C @@ -39,7 +39,6 @@ ZalesakSphereMassFlowRateEdgeAlg::ZalesakSphereMassFlowRateEdgeAlg( void ZalesakSphereMassFlowRateEdgeAlg::execute() { - const double eps = 1.0e-16; const int ndim = realm_.meta_data().spatial_dimension(); // Defaults from AMR-Wind @@ -56,7 +55,7 @@ ZalesakSphereMassFlowRateEdgeAlg::execute() run_algorithm( realm_.bulk_data(), KOKKOS_LAMBDA( - ShmemDataType & smdata, const stk::mesh::FastMeshIndex& edge, + ShmemDataType& /*smdata*/, const stk::mesh::FastMeshIndex& edge, const stk::mesh::FastMeshIndex& nodeL, const stk::mesh::FastMeshIndex& nodeR) { // Scratch work array for edgeAreaVector diff --git a/src/utils/ComputeVectorDivergence.C b/src/utils/ComputeVectorDivergence.C index b7c546ed6..631c0c9e1 100644 --- a/src/utils/ComputeVectorDivergence.C +++ b/src/utils/ComputeVectorDivergence.C @@ -302,7 +302,7 @@ void compute_edge_scalar_divergence( stk::mesh::BulkData& bulk, stk::mesh::PartVector& partVec, - stk::mesh::PartVector& bndyPartVec, + stk::mesh::PartVector& /*bndyPartVec*/, GenericFieldType* faceField, stk::mesh::FieldBase* scalarField) { diff --git a/unit_tests/CMakeLists.txt b/unit_tests/CMakeLists.txt index 1c19d298a..2e4ccf886 100644 --- a/unit_tests/CMakeLists.txt +++ b/unit_tests/CMakeLists.txt @@ -38,7 +38,10 @@ target_sources(${utest_ex_name} PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestSimdBasic.C ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestSingleHexPromotion.C ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestSpinnerLidarPattern.C + ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestStringTimeCoordFunction.C + ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestStringTimeCoordTemperatureAuxFunction.C ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestSuppAlgDataSharing.C + ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestTabulatedTemperatureAuxFunction.C ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestUtils.C ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestVSpace.C ) diff --git a/unit_tests/UnitTestStringTimeCoordFunction.C b/unit_tests/UnitTestStringTimeCoordFunction.C new file mode 100644 index 000000000..49e879b2d --- /dev/null +++ b/unit_tests/UnitTestStringTimeCoordFunction.C @@ -0,0 +1,95 @@ +#include "gtest/gtest.h" + +#include "user_functions/StringTimeCoordFunction.h" + +namespace sierra::nalu { +TEST(StringTimeCoordFunction, same_on_host_as_string_function) +{ + std::string func = "1 + 3 * t + x + y + z*y"; + const double val_at_1234 = 21; + ASSERT_DOUBLE_EQ(StringTimeCoordFunction(func)(1, 2, 3, 4), val_at_1234); +} + +double +val_executed_in_kernel( + StringTimeCoordFunction func, double t, double x, double y, double z) +{ + double val = 0; + Kokkos::parallel_reduce( + 1, KOKKOS_LAMBDA(int, double& v) { v += func(t, x, y, z); }, + Kokkos::Sum(val)); + return val; +} + +TEST(StringTimeCoordFunction, runs_on_device) +{ + std::string func = + "1 + 3 * sin(t) + exp(x) + log(y) + ((z>0 && z < 1) ? z*y : x/z)"; + + StringTimeCoordFunction f(func); + auto val_1 = f(1, 2, 3, 4); + + auto val_2 = val_executed_in_kernel(f, 1, 2, 3, 4); + ASSERT_DOUBLE_EQ(val_1, val_2); +} + +TEST(StringTimeCoordFunction, empty) +{ + EXPECT_ANY_THROW(StringTimeCoordFunction f("")); +} + +TEST(StringTimeCoordFunction, unparsable) +{ + EXPECT_ANY_THROW(StringTimeCoordFunction f("c++/45g")); +} + +TEST(StringTimeCoordFunction, undefined) +{ + EXPECT_ANY_THROW(StringTimeCoordFunction f("h(t) + 34")); +} + +TEST(StringTimeCoordFunction, copy) +{ + std::string s = "3*t + x*y - z"; + StringTimeCoordFunction f(s); + + StringTimeCoordFunction f2(s); + StringTimeCoordFunction f3(f2); + StringTimeCoordFunction f4 = f2; + + const double t = 2.5; + const double x = 1.1; + const double y = -1.5; + const double z = 12.5; + const double ans = 7.5 - 1.1 * 1.5 - 12.5; + EXPECT_DOUBLE_EQ(val_executed_in_kernel(f2, t, x, y, z), ans); + EXPECT_DOUBLE_EQ(val_executed_in_kernel(f3, t, x, y, z), ans); + EXPECT_DOUBLE_EQ(val_executed_in_kernel(f2, t, x, y, z), ans); + EXPECT_DOUBLE_EQ(val_executed_in_kernel(f4, t, x, y, z), ans); + + EXPECT_TRUE(f2.is_spatial()); +} + +TEST(StringTimeCoordFunction, time_only) +{ + StringTimeCoordFunction f = {"3*t"}; + + const double t = 2.5; + const double x = 1.1; + const double y = -1.5; + const double z = 12.5; + const double ans = 7.5; + + EXPECT_DOUBLE_EQ(val_executed_in_kernel(f, t, x, y, z), ans); + EXPECT_FALSE(f.is_spatial()); +} + +TEST(StringTimeCoordFunction, spatial_dim) +{ + ASSERT_EQ(StringTimeCoordFunction{"3*t + x"}.spatial_dim(), 1); + ASSERT_EQ(StringTimeCoordFunction{"y"}.spatial_dim(), 2); + ASSERT_EQ(StringTimeCoordFunction{"x + z"}.spatial_dim(), 3); + ASSERT_EQ(StringTimeCoordFunction{"y+z"}.spatial_dim(), 3); +} + +} // namespace sierra::nalu \ No newline at end of file diff --git a/unit_tests/UnitTestStringTimeCoordTemperatureAuxFunction.C b/unit_tests/UnitTestStringTimeCoordTemperatureAuxFunction.C new file mode 100644 index 000000000..64c3696db --- /dev/null +++ b/unit_tests/UnitTestStringTimeCoordTemperatureAuxFunction.C @@ -0,0 +1,36 @@ +#include "gtest/gtest.h" + +#include "user_functions/StringTimeCoordTemperatureAuxFunction.h" + +namespace sierra::nalu { +TEST(StringTimeCoordTemperatureAuxFunction, only_z) +{ + StringTimeCoordTemperatureAuxFunction func("z"); + + double coords[3] = {1, 2, 3}; + double temperature; + func.do_evaluate(coords, 0., 3, 1, &temperature, 1, 0, 1); + ASSERT_DOUBLE_EQ(temperature, 3.); +} + +TEST(StringTimeCoordTemperatureAuxFunction, only_y) +{ + StringTimeCoordTemperatureAuxFunction func("y"); + + double coords[3] = {1, 2, 1}; + double temperature; + func.do_evaluate(coords, 0., 3, 1, &temperature, 1, 0, 1); + ASSERT_DOUBLE_EQ(temperature, 2.); +} + +TEST(StringTimeCoordTemperatureAuxFunction, full) +{ + StringTimeCoordTemperatureAuxFunction func("log(1+t) + exp(x)*log(y) + z"); + + double coords[3] = {1, 2, 1}; + double temperature; + func.do_evaluate(coords, 0., 3, 1, &temperature, 1, 0, 1); + ASSERT_DOUBLE_EQ(temperature, 1 + exp(1) * log(2)); +} + +} // namespace sierra::nalu \ No newline at end of file diff --git a/unit_tests/UnitTestTabulatedTemperatureAuxFunction.C b/unit_tests/UnitTestTabulatedTemperatureAuxFunction.C new file mode 100644 index 000000000..d87ada064 --- /dev/null +++ b/unit_tests/UnitTestTabulatedTemperatureAuxFunction.C @@ -0,0 +1,46 @@ +#include "gtest/gtest.h" + +#include "user_functions/TabulatedTemperatureAuxFunction.h" + +namespace sierra::nalu { +TEST(TabulatedTemperatureAuxFunction, end_of_range) +{ + TabulatedTemperatureAuxFunction func({0, 1, 2}, {2, 3, 4}); + + double coords[3] = {1, 2, 3}; + double temperature; + func.do_evaluate(coords, 0., 3, 1, &temperature, 1, 0, 1); + ASSERT_DOUBLE_EQ(temperature, 4.); +} + +TEST(TabulatedTemperatureAuxFunction, matches_at_table_entry) +{ + TabulatedTemperatureAuxFunction func({0, 1, 2}, {2, 3, 4}); + + double coords[3] = {1, 2, 1}; + double temperature; + func.do_evaluate(coords, 0., 3, 1, &temperature, 1, 0, 1); + ASSERT_DOUBLE_EQ(temperature, 3.); +} + +TEST(TabulatedTemperatureAuxFunction, averages_halfway_between_points) +{ + TabulatedTemperatureAuxFunction func({0, 1, 2}, {2, 3, 4}); + + double coords[3] = {1, 2, 1.5}; + double temperature; + func.do_evaluate(coords, 0., 3, 1, &temperature, 1, 0, 1); + ASSERT_DOUBLE_EQ(temperature, (3 + 4) / 2.); +} + +TEST(TabulatedTemperatureAuxFunction, before_ranges) +{ + TabulatedTemperatureAuxFunction func({0, 1, 2}, {2, 3, 4}); + + double coords[3] = {1, 2, -1}; + double temperature; + func.do_evaluate(coords, 0., 3, 1, &temperature, 1, 0, 1); + ASSERT_DOUBLE_EQ(temperature, 2); +} + +} // namespace sierra::nalu \ No newline at end of file