From e1084bf3cf84338626e8c0cd5da647c69e273595 Mon Sep 17 00:00:00 2001 From: Marc Henry de Frahan Date: Wed, 13 Jul 2022 11:42:48 -0600 Subject: [PATCH] Add k-eps and k-omega RANS models --- docs/references/references.bib | 18 + docs/source/theory/supportedEquationSet.rst | 24 +- docs/source/theory/turbulenceModeling.rst | 10 +- include/ChienKEpsilonEquationSystem.h | 89 +++ include/Enums.h | 74 ++- include/NaluParsedTypes.h | 6 + include/NaluParsing.h | 32 +- include/ShearStressTransportEquationSystem.h | 6 +- include/TotalDissipationRateEquationSystem.h | 93 +++ include/WilcoxKOmegaEquationSystem.h | 87 +++ include/ngp_algorithms/TurbViscKEAlg.h | 53 ++ include/ngp_algorithms/TurbViscKOAlg.h | 55 ++ include/node_kernels/SDRKONodeKernel.h | 79 +++ include/node_kernels/SDRSSTDESNodeKernel.h | 1 - include/node_kernels/SDRSSTNodeKernel.h | 1 - include/node_kernels/TDRKENodeKernel.h | 75 +++ include/node_kernels/TKEKENodeKernel.h | 71 ++ include/node_kernels/TKEKONodeKernel.h | 69 ++ include/node_kernels/TKEKsgsNodeKernel.h | 1 - include/node_kernels/TKESSTDESNodeKernel.h | 1 - include/node_kernels/TKESSTIDDESNodeKernel.h | 1 - include/node_kernels/TKESSTNodeKernel.h | 1 - reg_tests/CTestList.cmake | 2 + .../KEChannelEdge/KEChannelEdge.yaml | 200 ++++++ .../KOChannelEdge/KOChannelEdge.yaml | 203 ++++++ src/CMakeLists.txt | 3 + src/ChienKEpsilonEquationSystem.C | 437 +++++++++++++ src/EquationSystems.C | 12 + src/LowMachEquationSystem.C | 10 + src/NaluParsing.C | 48 +- src/ShearStressTransportEquationSystem.C | 20 +- src/SolutionOptions.C | 7 + src/SpecificDissipationRateEquationSystem.C | 26 +- src/TotalDissipationRateEquationSystem.C | 614 ++++++++++++++++++ src/TurbKineticEnergyEquationSystem.C | 30 +- src/WilcoxKOmegaEquationSystem.C | 394 +++++++++++ src/ngp_algorithms/CMakeLists.txt | 2 + src/ngp_algorithms/TurbViscKEAlg.C | 75 +++ src/ngp_algorithms/TurbViscKOAlg.C | 99 +++ src/node_kernels/CMakeLists.txt | 4 + src/node_kernels/SDRKONodeKernel.C | 142 ++++ src/node_kernels/SDRSSTDESNodeKernel.C | 3 - src/node_kernels/SDRSSTNodeKernel.C | 3 - src/node_kernels/TDRKENodeKernel.C | 102 +++ src/node_kernels/TKEKENodeKernel.C | 93 +++ src/node_kernels/TKEKONodeKernel.C | 91 +++ src/node_kernels/TKEKsgsNodeKernel.C | 3 - src/node_kernels/TKESSTDESNodeKernel.C | 3 - src/node_kernels/TKESSTIDDESNodeKernel.C | 3 - src/node_kernels/TKESSTNodeKernel.C | 3 - unit_tests/kernels/UnitTestKernelUtils.C | 41 ++ unit_tests/kernels/UnitTestKernelUtils.h | 135 ++++ unit_tests/ngp_algorithms/CMakeLists.txt | 2 + .../ngp_algorithms/UnitTestTurbViscKEAlg.C | 63 ++ .../ngp_algorithms/UnitTestTurbViscKOAlg.C | 63 ++ unit_tests/node_kernels/CMakeLists.txt | 2 + .../node_kernels/UnitTestKENodeKernel.C | 124 ++++ .../node_kernels/UnitTestKONodeKernel.C | 124 ++++ 58 files changed, 3932 insertions(+), 101 deletions(-) create mode 100644 include/ChienKEpsilonEquationSystem.h create mode 100644 include/TotalDissipationRateEquationSystem.h create mode 100644 include/WilcoxKOmegaEquationSystem.h create mode 100644 include/ngp_algorithms/TurbViscKEAlg.h create mode 100644 include/ngp_algorithms/TurbViscKOAlg.h create mode 100644 include/node_kernels/SDRKONodeKernel.h create mode 100644 include/node_kernels/TDRKENodeKernel.h create mode 100644 include/node_kernels/TKEKENodeKernel.h create mode 100644 include/node_kernels/TKEKONodeKernel.h create mode 100644 reg_tests/test_files/KEChannelEdge/KEChannelEdge.yaml create mode 100644 reg_tests/test_files/KOChannelEdge/KOChannelEdge.yaml create mode 100644 src/ChienKEpsilonEquationSystem.C create mode 100644 src/TotalDissipationRateEquationSystem.C create mode 100644 src/WilcoxKOmegaEquationSystem.C create mode 100644 src/ngp_algorithms/TurbViscKEAlg.C create mode 100644 src/ngp_algorithms/TurbViscKOAlg.C create mode 100644 src/node_kernels/SDRKONodeKernel.C create mode 100644 src/node_kernels/TDRKENodeKernel.C create mode 100644 src/node_kernels/TKEKENodeKernel.C create mode 100644 src/node_kernels/TKEKONodeKernel.C create mode 100644 unit_tests/ngp_algorithms/UnitTestTurbViscKEAlg.C create mode 100644 unit_tests/ngp_algorithms/UnitTestTurbViscKOAlg.C create mode 100644 unit_tests/node_kernels/UnitTestKENodeKernel.C create mode 100644 unit_tests/node_kernels/UnitTestKONodeKernel.C diff --git a/docs/references/references.bib b/docs/references/references.bib index bbb283ff56..1ed40a513f 100644 --- a/docs/references/references.bib +++ b/docs/references/references.bib @@ -395,3 +395,21 @@ @misc{Adcock:2021 howpublished = {APS DFD}, year = {2021}, } + +@article{chien1982predictions, + title={Predictions of channel and boundary-layer flows with a low-Reynolds-number turbulence model}, + author={Chien, Kuei-Yuan}, + journal={AIAA journal}, + volume={20}, + number={1}, + pages={33--38}, + year={1982} +} + +@book{wilcox1998turbulence, + title={Turbulence modeling for CFD}, + author={Wilcox, David C and others}, + volume={2}, + year={1998}, + publisher={DCW industries La Canada, CA} +} \ No newline at end of file diff --git a/docs/source/theory/supportedEquationSet.rst b/docs/source/theory/supportedEquationSet.rst index f6df0f442c..430586fe87 100644 --- a/docs/source/theory/supportedEquationSet.rst +++ b/docs/source/theory/supportedEquationSet.rst @@ -685,14 +685,22 @@ For simulations in which a buoyancy source term is desired, the code supports th .. _eqn_komega_sst: -Shear Stress Transport (SST) RANS Model Suite -+++++++++++++++++++++++++++++++++++++++++++++ - -SST Formulation -~~~~~~~~~~~~~~~ - -Although Nalu-Wind is primarily expected to be a LES simulation tool, RANS -modeling is supported through the activation of the SST equation set. +RANS Model Suite +++++++++++++++++ + +Although Nalu-Wind is primarily expected to be a LES simulation tool, +RANS modeling is supported through the activation of different +two-equation RANS models: the Chien :math:`k-\epsilon` model +:cite:`chien1982predictions`, the Wilcox 1998 :math:`k-\omega` model +:cite:`wilcox1998turbulence`, and the SST model. For the first two +models, the reader is referred to the reference papers and the NASA +Turbulence Modeling Resource for the `Chien +`_ and `Wilcox +`_ models. The SST model +is explained in more details below. + +Shear Stress Transport (SST) Formulation +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ It has been observed that standard 1998 :math:`k-\omega` models display a strong sensitivity to the free stream value of :math:`\omega` (see diff --git a/docs/source/theory/turbulenceModeling.rst b/docs/source/theory/turbulenceModeling.rst index 55ac918e8e..33fb0ff0c6 100644 --- a/docs/source/theory/turbulenceModeling.rst +++ b/docs/source/theory/turbulenceModeling.rst @@ -166,14 +166,14 @@ One Equation :math:`k^{sgs}` See :math:`k^{sgs}` :ref:`PDE ` section. -SST RANS Model -++++++++++++++ +RANS Models ++++++++++++ -As noted, Nalu-Wind does support a SST RANS-based model (the reader is -referred to the SST equation set description). +As noted, Nalu-Wind supports several different RANS-based model (the +reader is referred to the RANS model equation set description). Hybrid RANS/LES Models -+++++++++++++++++++++++++++++++++++ +++++++++++++++++++++++ Nalu-Wind supports the Active Model Split (AMS) hybrid RANS/LES turbulence model :cite:`Haering-etal:2020`. The reader is referred to the AMS equation set for more details. diff --git a/include/ChienKEpsilonEquationSystem.h b/include/ChienKEpsilonEquationSystem.h new file mode 100644 index 0000000000..03fda6eacb --- /dev/null +++ b/include/ChienKEpsilonEquationSystem.h @@ -0,0 +1,89 @@ +// 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 ChienKEpsilonEquationSystem_h +#define ChienKEpsilonEquationSystem_h + +#include +#include +#include + +namespace stk { +struct topology; +namespace mesh { +class Part; +} +} // namespace stk + +namespace sierra { +namespace nalu { + +class EquationSystems; +class AlgorithmDriver; +class TurbKineticEnergyEquationSystem; +class TotalDissipationRateEquationSystem; + +class ChienKEpsilonEquationSystem : public EquationSystem +{ + +public: + ChienKEpsilonEquationSystem(EquationSystems& equationSystems); + virtual ~ChienKEpsilonEquationSystem(); + + virtual void load(const YAML::Node&); + + virtual void initialize(); + + virtual void register_nodal_fields(stk::mesh::Part* part); + + virtual void register_wall_bc( + stk::mesh::Part* part, + const stk::topology& theTopo, + const WallBoundaryConditionData& wallBCData); + + virtual void register_interior_algorithm(stk::mesh::Part* part); + + virtual void solve_and_update(); + + void initial_work(); + virtual void post_external_data_transfer_work(); + virtual void post_iter_work(); + + void clip_min_distance_to_wall(); + void compute_dplus_function(); + void update_and_clip(); + void clip_ke( + const stk::mesh::NgpMesh& ngpMesh, + const stk::mesh::Selector& sel, + stk::mesh::NgpField& tke, + stk::mesh::NgpField& tdr); + + std::unique_ptr tkeEqSys_; + std::unique_ptr tdrEqSys_; + + ScalarFieldType* tke_; + ScalarFieldType* tdr_; + ScalarFieldType* minDistanceToWall_; + ScalarFieldType* dplus_; + + bool isInit_; + + // saved of mesh parts that are for wall bcs + std::vector wallBcPart_; + + bool resetAMSAverages_; + + const double tkeMinValue_{1.0e-8}; + const double tdrMinValue_{1.0e-8}; +}; + +} // namespace nalu +} // namespace sierra + +#endif diff --git a/include/Enums.h b/include/Enums.h index 335b4908b0..e856079095 100644 --- a/include/Enums.h +++ b/include/Enums.h @@ -74,15 +74,16 @@ enum EquationType { EQ_ENTHALPY = 6, EQ_MESH_DISPLACEMENT = 7, EQ_SPEC_DISS_RATE = 8, - EQ_MASS_FRACTION = 9, - EQ_PNG = 10, - EQ_PNG_P = 11, - EQ_PNG_Z = 12, - EQ_PNG_H = 13, - EQ_PNG_U = 14, - EQ_PNG_TKE = 15, // FIXME... Last PNG managed like this.. - EQ_WALL_DISTANCE = 16, - EQ_GAMMA_TRANS = 17, + EQ_TOT_DISS_RATE = 9, + EQ_MASS_FRACTION = 10, + EQ_PNG = 11, + EQ_PNG_P = 12, + EQ_PNG_Z = 13, + EQ_PNG_H = 14, + EQ_PNG_U = 15, + EQ_PNG_TKE = 16, // FIXME... Last PNG managed like this.. + EQ_WALL_DISTANCE = 17, + EQ_GAMMA_TRANS = 18, EquationSystemType_END }; @@ -96,6 +97,7 @@ static const std::string EquationTypeMap[] = { "Enthalpy", "MeshVelocity", "Specific_Dissipation_Rate", + "Total_Dissipation_Rate", "Mass_Fraction", "PNG", "PNG_P", @@ -103,8 +105,7 @@ static const std::string EquationTypeMap[] = { "PNG_H", "PNG_U", "PNG_TKE", - "Wall_Distance" -}; + "Wall_Distance"}; enum UserDataType { CONSTANT_UD = 0, @@ -166,19 +167,15 @@ enum class TurbulenceModel { SST_DES, SST_AMS, SST_IDDES, + KE, + KO, TurbulenceModel_END }; // matching string name index into above enums (must match PERFECTLY) static const std::string TurbulenceModelNames[] = { - "laminar", - "ksgs", - "smagorinsky", - "wale", - "sst", - "sst_des", - "sst_ams", - "sst_iddes"}; + "laminar", "ksgs", "smagorinsky", "wale", "sst", "sst_des", + "sst_ams", "sst_iddes", "ke", "ko"}; enum TurbulenceModelConstant { TM_cMu = 0, @@ -234,15 +231,23 @@ enum TurbulenceModelConstant { TM_ams_peclet_offset = 50, TM_ams_peclet_slope = 51, TM_ams_peclet_scale = 52, - TM_tkeAmb = 53, - TM_sdrAmb = 54, - TM_caOne = 55, - TM_caTwo = 56, - TM_ceOne = 57, - TM_ceTwo = 58, - TM_c0t = 59, - TM_avgTimeCoeff = 60, - TM_END = 61 + TM_fMuExp = 53, + TM_utau = 54, + TM_cEpsOne = 55, + TM_cEpsTwo = 56, + TM_fOne = 57, + TM_sigmaK = 58, + TM_sigmaEps = 59, + TM_tkeAmb = 60, + TM_sdrAmb = 61, + TM_avgTimeCoeff = 62, + TM_alphaInf = 63, + TM_caOne = 64, + TM_caTwo = 65, + TM_ceOne = 66, + TM_ceTwo = 67, + TM_c0t = 68, + TM_END = 69 }; static const std::string TurbulenceModelConstantNames[] = { @@ -299,9 +304,22 @@ static const std::string TurbulenceModelConstantNames[] = { "ams_peclet_offset", "ams_peclet_slope", "ams_peclet_scale", + "fMuExp", + "utau", + "cEpsOne", + "cEpsTwo", + "fOne", + "sigmaK", + "sigmaEps", "tke_amb", "sdr_amb", "avgTimeCoeff", + "alphaInf", + "caOne", + "caTwo", + "ceOne", + "ceTwo", + "c0t", "END"}; enum ActuatorType { diff --git a/include/NaluParsedTypes.h b/include/NaluParsedTypes.h index 6e130c87e4..3f9cebcffa 100644 --- a/include/NaluParsedTypes.h +++ b/include/NaluParsedTypes.h @@ -54,6 +54,12 @@ struct SpecDissRate { {} }; +struct TotDissRate +{ + double totDissRate_; + TotDissRate() : totDissRate_(0.0) {} +}; + struct GammaInf { double gamma_; GammaInf() diff --git a/include/NaluParsing.h b/include/NaluParsing.h index 290f8f89af..896254ce5a 100644 --- a/include/NaluParsing.h +++ b/include/NaluParsing.h @@ -53,6 +53,7 @@ struct WallUserData : public UserData { Velocity u_; Velocity dx_; TurbKinEnergy tke_; + TotDissRate tdr_; MixtureFraction mixFrac_; MassFraction massFraction_; NormalHeatFlux q_; @@ -96,6 +97,7 @@ struct InflowUserData : public UserData { Velocity u_; TurbKinEnergy tke_; SpecDissRate sdr_; + TotDissRate tdr_; MixtureFraction mixFrac_; MassFraction massFraction_; GammaInf gamma_; @@ -103,12 +105,19 @@ struct InflowUserData : public UserData { bool uSpec_; bool tkeSpec_; bool sdrSpec_; + bool tdrSpec_; bool mixFracSpec_; bool massFractionSpec_; bool gammaSpec_; InflowUserData() : UserData(), - uSpec_(false), tkeSpec_(false), sdrSpec_(false), mixFracSpec_(false), massFractionSpec_(false), gammaSpec_(false) + uSpec_(false), + tkeSpec_(false), + sdrSpec_(false), + tdrSpec_(false), + mixFracSpec_(false), + massFractionSpec_(false), + gammaSpec_(false) {} }; @@ -117,6 +126,7 @@ struct OpenUserData : public UserData { Pressure p_; TurbKinEnergy tke_; SpecDissRate sdr_; + TotDissRate tdr_; MixtureFraction mixFrac_; MassFraction massFraction_; GammaOpen gamma_; @@ -125,6 +135,7 @@ struct OpenUserData : public UserData { bool pSpec_; bool tkeSpec_; bool sdrSpec_; + bool tdrSpec_; bool mixFracSpec_; bool massFractionSpec_; bool totalP_; @@ -133,9 +144,16 @@ struct OpenUserData : public UserData { OpenUserData() : UserData(), - uSpec_(false), pSpec_(false), tkeSpec_(false), - sdrSpec_(false), mixFracSpec_(false), massFractionSpec_(false), - totalP_{false}, gammaSpec_(false), entrainMethod_{EntrainmentMethod::COMPUTED} + uSpec_(false), + pSpec_(false), + tkeSpec_(false), + sdrSpec_(false), + tdrSpec_(false), + mixFracSpec_(false), + massFractionSpec_(false), + totalP_{false}, + gammaSpec_(false), + entrainMethod_{EntrainmentMethod::COMPUTED} {} }; @@ -450,6 +468,12 @@ template<> struct convert { static bool decode(const Node& node, sierra::nalu::SpecDissRate& rhs) ; }; +template <> +struct convert +{ + static bool decode(const Node& node, sierra::nalu::TotDissRate& rhs); +}; + template<> struct convert { static bool decode(const Node& node, sierra::nalu::GammaInf& rhs) ; }; diff --git a/include/ShearStressTransportEquationSystem.h b/include/ShearStressTransportEquationSystem.h index 00b14164de..4952264a38 100644 --- a/include/ShearStressTransportEquationSystem.h +++ b/include/ShearStressTransportEquationSystem.h @@ -71,9 +71,9 @@ class ShearStressTransportEquationSystem : public EquationSystem const stk::mesh::Selector& sel, stk::mesh::NgpField& gamma); - TurbKineticEnergyEquationSystem* tkeEqSys_; - SpecificDissipationRateEquationSystem* sdrEqSys_; - GammaEquationSystem *gammaEqSys_; + std::unique_ptr tkeEqSys_; + std::unique_ptr sdrEqSys_; + std::unique_ptr gammaEqSys_; ScalarFieldType* tke_; ScalarFieldType* sdr_; diff --git a/include/TotalDissipationRateEquationSystem.h b/include/TotalDissipationRateEquationSystem.h new file mode 100644 index 0000000000..6f715f5f5d --- /dev/null +++ b/include/TotalDissipationRateEquationSystem.h @@ -0,0 +1,93 @@ +// 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 TotalDissipationRateEquationSystem_h +#define TotalDissipationRateEquationSystem_h + +#include +#include +#include + +#include "ngp_algorithms/NodalGradAlgDriver.h" + +namespace stk { +struct topology; +} + +namespace sierra { +namespace nalu { + +class Realm; +class LinearSystem; +class EquationSystems; + +class TotalDissipationRateEquationSystem : public EquationSystem +{ + +public: + TotalDissipationRateEquationSystem(EquationSystems& equationSystems); + virtual ~TotalDissipationRateEquationSystem(); + + virtual void register_nodal_fields(stk::mesh::Part* part); + + void register_interior_algorithm(stk::mesh::Part* part); + + void register_inflow_bc( + stk::mesh::Part* part, + const stk::topology& theTopo, + const InflowBoundaryConditionData& inflowBCData); + + void register_open_bc( + stk::mesh::Part* part, + const stk::topology& theTopo, + const OpenBoundaryConditionData& openBCData); + + void register_wall_bc( + stk::mesh::Part* part, + const stk::topology& theTopo, + const WallBoundaryConditionData& wallBCData); + + virtual void register_symmetry_bc( + stk::mesh::Part* part, + const stk::topology& theTopo, + const SymmetryBoundaryConditionData& symmetryBCData); + + virtual void register_non_conformal_bc( + stk::mesh::Part* part, const stk::topology& theTopo); + + virtual void register_overset_bc(); + + void initialize(); + void reinitialize_linear_system(); + + void predict_state(); + void assemble_nodal_gradient(); + void compute_effective_diff_flux_coeff(); + + const bool managePNG_; + + ScalarFieldType* tdr_; + VectorFieldType* dedx_; + ScalarFieldType* eTmp_; + ScalarFieldType* visc_; + ScalarFieldType* tvisc_; + ScalarFieldType* evisc_; + ScalarFieldType* tdrWallBc_; + ScalarFieldType* assembledWallTdr_; + ScalarFieldType* assembledWallArea_; + + ScalarNodalGradAlgDriver nodalGradAlgDriver_; + std::unique_ptr effDiffFluxAlg_; + std::unique_ptr wallModelAlgDriver_; +}; + +} // namespace nalu +} // namespace sierra + +#endif diff --git a/include/WilcoxKOmegaEquationSystem.h b/include/WilcoxKOmegaEquationSystem.h new file mode 100644 index 0000000000..f2f64ca485 --- /dev/null +++ b/include/WilcoxKOmegaEquationSystem.h @@ -0,0 +1,87 @@ +// 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 WilcoxKOmegaEquationSystem_h +#define WilcoxKOmegaEquationSystem_h + +#include +#include +#include + +namespace stk { +struct topology; +namespace mesh { +class Part; +} +} // namespace stk + +namespace sierra { +namespace nalu { + +class EquationSystems; +class AlgorithmDriver; +class TurbKineticEnergyEquationSystem; +class SpecificDissipationRateEquationSystem; + +class WilcoxKOmegaEquationSystem : public EquationSystem +{ + +public: + WilcoxKOmegaEquationSystem(EquationSystems& equationSystems); + virtual ~WilcoxKOmegaEquationSystem(); + + virtual void load(const YAML::Node&); + + virtual void initialize(); + + virtual void register_nodal_fields(stk::mesh::Part* part); + + virtual void register_wall_bc( + stk::mesh::Part* part, + const stk::topology& theTopo, + const WallBoundaryConditionData& wallBCData); + + virtual void register_interior_algorithm(stk::mesh::Part* part); + + virtual void solve_and_update(); + + void initial_work(); + virtual void post_external_data_transfer_work(); + virtual void post_iter_work(); + + void clip_min_distance_to_wall(); + void update_and_clip(); + void clip_ko( + const stk::mesh::NgpMesh& ngpMesh, + const stk::mesh::Selector& sel, + stk::mesh::NgpField& tke, + stk::mesh::NgpField& sdr); + + std::unique_ptr tkeEqSys_; + std::unique_ptr sdrEqSys_; + + ScalarFieldType* tke_; + ScalarFieldType* sdr_; + ScalarFieldType* minDistanceToWall_; + + bool isInit_; + + // saved of mesh parts that are for wall bcs + std::vector wallBcPart_; + + bool resetAMSAverages_; + + const double tkeMinValue_{1.0e-8}; + const double sdrMinValue_{1.0e-8}; +}; + +} // namespace nalu +} // namespace sierra + +#endif diff --git a/include/ngp_algorithms/TurbViscKEAlg.h b/include/ngp_algorithms/TurbViscKEAlg.h new file mode 100644 index 0000000000..9773933c58 --- /dev/null +++ b/include/ngp_algorithms/TurbViscKEAlg.h @@ -0,0 +1,53 @@ +// 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 TurbViscKEAlg_h +#define TurbViscKEAlg_h + +#include "Algorithm.h" +#include "FieldTypeDef.h" + +#include "stk_mesh/base/Types.hpp" + +namespace sierra { +namespace nalu { + +class Realm; + +class TurbViscKEAlg : public Algorithm +{ +public: + using DblType = double; + + TurbViscKEAlg( + Realm& realm, + stk::mesh::Part* part, + ScalarFieldType* tvisc, + const bool = false); + + virtual ~TurbViscKEAlg() = default; + + virtual void execute() override; + +private: + ScalarFieldType* tviscField_{nullptr}; + unsigned density_{stk::mesh::InvalidOrdinal}; + unsigned tke_{stk::mesh::InvalidOrdinal}; + unsigned tdr_{stk::mesh::InvalidOrdinal}; + unsigned tvisc_{stk::mesh::InvalidOrdinal}; + unsigned dplus_{stk::mesh::InvalidOrdinal}; + + const DblType cMu_; + const DblType fMuExp_; +}; + +} // namespace nalu +} // namespace sierra + +#endif diff --git a/include/ngp_algorithms/TurbViscKOAlg.h b/include/ngp_algorithms/TurbViscKOAlg.h new file mode 100644 index 0000000000..6479cc3f69 --- /dev/null +++ b/include/ngp_algorithms/TurbViscKOAlg.h @@ -0,0 +1,55 @@ +// 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 TurbViscKOAlg_h +#define TurbViscKOAlg_h + +#include "Algorithm.h" +#include "FieldTypeDef.h" + +#include "stk_mesh/base/Types.hpp" + +namespace sierra { +namespace nalu { + +class Realm; + +class TurbViscKOAlg : public Algorithm +{ +public: + using DblType = double; + + TurbViscKOAlg( + Realm& realm, + stk::mesh::Part* part, + ScalarFieldType* tvisc, + const bool = false); + + virtual ~TurbViscKOAlg() = default; + + virtual void execute() override; + +private: + ScalarFieldType* tviscField_{nullptr}; + unsigned density_{stk::mesh::InvalidOrdinal}; + unsigned viscosity_{stk::mesh::InvalidOrdinal}; + unsigned tke_{stk::mesh::InvalidOrdinal}; + unsigned sdr_{stk::mesh::InvalidOrdinal}; + unsigned minDistance_{stk::mesh::InvalidOrdinal}; + unsigned dudx_{stk::mesh::InvalidOrdinal}; + unsigned tvisc_{stk::mesh::InvalidOrdinal}; + + const DblType aOne_; + const DblType betaStar_; +}; + +} // namespace nalu +} // namespace sierra + +#endif diff --git a/include/node_kernels/SDRKONodeKernel.h b/include/node_kernels/SDRKONodeKernel.h new file mode 100644 index 0000000000..337218f3c9 --- /dev/null +++ b/include/node_kernels/SDRKONodeKernel.h @@ -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. +// + +#ifndef SDRKONODEKERNEL_H +#define SDRKONODEKERNEL_H + +#include "node_kernels/NodeKernel.h" +#include "FieldTypeDef.h" + +#include "stk_mesh/base/BulkData.hpp" +#include "stk_mesh/base/Ngp.hpp" +#include "stk_mesh/base/NgpField.hpp" +#include "stk_mesh/base/Types.hpp" + +namespace sierra { +namespace nalu { + +class Realm; + +class SDRKONodeKernel : public NGPNodeKernel +{ +public: + SDRKONodeKernel(const stk::mesh::MetaData&); + + SDRKONodeKernel() = delete; + + KOKKOS_DEFAULTED_FUNCTION + virtual ~SDRKONodeKernel() = default; + + virtual void setup(Realm&) override; + + KOKKOS_FUNCTION + virtual void execute( + NodeKernelTraits::LhsType&, + NodeKernelTraits::RhsType&, + const stk::mesh::FastMeshIndex&) override; + +private: + stk::mesh::NgpField tke_; + stk::mesh::NgpField sdr_; + stk::mesh::NgpField density_; + stk::mesh::NgpField tvisc_; + stk::mesh::NgpField visc_; + stk::mesh::NgpField dudx_; + stk::mesh::NgpField dkdx_; + stk::mesh::NgpField dwdx_; + stk::mesh::NgpField dualNodalVolume_; + + unsigned tkeID_{stk::mesh::InvalidOrdinal}; + unsigned sdrID_{stk::mesh::InvalidOrdinal}; + unsigned densityID_{stk::mesh::InvalidOrdinal}; + unsigned tviscID_{stk::mesh::InvalidOrdinal}; + unsigned viscID_{stk::mesh::InvalidOrdinal}; + unsigned dudxID_{stk::mesh::InvalidOrdinal}; + unsigned dkdxID_{stk::mesh::InvalidOrdinal}; + unsigned dwdxID_{stk::mesh::InvalidOrdinal}; + unsigned dualNodalVolumeID_{stk::mesh::InvalidOrdinal}; + + NodeKernelTraits::DblType betaStar_; + NodeKernelTraits::DblType tkeProdLimitRatio_; + NodeKernelTraits::DblType sigmaWTwo_; + NodeKernelTraits::DblType betaOne_; + NodeKernelTraits::DblType betaTwo_; + NodeKernelTraits::DblType gammaOne_; + NodeKernelTraits::DblType gammaTwo_; + + const int nDim_; +}; + +} // namespace nalu +} // namespace sierra + +#endif /* SDRKONODEKERNEL_H */ diff --git a/include/node_kernels/SDRSSTDESNodeKernel.h b/include/node_kernels/SDRSSTDESNodeKernel.h index 9911683f3d..5aa50bbd54 100644 --- a/include/node_kernels/SDRSSTDESNodeKernel.h +++ b/include/node_kernels/SDRSSTDESNodeKernel.h @@ -75,7 +75,6 @@ class SDRSSTDESNodeKernel : public NGPNodeKernel NodeKernelTraits::DblType betaTwo_; NodeKernelTraits::DblType gammaOne_; NodeKernelTraits::DblType gammaTwo_; - NodeKernelTraits::DblType relaxFac_; NodeKernelTraits::DblType cDESke_; NodeKernelTraits::DblType cDESkw_; NodeKernelTraits::DblType sdrAmb_; diff --git a/include/node_kernels/SDRSSTNodeKernel.h b/include/node_kernels/SDRSSTNodeKernel.h index 54008e7247..70cf0bd275 100644 --- a/include/node_kernels/SDRSSTNodeKernel.h +++ b/include/node_kernels/SDRSSTNodeKernel.h @@ -71,7 +71,6 @@ class SDRSSTNodeKernel : public NGPNodeKernel NodeKernelTraits::DblType gammaOne_; NodeKernelTraits::DblType gammaTwo_; NodeKernelTraits::DblType sdrAmb_; - NodeKernelTraits::DblType relaxFac_; const int nDim_; diff --git a/include/node_kernels/TDRKENodeKernel.h b/include/node_kernels/TDRKENodeKernel.h new file mode 100644 index 0000000000..56d0d03c76 --- /dev/null +++ b/include/node_kernels/TDRKENodeKernel.h @@ -0,0 +1,75 @@ +// 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 TDRKENODEKERNEL_H +#define TDRKENODEKERNEL_H + +#include "node_kernels/NodeKernel.h" +#include "FieldTypeDef.h" + +#include "stk_mesh/base/BulkData.hpp" +#include "stk_mesh/base/Ngp.hpp" +#include "stk_mesh/base/NgpField.hpp" +#include "stk_mesh/base/Types.hpp" + +namespace sierra { +namespace nalu { + +class Realm; + +class TDRKENodeKernel : public NGPNodeKernel +{ +public: + TDRKENodeKernel(const stk::mesh::MetaData&); + + TDRKENodeKernel() = delete; + + KOKKOS_DEFAULTED_FUNCTION + virtual ~TDRKENodeKernel() = default; + + virtual void setup(Realm&) override; + + KOKKOS_FUNCTION + virtual void execute( + NodeKernelTraits::LhsType&, + NodeKernelTraits::RhsType&, + const stk::mesh::FastMeshIndex&) override; + +private: + stk::mesh::NgpField tke_; + stk::mesh::NgpField tdr_; + stk::mesh::NgpField density_; + stk::mesh::NgpField tvisc_; + stk::mesh::NgpField visc_; + stk::mesh::NgpField wallDist_; + stk::mesh::NgpField dplus_; + stk::mesh::NgpField dudx_; + stk::mesh::NgpField dualNodalVolume_; + + unsigned tkeID_{stk::mesh::InvalidOrdinal}; + unsigned tdrID_{stk::mesh::InvalidOrdinal}; + unsigned densityID_{stk::mesh::InvalidOrdinal}; + unsigned tviscID_{stk::mesh::InvalidOrdinal}; + unsigned viscID_{stk::mesh::InvalidOrdinal}; + unsigned wallDistID_{stk::mesh::InvalidOrdinal}; + unsigned dplusID_{stk::mesh::InvalidOrdinal}; + unsigned dudxID_{stk::mesh::InvalidOrdinal}; + unsigned dualNodalVolumeID_{stk::mesh::InvalidOrdinal}; + + NodeKernelTraits::DblType cEpsOne_; + NodeKernelTraits::DblType cEpsTwo_; + NodeKernelTraits::DblType fOne_; + + const int nDim_; +}; + +} // namespace nalu +} // namespace sierra + +#endif /* TDRKENODEKERNEL_H */ diff --git a/include/node_kernels/TKEKENodeKernel.h b/include/node_kernels/TKEKENodeKernel.h new file mode 100644 index 0000000000..d81c636853 --- /dev/null +++ b/include/node_kernels/TKEKENodeKernel.h @@ -0,0 +1,71 @@ +// 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 TKEKENODEKERNEL_H +#define TKEKENODEKERNEL_H + +#include "node_kernels/NodeKernel.h" +#include "FieldTypeDef.h" +#include "stk_mesh/base/BulkData.hpp" +#include "stk_mesh/base/Ngp.hpp" +#include "stk_mesh/base/NgpField.hpp" +#include "stk_mesh/base/Types.hpp" + +namespace sierra { +namespace nalu { + +class Realm; + +class TKEKENodeKernel : public NGPNodeKernel +{ +public: + TKEKENodeKernel(const stk::mesh::MetaData&); + + TKEKENodeKernel() = delete; + + KOKKOS_DEFAULTED_FUNCTION + virtual ~TKEKENodeKernel() = default; + + virtual void setup(Realm&) override; + + KOKKOS_FUNCTION + virtual void execute( + NodeKernelTraits::LhsType&, + NodeKernelTraits::RhsType&, + const stk::mesh::FastMeshIndex&) override; + +private: + stk::mesh::NgpField tke_; + stk::mesh::NgpField tdr_; + stk::mesh::NgpField density_; + stk::mesh::NgpField visc_; + stk::mesh::NgpField tvisc_; + stk::mesh::NgpField dudx_; + stk::mesh::NgpField wallDist_; + stk::mesh::NgpField dualNodalVolume_; + + unsigned tkeID_{stk::mesh::InvalidOrdinal}; + unsigned tdrID_{stk::mesh::InvalidOrdinal}; + unsigned densityID_{stk::mesh::InvalidOrdinal}; + unsigned tviscID_{stk::mesh::InvalidOrdinal}; + unsigned viscID_{stk::mesh::InvalidOrdinal}; + unsigned dudxID_{stk::mesh::InvalidOrdinal}; + unsigned wallDistID_{stk::mesh::InvalidOrdinal}; + unsigned dualNodalVolumeID_{stk::mesh::InvalidOrdinal}; + + NodeKernelTraits::DblType betaStar_; + NodeKernelTraits::DblType tkeProdLimitRatio_; + + const int nDim_; +}; + +} // namespace nalu +} // namespace sierra + +#endif /* TKEKENODEKERNEL_H */ diff --git a/include/node_kernels/TKEKONodeKernel.h b/include/node_kernels/TKEKONodeKernel.h new file mode 100644 index 0000000000..44fe9c3888 --- /dev/null +++ b/include/node_kernels/TKEKONodeKernel.h @@ -0,0 +1,69 @@ +// 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 TKEKONODEKERNEL_H +#define TKEKONODEKERNEL_H + +#include "node_kernels/NodeKernel.h" +#include "FieldTypeDef.h" +#include "stk_mesh/base/BulkData.hpp" +#include "stk_mesh/base/Ngp.hpp" +#include "stk_mesh/base/NgpField.hpp" +#include "stk_mesh/base/Types.hpp" + +namespace sierra { +namespace nalu { + +class Realm; + +class TKEKONodeKernel : public NGPNodeKernel +{ +public: + TKEKONodeKernel(const stk::mesh::MetaData&); + + TKEKONodeKernel() = delete; + + KOKKOS_DEFAULTED_FUNCTION + virtual ~TKEKONodeKernel() = default; + + virtual void setup(Realm&) override; + + KOKKOS_FUNCTION + virtual void execute( + NodeKernelTraits::LhsType&, + NodeKernelTraits::RhsType&, + const stk::mesh::FastMeshIndex&) override; + +private: + stk::mesh::NgpField tke_; + stk::mesh::NgpField sdr_; + stk::mesh::NgpField density_; + stk::mesh::NgpField tvisc_; + stk::mesh::NgpField visc_; + stk::mesh::NgpField dudx_; + stk::mesh::NgpField dualNodalVolume_; + + unsigned tkeID_{stk::mesh::InvalidOrdinal}; + unsigned sdrID_{stk::mesh::InvalidOrdinal}; + unsigned densityID_{stk::mesh::InvalidOrdinal}; + unsigned tviscID_{stk::mesh::InvalidOrdinal}; + unsigned viscID_{stk::mesh::InvalidOrdinal}; + unsigned dudxID_{stk::mesh::InvalidOrdinal}; + unsigned dualNodalVolumeID_{stk::mesh::InvalidOrdinal}; + + NodeKernelTraits::DblType betaStar_; + NodeKernelTraits::DblType tkeProdLimitRatio_; + + const int nDim_; +}; + +} // namespace nalu +} // namespace sierra + +#endif /* TKEKONODEKERNEL_H */ diff --git a/include/node_kernels/TKEKsgsNodeKernel.h b/include/node_kernels/TKEKsgsNodeKernel.h index bd0dd0b55c..2e8fd8a8cd 100644 --- a/include/node_kernels/TKEKsgsNodeKernel.h +++ b/include/node_kernels/TKEKsgsNodeKernel.h @@ -59,7 +59,6 @@ class TKEKsgsNodeKernel : public NGPNodeKernel NodeKernelTraits::DblType cEps_; NodeKernelTraits::DblType tkeProdLimitRatio_; - NodeKernelTraits::DblType relaxFac_; const int nDim_; }; diff --git a/include/node_kernels/TKESSTDESNodeKernel.h b/include/node_kernels/TKESSTDESNodeKernel.h index fc6360a4eb..f4c1d0ec82 100644 --- a/include/node_kernels/TKESSTDESNodeKernel.h +++ b/include/node_kernels/TKESSTDESNodeKernel.h @@ -67,7 +67,6 @@ class TKESSTDESNodeKernel : public NGPNodeKernel NodeKernelTraits::DblType cDESkw_; NodeKernelTraits::DblType tkeAmb_; NodeKernelTraits::DblType sdrAmb_; - NodeKernelTraits::DblType relaxFac_; const int nDim_; }; diff --git a/include/node_kernels/TKESSTIDDESNodeKernel.h b/include/node_kernels/TKESSTIDDESNodeKernel.h index 8a5a520026..b8d226931d 100644 --- a/include/node_kernels/TKESSTIDDESNodeKernel.h +++ b/include/node_kernels/TKESSTIDDESNodeKernel.h @@ -74,7 +74,6 @@ class TKESSTIDDESNodeKernel : public NGPNodeKernel NodeKernelTraits::DblType iddes_Ct_; NodeKernelTraits::DblType tkeAmb_; NodeKernelTraits::DblType sdrAmb_; - NodeKernelTraits::DblType relaxFac_; const int nDim_; }; diff --git a/include/node_kernels/TKESSTNodeKernel.h b/include/node_kernels/TKESSTNodeKernel.h index 4c1a228b8c..c85f20e04f 100644 --- a/include/node_kernels/TKESSTNodeKernel.h +++ b/include/node_kernels/TKESSTNodeKernel.h @@ -59,7 +59,6 @@ class TKESSTNodeKernel : public NGPNodeKernel NodeKernelTraits::DblType tkeProdLimitRatio_; NodeKernelTraits::DblType tkeAmb_; NodeKernelTraits::DblType sdrAmb_; - NodeKernelTraits::DblType relaxFac_; const int nDim_; }; diff --git a/reg_tests/CTestList.cmake b/reg_tests/CTestList.cmake index 6679ae4f98..0c55840906 100644 --- a/reg_tests/CTestList.cmake +++ b/reg_tests/CTestList.cmake @@ -199,6 +199,8 @@ if(NOT ENABLE_CUDA) add_test_r(ablNeutralEdgeAMS 4) add_test_r_rst(amsChannelEdge 4) add_test_r(SSTChannelEdge 4) + add_test_r(KOChannelEdge 4) + add_test_r(KEChannelEdge 4) add_test_r_rst(SSTAMSChannelEdge 4) #add_test_r_rst(SSTAMSOversetRotCylinder 2) add_test_r(ablNeutralNGPHypre 2) diff --git a/reg_tests/test_files/KEChannelEdge/KEChannelEdge.yaml b/reg_tests/test_files/KEChannelEdge/KEChannelEdge.yaml new file mode 100644 index 0000000000..afdf67534d --- /dev/null +++ b/reg_tests/test_files/KEChannelEdge/KEChannelEdge.yaml @@ -0,0 +1,200 @@ +Simulations: + - name: sim1 + time_integrator: ti_1 + optimizer: opt1 + +linear_solvers: + + - name: solve_scalar + type: tpetra + method: gmres + preconditioner: riluk + tolerance: 1e-5 + max_iterations: 200 + kspace: 200 + output_level: 0 + + - name: solve_cont + type: hypre + method: hypre_gmres + preconditioner: boomerAMG + tolerance: 1e-5 + max_iterations: 50 + kspace: 75 + output_level: 0 + bamg_coarsen_type: 8 + bamg_interp_type: 6 + bamg_cycle_type: 1 + +realms: + + - name: realm_1 + mesh: ../../mesh/periodicChannel.g + use_edges: yes + check_for_missing_bcs: yes + automatic_decomposition_type: rcb + + time_step_control: + target_courant: 1.0 + time_step_change_factor: 1.2 + + equation_systems: + name: theEqSys + max_iterations: 4 + + solver_system_specification: + velocity: solve_scalar + turbulent_ke: solve_scalar + total_dissipation_rate: solve_scalar + pressure: solve_cont + ndtw: solve_cont + + systems: + - WallDistance: + name: myNDTW + max_iterations: 1 + convergence_tolerance: 1.0e-8 + + - LowMachEOM: + name: myLowMach + max_iterations: 1 + convergence_tolerance: 1e-8 + + - ChienKEpsilon: + name: myKE + max_iterations: 1 + convergence_tolerance: 1e-8 + + initial_conditions: + - constant: ic_1 + target_name: Unspecified-2-HEX + value: + pressure: 0 + velocity: [22.0,0.0,0.0] + turbulent_ke: 0.0005 + total_dissipation_rate: 0.0005 + + material_properties: + target_name: Unspecified-2-HEX + specifications: + - name: density + type: constant + value: 1.0 + - name: viscosity + type: constant + value: 9.99488e-4 + + boundary_conditions: + + - wall_boundary_condition: bc_bot + target_name: bottom + wall_user_data: + velocity: [0,0,0] + turbulent_ke: 0.0 + use_wall_function: no + + - wall_boundary_condition: bc_top + target_name: top + wall_user_data: + velocity: [0,0,0] + turbulent_ke: 0.0 + use_wall_function: no + + - periodic_boundary_condition: bc_inlet_outlet + target_name: [inlet, outlet] + periodic_user_data: + search_tolerance: 0.001 + + - periodic_boundary_condition: bc_front_back + target_name: [front, back] + periodic_user_data: + search_tolerance: 0.001 + + solution_options: + name: myOptions + turbulence_model: ke + projected_timescale_type: momentum_diag_inv + + fix_pressure_at_node: + value: 0.0 + node_lookup_type: spatial_location + location: [ 1.0, 1.0, 1.0 ] + search_target_part: [Unspecified-2-HEX] + search_method: stk_kdtree + + options: + - hybrid_factor: + velocity: 1.0 + turbulent_ke: 1.0 + total_dissipation_rate: 1.0 + + - alpha_upw: + velocity: 1.0 + turbulent_ke: 1.0 + total_dissipation_rate: 1.0 + + - upw_factor: + velocity: 1.0 + turbulent_ke: 0.0 + total_dissipation_rate: 0.0 + + - noc_correction: + pressure: yes + + - limiter: + pressure: no + velocity: yes + turbulent_ke: yes + total_dissipation_rate: yes + + - projected_nodal_gradient: + velocity: element + pressure: element + turbulent_ke: element + total_dissipation_rate: element + ndtw: element + + - relaxation_factor: + velocity: 1.0 + pressure: 1.0 + turbulent_ke: 1.0 + total_dissipation_rate: 1.0 + + - source_terms: + momentum: body_force + + - source_term_parameters: + momentum: [1.00, 0.0, 0.0] + + restart: + restart_data_base_name: KEChannelEdge.rst + restart_frequency: 10 + restart_start: 0 + + output: + output_data_base_name: KEChannelEdge.e + output_frequency: 10 + output_node_set: no + output_variables: + - velocity + - pressure + - pressure_force + - tau_wall + - turbulent_ke + - total_dissipation_rate + - minimum_distance_to_wall + - turbulent_viscosity + - element_courant + +Time_Integrators: + - StandardTimeIntegrator: + name: ti_1 + start_time: 0 + time_step: 2.0e-3 + termination_step_count: 20 + time_stepping_type: fixed + time_step_count: 0 + second_order_accuracy: yes + + realms: + - realm_1 diff --git a/reg_tests/test_files/KOChannelEdge/KOChannelEdge.yaml b/reg_tests/test_files/KOChannelEdge/KOChannelEdge.yaml new file mode 100644 index 0000000000..af629dc3f1 --- /dev/null +++ b/reg_tests/test_files/KOChannelEdge/KOChannelEdge.yaml @@ -0,0 +1,203 @@ +Simulations: + - name: sim1 + time_integrator: ti_1 + optimizer: opt1 + +linear_solvers: + + - name: solve_scalar + type: tpetra + method: gmres + preconditioner: riluk + tolerance: 1e-5 + max_iterations: 200 + kspace: 200 + output_level: 0 + + - name: solve_cont + type: hypre + method: hypre_gmres + preconditioner: boomerAMG + tolerance: 1e-5 + max_iterations: 50 + kspace: 75 + output_level: 0 + bamg_coarsen_type: 8 + bamg_interp_type: 6 + bamg_cycle_type: 1 + +realms: + + - name: realm_1 + mesh: ../../mesh/periodicChannel.g + use_edges: yes + check_for_missing_bcs: yes + automatic_decomposition_type: rcb + + time_step_control: + target_courant: 1.0 + time_step_change_factor: 1.2 + + equation_systems: + name: theEqSys + max_iterations: 4 + + solver_system_specification: + velocity: solve_scalar + turbulent_ke: solve_scalar + specific_dissipation_rate: solve_scalar + pressure: solve_cont + ndtw: solve_cont + + systems: + - WallDistance: + name: myNDTW + max_iterations: 1 + convergence_tolerance: 1.0e-8 + + - LowMachEOM: + name: myLowMach + max_iterations: 1 + convergence_tolerance: 1e-8 + + - WilcoxKOmega: + name: myKO + max_iterations: 1 + convergence_tolerance: 1e-8 + + initial_conditions: + - constant: ic_1 + target_name: Unspecified-2-HEX + value: + pressure: 0 + velocity: [22.0,0.0,0.0] + turbulent_ke: 0.05 + specific_dissipation_rate: 3528.0 + + material_properties: + target_name: Unspecified-2-HEX + specifications: + - name: density + type: constant + value: 1.0 + - name: viscosity + type: constant + value: 9.99488e-4 + + boundary_conditions: + + - wall_boundary_condition: bc_bot + target_name: bottom + wall_user_data: + velocity: [0,0,0] + turbulent_ke: 0.0 + use_wall_function: no + + - wall_boundary_condition: bc_top + target_name: top + wall_user_data: + velocity: [0,0,0] + turbulent_ke: 0.0 + use_wall_function: no + + - periodic_boundary_condition: bc_inlet_outlet + target_name: [inlet, outlet] + periodic_user_data: + search_tolerance: 0.001 + + - periodic_boundary_condition: bc_front_back + target_name: [front, back] + periodic_user_data: + search_tolerance: 0.001 + + solution_options: + name: myOptions + turbulence_model: ko + projected_timescale_type: momentum_diag_inv + + fix_pressure_at_node: + value: 0.0 + node_lookup_type: spatial_location + location: [ 1.0, 1.0, 1.0 ] + search_target_part: [Unspecified-2-HEX] + search_method: stk_kdtree + + options: + - hybrid_factor: + velocity: 1.0 + turbulent_ke: 1.0 + specific_dissipation_rate: 1.0 + + - alpha_upw: + velocity: 1.0 + turbulent_ke: 1.0 + specific_dissipation_rate: 1.0 + + - upw_factor: + velocity: 1.0 + turbulent_ke: 0.0 + specific_dissipation_rate: 0.0 + + - noc_correction: + pressure: yes + + - limiter: + pressure: no + velocity: yes + turbulent_ke: yes + specific_dissipation_rate: yes + + - projected_nodal_gradient: + velocity: element + pressure: element + turbulent_ke: element + specific_dissipation_rate: element + ndtw: element + + - relaxation_factor: + velocity: 1.0 + pressure: 1.0 + turbulent_ke: 1.0 + specific_dissipation_rate: 1.0 + + - turbulence_model_constants: + SDRWallFactor: 0.625 + + - source_terms: + momentum: body_force + + - source_term_parameters: + momentum: [1.00, 0.0, 0.0] + + restart: + restart_data_base_name: KOChannelEdge.rst + restart_frequency: 10 + restart_start: 0 + + output: + output_data_base_name: KOChannelEdge.e + output_frequency: 10 + output_node_set: no + output_variables: + - velocity + - pressure + - pressure_force + - tau_wall + - turbulent_ke + - specific_dissipation_rate + - minimum_distance_to_wall + - turbulent_viscosity + - element_courant + +Time_Integrators: + - StandardTimeIntegrator: + name: ti_1 + start_time: 0 + time_step: 2.0e-3 + termination_step_count: 20 + time_stepping_type: fixed + time_step_count: 0 + second_order_accuracy: yes + + realms: + - realm_1 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index e0e9c37b53..66c4c16264 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -23,6 +23,7 @@ target_sources(nalu PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/AuxFunctionAlgorithm.C ${CMAKE_CURRENT_SOURCE_DIR}/AveragingInfo.C ${CMAKE_CURRENT_SOURCE_DIR}/BoundaryConditions.C + ${CMAKE_CURRENT_SOURCE_DIR}/ChienKEpsilonEquationSystem.C ${CMAKE_CURRENT_SOURCE_DIR}/ComputeHeatTransferEdgeWallAlgorithm.C ${CMAKE_CURRENT_SOURCE_DIR}/ComputeMdotNonConformalAlgorithm.C ${CMAKE_CURRENT_SOURCE_DIR}/ComputeSSTMaxLengthScaleElemAlgorithm.C @@ -86,6 +87,7 @@ target_sources(nalu PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/SurfaceForceAndMomentAlgorithmDriver.C ${CMAKE_CURRENT_SOURCE_DIR}/SurfaceForceAndMomentWallFunctionAlgorithm.C ${CMAKE_CURRENT_SOURCE_DIR}/TimeIntegrator.C + ${CMAKE_CURRENT_SOURCE_DIR}/TotalDissipationRateEquationSystem.C ${CMAKE_CURRENT_SOURCE_DIR}/TpetraLinearSystem.C ${CMAKE_CURRENT_SOURCE_DIR}/TpetraLinearSystemHelpers.C ${CMAKE_CURRENT_SOURCE_DIR}/TpetraSegregatedLinearSystem.C @@ -94,6 +96,7 @@ target_sources(nalu PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/TurbViscWaleAlgorithm.C ${CMAKE_CURRENT_SOURCE_DIR}/TurbulenceAveragingPostProcessing.C ${CMAKE_CURRENT_SOURCE_DIR}/WallDistEquationSystem.C + ${CMAKE_CURRENT_SOURCE_DIR}/WilcoxKOmegaEquationSystem.C ) if(ENABLE_FFTW) diff --git a/src/ChienKEpsilonEquationSystem.C b/src/ChienKEpsilonEquationSystem.C new file mode 100644 index 0000000000..6344769d78 --- /dev/null +++ b/src/ChienKEpsilonEquationSystem.C @@ -0,0 +1,437 @@ +// 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 +#include +#include +#include +#include +#include +#include +#include + +// ngp +#include "FieldTypeDef.h" +#include "ngp_algorithms/GeometryAlgDriver.h" +#include "ngp_algorithms/WallFuncGeometryAlg.h" +#include "ngp_utils/NgpLoopUtils.h" +#include "ngp_utils/NgpFieldUtils.h" +#include "ngp_utils/NgpFieldBLAS.h" + +// stk_util +#include +#include "utils/StkHelpers.h" + +// stk_mesh/base/fem +#include +#include +#include + +#include + +// stk_io +#include + +// basic c++ +#include +#include +#include + +// ngp +#include "ngp_utils/NgpFieldBLAS.h" + +namespace sierra { +namespace nalu { + +//========================================================================== +// Class Definition +//========================================================================== +// ChienKEpsilonEquationSystem - manage SST +//========================================================================== +//-------------------------------------------------------------------------- +//-------- constructor ----------------------------------------------------- +//-------------------------------------------------------------------------- +ChienKEpsilonEquationSystem::ChienKEpsilonEquationSystem( + EquationSystems& eqSystems) + : EquationSystem(eqSystems, "ChienKEpsilonWrap"), + tkeEqSys_(std::make_unique(eqSystems)), + tdrEqSys_(std::make_unique(eqSystems)), + tke_(NULL), + tdr_(NULL), + minDistanceToWall_(NULL), + dplus_(NULL), + isInit_(true), + resetAMSAverages_(realm_.solutionOptions_->resetAMSAverages_) +{ + // push back EQ to manager + realm_.push_equation_to_systems(this); +} + +//-------------------------------------------------------------------------- +//-------- destructor ------------------------------------------------------ +//-------------------------------------------------------------------------- +ChienKEpsilonEquationSystem::~ChienKEpsilonEquationSystem() {} + +void +ChienKEpsilonEquationSystem::load(const YAML::Node& node) +{ + EquationSystem::load(node); + + if (realm_.query_for_overset()) { + tkeEqSys_->decoupledOverset_ = decoupledOverset_; + tkeEqSys_->numOversetIters_ = numOversetIters_; + tdrEqSys_->decoupledOverset_ = decoupledOverset_; + tdrEqSys_->numOversetIters_ = numOversetIters_; + } +} + +//-------------------------------------------------------------------------- +//-------- initialize ------------------------------------------------------ +//-------------------------------------------------------------------------- +void +ChienKEpsilonEquationSystem::initialize() +{ + // let equation systems that are owned some information + tkeEqSys_->convergenceTolerance_ = convergenceTolerance_; + tdrEqSys_->convergenceTolerance_ = convergenceTolerance_; +} + +//-------------------------------------------------------------------------- +//-------- register_nodal_fields ------------------------------------------- +//-------------------------------------------------------------------------- +void +ChienKEpsilonEquationSystem::register_nodal_fields(stk::mesh::Part* part) +{ + + stk::mesh::MetaData& meta_data = realm_.meta_data(); + const int numStates = realm_.number_of_states(); + + // re-register tke and tdr for convenience + tke_ = &(meta_data.declare_field( + stk::topology::NODE_RANK, "turbulent_ke", numStates)); + stk::mesh::put_field_on_mesh(*tke_, *part, nullptr); + tdr_ = &(meta_data.declare_field( + stk::topology::NODE_RANK, "total_dissipation_rate", numStates)); + stk::mesh::put_field_on_mesh(*tdr_, *part, nullptr); + + // SST parameters that everyone needs + minDistanceToWall_ = &(meta_data.declare_field( + stk::topology::NODE_RANK, "minimum_distance_to_wall")); + stk::mesh::put_field_on_mesh(*minDistanceToWall_, *part, nullptr); + dplus_ = &(meta_data.declare_field( + stk::topology::NODE_RANK, "dplus_wall_function")); + stk::mesh::put_field_on_mesh(*dplus_, *part, nullptr); + + // add to restart field + realm_.augment_restart_variable_list("minimum_distance_to_wall"); + realm_.augment_restart_variable_list("dplus_wall_function"); +} + +//-------------------------------------------------------------------------- +//-------- register_interior_algorithm ------------------------------------- +//-------------------------------------------------------------------------- +void +ChienKEpsilonEquationSystem::register_interior_algorithm( + stk::mesh::Part* /*part*/) +{ + // nothing to do here... +} + +//-------------------------------------------------------------------------- +//-------- register_wall_bc ------------------------------------------------ +//-------------------------------------------------------------------------- +void +ChienKEpsilonEquationSystem::register_wall_bc( + stk::mesh::Part* part, + const stk::topology& partTopo, + const WallBoundaryConditionData& wallBCData) +{ + // determine if using RANS for ABL + WallUserData userData = wallBCData.userData_; + bool RANSAblBcApproach = userData.RANSAblBcApproach_; + + // push mesh part + wallBcPart_.push_back(part); + + auto& meta = realm_.meta_data(); + auto& assembledWallArea = meta.declare_field( + stk::topology::NODE_RANK, "assembled_wall_area_wf"); + stk::mesh::put_field_on_mesh(assembledWallArea, *part, nullptr); + auto& assembledWallNormDist = meta.declare_field( + stk::topology::NODE_RANK, "assembled_wall_normal_distance"); + stk::mesh::put_field_on_mesh(assembledWallNormDist, *part, nullptr); + auto& wallNormDistBip = meta.declare_field( + meta.side_rank(), "wall_normal_distance_bip"); + auto* meFC = MasterElementRepo::get_surface_master_element(partTopo); + const int numScsBip = meFC->num_integration_points(); + stk::mesh::put_field_on_mesh(wallNormDistBip, *part, numScsBip, nullptr); +} + +//-------------------------------------------------------------------------- +//-------- solve_and_update ------------------------------------------------ +//-------------------------------------------------------------------------- +void +ChienKEpsilonEquationSystem::solve_and_update() +{ + // wrap timing + // SST_FIXME: deal with timers; all on misc for SSTEqs double timeA, timeB; + if (isInit_) { + // compute projected nodal gradients + tkeEqSys_->compute_projected_nodal_gradient(); + tdrEqSys_->assemble_nodal_gradient(); + clip_min_distance_to_wall(); + compute_dplus_function(); + + isInit_ = false; + } else if (realm_.has_mesh_motion()) { + if (realm_.currentNonlinearIteration_ == 1) + clip_min_distance_to_wall(); + compute_dplus_function(); + } + + // SST effective viscosity for k and omega + tkeEqSys_->compute_effective_diff_flux_coeff(); + tdrEqSys_->compute_effective_diff_flux_coeff(); + + // start the iteration loop + for (int k = 0; k < maxIterations_; ++k) { + + NaluEnv::self().naluOutputP0() + << " " << k + 1 << "/" << maxIterations_ << std::setw(15) << std::right + << name_ << std::endl; + + for (int oi = 0; oi < numOversetIters_; ++oi) { + // tke and tdr assemble, load_complete and solve; Jacobi iteration + tkeEqSys_->assemble_and_solve(tkeEqSys_->kTmp_); + tdrEqSys_->assemble_and_solve(tdrEqSys_->eTmp_); + + update_and_clip(); + + if (decoupledOverset_ && realm_.hasOverset_) { + realm_.overset_field_update(tkeEqSys_->tke_, 1, 1); + realm_.overset_field_update(tdrEqSys_->tdr_, 1, 1); + } + } + // compute projected nodal gradients + tkeEqSys_->compute_projected_nodal_gradient(); + tdrEqSys_->assemble_nodal_gradient(); + } +} + +/** Perform sanity checks on TKE/TDR fields + */ +void +ChienKEpsilonEquationSystem::initial_work() +{ + const auto& meshInfo = realm_.mesh_info(); + const auto& meta = meshInfo.meta(); + const auto& ngpMesh = meshInfo.ngp_mesh(); + const auto& fieldMgr = meshInfo.ngp_field_manager(); + + auto& tkeNp1 = fieldMgr.get_field(tke_->mesh_meta_data_ordinal()); + auto& tdrNp1 = fieldMgr.get_field(tdr_->mesh_meta_data_ordinal()); + const stk::mesh::Selector sel = + (meta.locally_owned_part() | meta.globally_shared_part()) & + stk::mesh::selectField(*tdr_); + clip_ke(ngpMesh, sel, tkeNp1, tdrNp1); +} + +void +ChienKEpsilonEquationSystem::post_external_data_transfer_work() +{ + const auto& meshInfo = realm_.mesh_info(); + const auto& meta = meshInfo.meta(); + const auto& ngpMesh = meshInfo.ngp_mesh(); + const auto& fieldMgr = meshInfo.ngp_field_manager(); + + auto& tkeNp1 = fieldMgr.get_field(tke_->mesh_meta_data_ordinal()); + auto& tdrNp1 = fieldMgr.get_field(tdr_->mesh_meta_data_ordinal()); + + const stk::mesh::Selector owned_and_shared = + (meta.locally_owned_part() | meta.globally_shared_part()); + auto interior_sel = owned_and_shared & stk::mesh::selectField(*tdr_); + clip_ke(ngpMesh, interior_sel, tkeNp1, tdrNp1); + + auto tdrBCField = + meta.get_field(stk::topology::NODE_RANK, "tdr_bc"); + auto tkeBCField = + meta.get_field(stk::topology::NODE_RANK, "tke_bc"); + if (tdrBCField != nullptr) { + ThrowRequire(tkeBCField); + auto bc_sel = owned_and_shared & stk::mesh::selectField(*tdrBCField); + auto ngpTkeBC = + fieldMgr.get_field(tkeBCField->mesh_meta_data_ordinal()); + auto ngpTdrBC = + fieldMgr.get_field(tdrBCField->mesh_meta_data_ordinal()); + clip_ke(ngpMesh, bc_sel, ngpTkeBC, ngpTdrBC); + } +} + +/** Update solution but ensure that TKE and TDR are greater than zero + */ +void +ChienKEpsilonEquationSystem::update_and_clip() +{ + using MeshIndex = nalu_ngp::NGPMeshTraits<>::MeshIndex; + + const auto& meshInfo = realm_.mesh_info(); + const auto& meta = meshInfo.meta(); + const auto& ngpMesh = meshInfo.ngp_mesh(); + const auto& fieldMgr = meshInfo.ngp_field_manager(); + + auto& tkeNp1 = fieldMgr.get_field(tke_->mesh_meta_data_ordinal()); + auto& tdrNp1 = fieldMgr.get_field(tdr_->mesh_meta_data_ordinal()); + const auto& kTmp = + fieldMgr.get_field(tkeEqSys_->kTmp_->mesh_meta_data_ordinal()); + const auto& eTmp = + fieldMgr.get_field(tdrEqSys_->eTmp_->mesh_meta_data_ordinal()); + + auto* turbViscosity = meta.get_field( + stk::topology::NODE_RANK, "turbulent_viscosity"); + + const stk::mesh::Selector sel = + (meta.locally_owned_part() | meta.globally_shared_part()) & + stk::mesh::selectField(*turbViscosity); + + // Bring class variables to local scope for lambda capture + const double tkeMinVal = tkeMinValue_; + const double tdrMinVal = tdrMinValue_; + + tkeNp1.sync_to_device(); + tdrNp1.sync_to_device(); + + nalu_ngp::run_entity_algorithm( + "KE::update_and_clip", ngpMesh, stk::topology::NODE_RANK, sel, + KOKKOS_LAMBDA(const MeshIndex& mi) { + const double tkeNew = tkeNp1.get(mi, 0) + kTmp.get(mi, 0); + const double tdrNew = tdrNp1.get(mi, 0) + eTmp.get(mi, 0); + + tkeNp1.get(mi, 0) = (tkeNew < 0.0) ? tkeMinVal : tkeNew; + tdrNp1.get(mi, 0) = stk::math::max(tdrNew, tdrMinVal); + }); + + tkeNp1.modify_on_device(); + tdrNp1.modify_on_device(); +} + +void +ChienKEpsilonEquationSystem::clip_ke( + const stk::mesh::NgpMesh& ngpMesh, + const stk::mesh::Selector& sel, + stk::mesh::NgpField& tke, + stk::mesh::NgpField& tdr) +{ + tke.sync_to_device(); + tdr.sync_to_device(); + + // Bring class variables to local scope for lambda capture + const double tkeMinVal = tkeMinValue_; + const double tdrMinVal = tdrMinValue_; + + nalu_ngp::run_entity_algorithm( + "KE::clip", ngpMesh, stk::topology::NODE_RANK, sel, + KOKKOS_LAMBDA(const nalu_ngp::NGPMeshTraits<>::MeshIndex& mi) { + const double tkeNew = tke.get(mi, 0); + const double tdrNew = tdr.get(mi, 0); + + tke.get(mi, 0) = (tkeNew < 0.0) ? tkeMinVal : tkeNew; + tdr.get(mi, 0) = stk::math::max(tdrNew, tdrMinVal); + }); + tke.modify_on_device(); + tdr.modify_on_device(); +} + +//-------------------------------------------------------------------------- +//-------- clip_min_distance_to_wall --------------------------------------- +//-------------------------------------------------------------------------- +void +ChienKEpsilonEquationSystem::clip_min_distance_to_wall() +{ + using MeshIndex = nalu_ngp::NGPMeshTraits<>::MeshIndex; + const auto& meshInfo = realm_.mesh_info(); + const auto& ngpMesh = meshInfo.ngp_mesh(); + const auto& meta = meshInfo.meta(); + const auto& fieldMgr = meshInfo.ngp_field_manager(); + + auto& ndtw = + fieldMgr.get_field(minDistanceToWall_->mesh_meta_data_ordinal()); + const auto& wallNormDist = + nalu_ngp::get_ngp_field(meshInfo, "assembled_wall_normal_distance"); + + const stk::mesh::Selector sel = + (meta.locally_owned_part() | meta.globally_shared_part()) & + stk::mesh::selectUnion(wallBcPart_); + + ndtw.sync_to_device(); + + nalu_ngp::run_entity_algorithm( + "SST::clip_ndtw", ngpMesh, stk::topology::NODE_RANK, sel, + KOKKOS_LAMBDA(const MeshIndex& mi) { + const double minD = ndtw.get(mi, 0); + + ndtw.get(mi, 0) = stk::math::max(minD, wallNormDist.get(mi, 0)); + }); + ndtw.modify_on_device(); + + stk::mesh::parallel_max(realm_.bulk_data(), {minDistanceToWall_}); + if (realm_.hasPeriodic_) { + realm_.periodic_field_max(minDistanceToWall_, 1); + } +} + +/** Compute non-local function of distance to wall + */ +void +ChienKEpsilonEquationSystem::compute_dplus_function() +{ + using MeshIndex = nalu_ngp::NGPMeshTraits<>::MeshIndex; + + const auto& meshInfo = realm_.mesh_info(); + const auto& meta = meshInfo.meta(); + const auto& ngpMesh = meshInfo.ngp_mesh(); + const auto& fieldMgr = meshInfo.ngp_field_manager(); + + const double utau = realm_.get_turb_model_constant(TM_utau); + + const auto& density = nalu_ngp::get_ngp_field(meshInfo, "density"); + const auto& viscosity = nalu_ngp::get_ngp_field(meshInfo, "viscosity"); + const auto& ndtw = + fieldMgr.get_field(minDistanceToWall_->mesh_meta_data_ordinal()); + auto& dplus = fieldMgr.get_field(dplus_->mesh_meta_data_ordinal()); + + const stk::mesh::Selector sel = + (meta.locally_owned_part() | meta.globally_shared_part()) & + (stk::mesh::selectField(*dplus_)); + + dplus.sync_to_device(); + + nalu_ngp::run_entity_algorithm( + "KE::compute_dplus_function", ngpMesh, stk::topology::NODE_RANK, sel, + KOKKOS_LAMBDA(const MeshIndex& mi) { + const double rho = density.get(mi, 0); + const double mu = viscosity.get(mi, 0); + const double minD = ndtw.get(mi, 0); + + dplus.get(mi, 0) = minD * rho * utau / mu; + }); + + dplus.modify_on_device(); +} + +void +ChienKEpsilonEquationSystem::post_iter_work() +{ + // nothing to do here ... +} + +} // namespace nalu +} // namespace sierra diff --git a/src/EquationSystems.C b/src/EquationSystems.C index 36de374c33..4caae8bac4 100644 --- a/src/EquationSystems.C +++ b/src/EquationSystems.C @@ -27,6 +27,8 @@ #include #include #include +#include +#include #include #include @@ -124,6 +126,16 @@ void EquationSystems::load(const YAML::Node & y_node) if (root()->debug()) NaluEnv::self().naluOutputP0() << "eqSys = tke/sdr " << std::endl; eqSys = new ShearStressTransportEquationSystem(*this); + } else if (expect_map(y_system, "ChienKEpsilon", true)) { + y_eqsys = expect_map(y_system, "ChienKEpsilon", true); + if (root()->debug()) + NaluEnv::self().naluOutputP0() << "eqSys = tke/tdr " << std::endl; + eqSys = new ChienKEpsilonEquationSystem(*this); + } else if (expect_map(y_system, "WilcoxKOmega", true)) { + y_eqsys = expect_map(y_system, "WilcoxKOmega", true); + if (root()->debug()) + NaluEnv::self().naluOutputP0() << "eqSys = tke/sdr " << std::endl; + eqSys = new WilcoxKOmegaEquationSystem(*this); } else if (expect_map(y_system, "TurbKineticEnergy", true)) { y_eqsys = expect_map(y_system, "TurbKineticEnergy", true); if (root()->debug()) diff --git a/src/LowMachEquationSystem.C b/src/LowMachEquationSystem.C index a2a95d5096..abdbae4b15 100644 --- a/src/LowMachEquationSystem.C +++ b/src/LowMachEquationSystem.C @@ -116,6 +116,8 @@ #include "ngp_algorithms/EffDiffFluxCoeffAlg.h" #include "ngp_algorithms/TurbViscKsgsAlg.h" #include "ngp_algorithms/TurbViscSSTAlg.h" +#include "ngp_algorithms/TurbViscKEAlg.h" +#include "ngp_algorithms/TurbViscKOAlg.h" #include "ngp_algorithms/WallFuncGeometryAlg.h" #include "ngp_algorithms/DynamicPressureOpenAlg.h" #include "ngp_algorithms/MomentumABLWallFuncMaskUtil.h" @@ -1490,6 +1492,14 @@ MomentumEquationSystem::register_interior_algorithm( tviscAlg_.reset(new TurbViscSSTAlg(realm_, part, tvisc_, true)); break; + case TurbulenceModel::KE: + tviscAlg_.reset(new TurbViscKEAlg(realm_, part, tvisc_)); + break; + + case TurbulenceModel::KO: + tviscAlg_.reset(new TurbViscKOAlg(realm_, part, tvisc_)); + break; + default: throw std::runtime_error("Unsupported turbulence model provided"); } diff --git a/src/NaluParsing.C b/src/NaluParsing.C index 836f39c3f4..b61078b305 100644 --- a/src/NaluParsing.C +++ b/src/NaluParsing.C @@ -622,17 +622,30 @@ namespace YAML return true; } - bool convert::decode(const Node& node, - sierra::nalu::GammaInf& gamma) - { - if (!node.IsScalar()) - { - return false; - } + bool + convert::decode( + const Node& node, sierra::nalu::TotDissRate& tdr) + { + if (!node.IsScalar()) { + return false; + } + + tdr.totDissRate_ = node.as(); - gamma.gamma_ = node.as(); + return true; + } - return true; + bool + convert::decode( + const Node& node, sierra::nalu::GammaInf& gamma) + { + if (!node.IsScalar()) { + return false; + } + + gamma.gamma_ = node.as(); + + return true; } bool convert::decode(const Node& node, @@ -755,6 +768,13 @@ namespace YAML wallData.bcDataSpecifiedMap_["turbulent_ke"] = true; wallData.bcDataTypeMap_["turbulent_ke"] = sierra::nalu::CONSTANT_UD; } + if (node["total_dissipation_rate"]) { + wallData.tdr_ = + node["total_dissipation_rate"].as(); + wallData.bcDataSpecifiedMap_["total_dissipation_rate"] = true; + wallData.bcDataTypeMap_["total_dissipation_rate"] = + sierra::nalu::CONSTANT_UD; + } if (node["temperature"]) { wallData.temperature_ = @@ -932,6 +952,11 @@ namespace YAML sierra::nalu::SpecDissRate>(); inflowData.sdrSpec_ = true; } + if (node["total_dissipation_rate"]) { + inflowData.tdr_ = + node["total_dissipation_rate"].as(); + inflowData.tdrSpec_ = true; + } if (node["mixture_fraction"]) { inflowData.mixFrac_ = node["mixture_fraction"].as< @@ -1008,6 +1033,11 @@ namespace YAML sierra::nalu::SpecDissRate>(); openData.sdrSpec_ = true; } + if (node["total_dissipation_rate"]) { + openData.tdr_ = + node["total_dissipation_rate"].as(); + openData.tdrSpec_ = true; + } if (node["pressure"]) { openData.p_ = node["pressure"].as(); diff --git a/src/ShearStressTransportEquationSystem.C b/src/ShearStressTransportEquationSystem.C index fca8c79b4d..62d9402217 100644 --- a/src/ShearStressTransportEquationSystem.C +++ b/src/ShearStressTransportEquationSystem.C @@ -67,9 +67,9 @@ namespace nalu{ ShearStressTransportEquationSystem::ShearStressTransportEquationSystem( EquationSystems& eqSystems) : EquationSystem(eqSystems, "ShearStressTransportWrap"), - tkeEqSys_(NULL), - sdrEqSys_(NULL), - gammaEqSys_(NULL), + tkeEqSys_(std::make_unique(eqSystems)), + sdrEqSys_(std::make_unique(eqSystems)), + gammaEqSys_(realm_.solutionOptions_->gammaEqActive_? std::make_unique(eqSystems) : nullptr ), tke_(NULL), sdr_(NULL), gamma_(NULL), @@ -82,11 +82,6 @@ ShearStressTransportEquationSystem::ShearStressTransportEquationSystem( { // push back EQ to manager realm_.push_equation_to_systems(this); - - // create momentum and pressure - tkeEqSys_= new TurbKineticEnergyEquationSystem(eqSystems); - sdrEqSys_ = new SpecificDissipationRateEquationSystem(eqSystems); - if (realm_.solutionOptions_->gammaEqActive_) gammaEqSys_ = new GammaEquationSystem(eqSystems); } //-------------------------------------------------------------------------- @@ -409,6 +404,9 @@ ShearStressTransportEquationSystem::update_and_clip() const double tkeMinVal = tkeMinValue_; const double sdrMinVal = sdrMinValue_; + tkeNp1.sync_to_device(); + sdrNp1.sync_to_device(); + nalu_ngp::run_entity_algorithm( "SST::update_and_clip", ngpMesh, stk::topology::NODE_RANK, sel, KOKKOS_LAMBDA(const MeshIndex& mi) { @@ -444,7 +442,9 @@ ShearStressTransportEquationSystem::update_and_clip_gamma() stk::mesh::selectField(*turbViscosity); const double gammaMinVal = gammaMinValue_; const double gammaMaxVal = gammaMaxValue_; - + + gammaNp1.sync_to_device(); + nalu_ngp::run_entity_algorithm( "SST_GammaEQActive::update_and_clip", ngpMesh, stk::topology::NODE_RANK, sel, KOKKOS_LAMBDA(const MeshIndex& mi) { @@ -576,6 +576,8 @@ ShearStressTransportEquationSystem::compute_f_one_blending() (meta.locally_owned_part() | meta.globally_shared_part()) & (stk::mesh::selectField(*fOneBlending_)); + fOneBlend.sync_to_device(); + nalu_ngp::run_entity_algorithm( "SST::compute_fone_blending", ngpMesh, stk::topology::NODE_RANK, sel, KOKKOS_LAMBDA(const MeshIndex& mi) { diff --git a/src/SolutionOptions.C b/src/SolutionOptions.C index 26eae9b86c..b9988e599d 100644 --- a/src/SolutionOptions.C +++ b/src/SolutionOptions.C @@ -613,6 +613,13 @@ SolutionOptions::initialize_turbulence_constants() turbModelConstantMap_[TM_ams_peclet_offset] = 0.6; turbModelConstantMap_[TM_ams_peclet_slope] = 12.0; turbModelConstantMap_[TM_ams_peclet_scale] = 100.0; + turbModelConstantMap_[TM_fMuExp] = -0.0115; + turbModelConstantMap_[TM_utau] = 1.0; + turbModelConstantMap_[TM_cEpsOne] = 1.35; + turbModelConstantMap_[TM_cEpsTwo] = 1.80; + turbModelConstantMap_[TM_fOne] = 1.0; + turbModelConstantMap_[TM_sigmaK] = 1.0; + turbModelConstantMap_[TM_sigmaEps] = 1.3; turbModelConstantMap_[TM_tkeAmb] = 0.0; turbModelConstantMap_[TM_sdrAmb] = 0.0; turbModelConstantMap_[TM_avgTimeCoeff] = 1.0; diff --git a/src/SpecificDissipationRateEquationSystem.C b/src/SpecificDissipationRateEquationSystem.C index 8afebc5b9b..d5a8929420 100644 --- a/src/SpecificDissipationRateEquationSystem.C +++ b/src/SpecificDissipationRateEquationSystem.C @@ -52,6 +52,7 @@ #include #include #include +#include // ngp #include "ngp_utils/NgpFieldBLAS.h" @@ -59,6 +60,7 @@ #include "ngp_algorithms/NodalGradElemAlg.h" #include "ngp_algorithms/NodalGradBndryElemAlg.h" #include "ngp_algorithms/EffSSTDiffFluxCoeffAlg.h" +#include "ngp_algorithms/EffDiffFluxCoeffAlg.h" #include "ngp_algorithms/SDRWallFuncAlg.h" #include "ngp_algorithms/SDRLowReWallAlg.h" #include "ngp_algorithms/SDRWallFuncAlgDriver.h" @@ -258,7 +260,6 @@ SpecificDissipationRateEquationSystem::register_interior_algorithm( nodeAlg.add_kernel(realm_.bulk_data(), sdr_); if (TurbulenceModel::SST == realm_.solutionOptions_->turbulenceModel_) { - NaluEnv::self().naluOutputP0() << "call SDRSSTNodeKernel1: " <(realm_.meta_data()); } else if ( (TurbulenceModel::SST_DES == @@ -267,13 +268,14 @@ SpecificDissipationRateEquationSystem::register_interior_algorithm( realm_.solutionOptions_->turbulenceModel_)) { nodeAlg.add_kernel(realm_.meta_data()); } else if ( - TurbulenceModel::SST_AMS == realm_.solutionOptions_->turbulenceModel_) + TurbulenceModel::SST_AMS == realm_.solutionOptions_->turbulenceModel_){ nodeAlg.add_kernel( realm_.meta_data(), realm_.solutionOptions_->get_coordinates_name()); - else { - nodeAlg.add_kernel(realm_.meta_data()); - NaluEnv::self().naluOutputP0() << "call SDRSSTNodeKernel2: " <turbulenceModel_) { + nodeAlg.add_kernel(realm_.meta_data()); + } else { + throw std::runtime_error("Invalid turbulence model in SDR equation system: " + TurbulenceModelNames[static_cast(realm_.solutionOptions_->turbulenceModel_)]); } }, [&](AssembleNGPNodeSolverAlgorithm& nodeAlg, std::string& srcName) { @@ -292,10 +294,15 @@ SpecificDissipationRateEquationSystem::register_interior_algorithm( // effective diffusive flux coefficient alg for SST if (!effDiffFluxAlg_) { - const double sigmaWOne = realm_.get_turb_model_constant(TM_sigmaWOne); - const double sigmaWTwo = realm_.get_turb_model_constant(TM_sigmaWTwo); - effDiffFluxAlg_.reset(new EffSSTDiffFluxCoeffAlg( - realm_, part, visc_, tvisc_, evisc_, sigmaWOne, sigmaWTwo)); + if (TurbulenceModel::KO == realm_.solutionOptions_->turbulenceModel_) { + effDiffFluxAlg_.reset(new EffDiffFluxCoeffAlg( + realm_, part, visc_, tvisc_, evisc_, 1.0, 2.0, realm_.is_turbulent())); + } else { + const double sigmaWOne = realm_.get_turb_model_constant(TM_sigmaWOne); + const double sigmaWTwo = realm_.get_turb_model_constant(TM_sigmaWTwo); + effDiffFluxAlg_.reset(new EffSSTDiffFluxCoeffAlg( + realm_, part, visc_, tvisc_, evisc_, sigmaWOne, sigmaWTwo)); + } } else { effDiffFluxAlg_->partVec_.push_back(part); } @@ -663,7 +670,6 @@ SpecificDissipationRateEquationSystem::predict_state() (meta.locally_owned_part() | meta.globally_shared_part() | meta.aura_part()) & stk::mesh::selectField(*sdr_); nalu_ngp::field_copy(ngpMesh, sel, sdrNp1, sdrN); - sdrNp1.modify_on_device(); } } // namespace nalu diff --git a/src/TotalDissipationRateEquationSystem.C b/src/TotalDissipationRateEquationSystem.C new file mode 100644 index 0000000000..f8f6d91160 --- /dev/null +++ b/src/TotalDissipationRateEquationSystem.C @@ -0,0 +1,614 @@ +// 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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// template for supp algs +#include +#include +#include + +// consolidated +#include + +// edge kernels +#include +#include + +// node kernels +#include +#include +#include +#include + +// ngp +#include "ngp_utils/NgpFieldBLAS.h" +#include "ngp_algorithms/NodalGradEdgeAlg.h" +#include "ngp_algorithms/NodalGradElemAlg.h" +#include "ngp_algorithms/NodalGradBndryElemAlg.h" +#include "ngp_algorithms/EffDiffFluxCoeffAlg.h" +#include "utils/StkHelpers.h" + +#include + +// stk_util +#include + +// stk_mesh/base/fem +#include +#include +#include +#include +#include +#include +#include + +// stk_io +#include + +// stk_topo +#include + +// stk_util +#include + +namespace sierra { +namespace nalu { + +//========================================================================== +// Class Definition +//========================================================================== +// TotalDissipationRateEquationSystem - manages tdr pde system +//========================================================================== +//-------------------------------------------------------------------------- +//-------- constructor ----------------------------------------------------- +//-------------------------------------------------------------------------- +TotalDissipationRateEquationSystem::TotalDissipationRateEquationSystem( + EquationSystems& eqSystems) + : EquationSystem(eqSystems, "TotDissRateEQS", "total_dissipation_rate"), + managePNG_(realm_.get_consistent_mass_matrix_png("total_dissipation_rate")), + tdr_(NULL), + dedx_(NULL), + eTmp_(NULL), + visc_(NULL), + tvisc_(NULL), + evisc_(NULL), + tdrWallBc_(NULL), + assembledWallTdr_(NULL), + assembledWallArea_(NULL), + nodalGradAlgDriver_(realm_, "dedx") +{ + dofName_ = "total_dissipation_rate"; + + // extract solver name and solver object + std::string solverName = + realm_.equationSystems_.get_solver_block_name("total_dissipation_rate"); + LinearSolver* solver = realm_.root()->linearSolvers_->create_solver( + solverName, realm_.name(), EQ_TOT_DISS_RATE); + linsys_ = LinearSystem::create(realm_, 1, this, solver); + + // determine nodal gradient form + set_nodal_gradient("total_dissipation_rate"); + NaluEnv::self().naluOutputP0() + << "Edge projected nodal gradient for total_dissipation_rate: " + << edgeNodalGradient_ << std::endl; + + // push back EQ to manager + realm_.push_equation_to_systems(this); + + // create projected nodal gradient equation system + if (managePNG_) + throw std::runtime_error( + "TotalDissipationRateEquationSystem::Error managePNG is not complete"); +} + +//-------------------------------------------------------------------------- +//-------- destructor ------------------------------------------------------ +//-------------------------------------------------------------------------- +TotalDissipationRateEquationSystem::~TotalDissipationRateEquationSystem() = + default; + +//-------------------------------------------------------------------------- +//-------- register_nodal_fields ------------------------------------------- +//-------------------------------------------------------------------------- +void +TotalDissipationRateEquationSystem::register_nodal_fields(stk::mesh::Part* part) +{ + + stk::mesh::MetaData& meta_data = realm_.meta_data(); + + const int nDim = meta_data.spatial_dimension(); + const int numStates = realm_.number_of_states(); + + // register dof; set it as a restart variable + tdr_ = &(meta_data.declare_field( + stk::topology::NODE_RANK, "total_dissipation_rate", numStates)); + stk::mesh::put_field_on_mesh(*tdr_, *part, nullptr); + realm_.augment_restart_variable_list("total_dissipation_rate"); + + dedx_ = &( + meta_data.declare_field(stk::topology::NODE_RANK, "dedx")); + stk::mesh::put_field_on_mesh(*dedx_, *part, nDim, nullptr); + + // delta solution for linear solver; share delta since this is a split system + eTmp_ = &( + meta_data.declare_field(stk::topology::NODE_RANK, "eTmp")); + stk::mesh::put_field_on_mesh(*eTmp_, *part, nullptr); + + visc_ = &(meta_data.declare_field( + stk::topology::NODE_RANK, "viscosity")); + stk::mesh::put_field_on_mesh(*visc_, *part, nullptr); + + tvisc_ = &(meta_data.declare_field( + stk::topology::NODE_RANK, "turbulent_viscosity")); + stk::mesh::put_field_on_mesh(*tvisc_, *part, nullptr); + + evisc_ = &(meta_data.declare_field( + stk::topology::NODE_RANK, "effective_viscosity_tdr")); + stk::mesh::put_field_on_mesh(*evisc_, *part, nullptr); + + // make sure all states are properly populated (restart can handle this) + if ( + numStates > 2 && + (!realm_.restarted_simulation() || realm_.support_inconsistent_restart())) { + ScalarFieldType& tdrN = tdr_->field_of_state(stk::mesh::StateN); + ScalarFieldType& tdrNp1 = tdr_->field_of_state(stk::mesh::StateNP1); + + CopyFieldAlgorithm* theCopyAlg = new CopyFieldAlgorithm( + realm_, part, &tdrNp1, &tdrN, 0, 1, stk::topology::NODE_RANK); + copyStateAlg_.push_back(theCopyAlg); + } +} + +//-------------------------------------------------------------------------- +//-------- register_interior_algorithm ------------------------------------- +//-------------------------------------------------------------------------- +void +TotalDissipationRateEquationSystem::register_interior_algorithm( + stk::mesh::Part* part) +{ + + // types of algorithms + const AlgorithmType algType = INTERIOR; + + ScalarFieldType& tdrNp1 = tdr_->field_of_state(stk::mesh::StateNP1); + VectorFieldType& dedxNone = dedx_->field_of_state(stk::mesh::StateNone); + + if (edgeNodalGradient_ && realm_.realmUsesEdges_) + nodalGradAlgDriver_.register_edge_algorithm( + algType, part, "tdr_nodal_grad", &tdrNp1, &dedxNone); + else + nodalGradAlgDriver_.register_elem_algorithm( + algType, part, "tdr_nodal_grad", &tdrNp1, &dedxNone, edgeNodalGradient_); + + // solver; interior contribution (advection + diffusion) + if (!realm_.solutionOptions_->useConsolidatedSolverAlg_) { + + std::map::iterator itsi = + solverAlgDriver_->solverAlgMap_.find(algType); + if (itsi == solverAlgDriver_->solverAlgMap_.end()) { + SolverAlgorithm* theAlg = NULL; + if (realm_.realmUsesEdges_) { + const bool useAvgMdot = realm_.solutionOptions_->turbulenceModel_ == TurbulenceModel::SST_AMS; + theAlg = new ScalarEdgeSolverAlg( + realm_, part, this, tdr_, dedx_, evisc_, useAvgMdot); + } else { + throw std::runtime_error( + "TDREQS: Attempt to use non-NGP element solver algorithm"); + } + solverAlgDriver_->solverAlgMap_[algType] = theAlg; + + // look for fully integrated source terms + std::map>::iterator isrc = + realm_.solutionOptions_->elemSrcTermsMap_.find( + "total_dissipation_rate"); + if (isrc != realm_.solutionOptions_->elemSrcTermsMap_.end()) { + throw std::runtime_error( + "TotalDissipationElemSrcTerms::Error can not use element source " + "terms for an edge-based scheme"); + } + } else { + itsi->second->partVec_.push_back(part); + } + + // Check if the user has requested CMM or LMM algorithms; if so, do not + // include Nodal Mass algorithms + std::vector checkAlgNames = { + "total_dissipation_rate_time_derivative", + "lumped_total_dissipation_rate_time_derivative"}; + bool elementMassAlg = supp_alg_is_requested(checkAlgNames); + auto& solverAlgMap = solverAlgDriver_->solverAlgMap_; + process_ngp_node_kernels( + solverAlgMap, realm_, part, this, + [&](AssembleNGPNodeSolverAlgorithm& nodeAlg) { + if (!elementMassAlg) + nodeAlg.add_kernel(realm_.bulk_data(), tdr_); + + if (TurbulenceModel::KE == realm_.solutionOptions_->turbulenceModel_) { + nodeAlg.add_kernel(realm_.meta_data()); + } else { + nodeAlg.add_kernel(realm_.meta_data()); + } + }, + [&](AssembleNGPNodeSolverAlgorithm& nodeAlg, std::string& srcName) { + if (srcName == "gcl") { + nodeAlg.add_kernel(realm_.bulk_data(), tdr_); + NaluEnv::self().naluOutputP0() << " - " << srcName << std::endl; + } else + throw std::runtime_error("TDREqSys: Invalid source term: " + srcName); + }); + } else { + throw std::runtime_error("TDREQS: Element terms not supported"); + } + + // effective diffusive flux coefficient alg for SST + if (!effDiffFluxAlg_) { + const double sigmaEps = realm_.get_turb_model_constant(TM_sigmaEps); + effDiffFluxAlg_.reset(new EffDiffFluxCoeffAlg( + realm_, part, visc_, tvisc_, evisc_, 1.0, sigmaEps, + realm_.is_turbulent())); + } else { + effDiffFluxAlg_->partVec_.push_back(part); + } +} + +//-------------------------------------------------------------------------- +//-------- register_inflow_bc ---------------------------------------------- +//-------------------------------------------------------------------------- +void +TotalDissipationRateEquationSystem::register_inflow_bc( + stk::mesh::Part* part, + const stk::topology& /*theTopo*/, + const InflowBoundaryConditionData& inflowBCData) +{ + + // algorithm type + const AlgorithmType algType = INFLOW; + + ScalarFieldType& tdrNp1 = tdr_->field_of_state(stk::mesh::StateNP1); + VectorFieldType& dedxNone = dedx_->field_of_state(stk::mesh::StateNone); + + stk::mesh::MetaData& meta_data = realm_.meta_data(); + + // register boundary data; tdr_bc + ScalarFieldType* theBcField = &(meta_data.declare_field( + stk::topology::NODE_RANK, "tdr_bc")); + stk::mesh::put_field_on_mesh(*theBcField, *part, nullptr); + + // extract the value for user specified tke and save off the AuxFunction + InflowUserData userData = inflowBCData.userData_; + TotDissRate tdr = userData.tdr_; + std::vector userSpec(1); + userSpec[0] = tdr.totDissRate_; + + // new it + ConstantAuxFunction* theAuxFunc = new ConstantAuxFunction(0, 1, userSpec); + + // 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 tdr_bc to total_dissipation_rate np1... + CopyFieldAlgorithm* theCopyAlg = new CopyFieldAlgorithm( + realm_, part, theBcField, &tdrNp1, 0, 1, stk::topology::NODE_RANK); + bcDataMapAlg_.push_back(theCopyAlg); + + // non-solver; dedx; allow for element-based shifted + nodalGradAlgDriver_.register_face_algorithm( + algType, part, "tdr_nodal_grad", &tdrNp1, &dedxNone, edgeNodalGradient_); + + // Dirichlet bc + std::map::iterator itd = + solverAlgDriver_->solverDirichAlgMap_.find(algType); + if (itd == solverAlgDriver_->solverDirichAlgMap_.end()) { + DirichletBC* theAlg = + new DirichletBC(realm_, this, part, &tdrNp1, theBcField, 0, 1); + solverAlgDriver_->solverDirichAlgMap_[algType] = theAlg; + } else { + itd->second->partVec_.push_back(part); + } +} + +//-------------------------------------------------------------------------- +//-------- register_open_bc ------------------------------------------------ +//-------------------------------------------------------------------------- +void +TotalDissipationRateEquationSystem::register_open_bc( + stk::mesh::Part* part, + const stk::topology& partTopo, + const OpenBoundaryConditionData& openBCData) +{ + + // algorithm type + const AlgorithmType algType = OPEN; + + ScalarFieldType& tdrNp1 = tdr_->field_of_state(stk::mesh::StateNP1); + VectorFieldType& dedxNone = dedx_->field_of_state(stk::mesh::StateNone); + + stk::mesh::MetaData& meta_data = realm_.meta_data(); + + // register boundary data; tdr_bc + ScalarFieldType* theBcField = &(meta_data.declare_field( + stk::topology::NODE_RANK, "open_tdr_bc")); + stk::mesh::put_field_on_mesh(*theBcField, *part, nullptr); + + // extract the value for user specified tke and save off the AuxFunction + OpenUserData userData = openBCData.userData_; + TotDissRate tdr = userData.tdr_; + std::vector userSpec(1); + userSpec[0] = tdr.totDissRate_; + + // new it + ConstantAuxFunction* theAuxFunc = new ConstantAuxFunction(0, 1, userSpec); + + // bc data alg + AuxFunctionAlgorithm* auxAlg = new AuxFunctionAlgorithm( + realm_, part, theBcField, theAuxFunc, stk::topology::NODE_RANK); + bcDataAlg_.push_back(auxAlg); + + // non-solver; dedx; allow for element-based shifted + nodalGradAlgDriver_.register_face_algorithm( + algType, part, "tdr_nodal_grad", &tdrNp1, &dedxNone, edgeNodalGradient_); + + if (realm_.realmUsesEdges_) { + auto& solverAlgMap = solverAlgDriver_->solverAlgorithmMap_; + AssembleElemSolverAlgorithm* elemSolverAlg = nullptr; + bool solverAlgWasBuilt = false; + + std::tie(elemSolverAlg, solverAlgWasBuilt) = + build_or_add_part_to_face_bc_solver_alg( + *this, *part, solverAlgMap, "open"); + + auto& dataPreReqs = elemSolverAlg->dataNeededByKernels_; + auto& activeKernels = elemSolverAlg->activeKernels_; + + build_face_topo_kernel_automatic( + partTopo, *this, activeKernels, "tdr_open", realm_.meta_data(), + *realm_.solutionOptions_, tdr_, theBcField, dataPreReqs); + } else { + throw std::runtime_error( + "TDREQS: Attempt to use non-NGP element open algorithm"); + } +} + +//-------------------------------------------------------------------------- +//-------- register_wall_bc ------------------------------------------------ +//-------------------------------------------------------------------------- +void +TotalDissipationRateEquationSystem::register_wall_bc( + stk::mesh::Part* part, + const stk::topology& /*theTopo*/, + const WallBoundaryConditionData& /*wallBCData*/) +{ + + // algorithm type + const AlgorithmType algType = WALL; + + // np1 + ScalarFieldType& tdrNp1 = tdr_->field_of_state(stk::mesh::StateNP1); + VectorFieldType& dedxNone = dedx_->field_of_state(stk::mesh::StateNone); + + stk::mesh::MetaData& meta_data = realm_.meta_data(); + + // register boundary data; tdr_bc + tdrWallBc_ = &(meta_data.declare_field( + stk::topology::NODE_RANK, "tdr_bc")); + stk::mesh::put_field_on_mesh(*tdrWallBc_, *part, nullptr); + + // need to register the assembles wall value for tdr; can not share with + // tdr_bc + assembledWallTdr_ = &(meta_data.declare_field( + stk::topology::NODE_RANK, "wall_model_tdr_bc")); + stk::mesh::put_field_on_mesh(*assembledWallTdr_, *part, nullptr); + + assembledWallArea_ = &(meta_data.declare_field( + stk::topology::NODE_RANK, "assembled_wall_area_tdr")); + stk::mesh::put_field_on_mesh(*assembledWallArea_, *part, nullptr); + + // Dirichlet bc + std::map::iterator itd = + solverAlgDriver_->solverDirichAlgMap_.find(algType); + if (itd == solverAlgDriver_->solverDirichAlgMap_.end()) { + DirichletBC* theAlg = + new DirichletBC(realm_, this, part, &tdrNp1, tdrWallBc_, 0, 1); + solverAlgDriver_->solverDirichAlgMap_[algType] = theAlg; + } else { + itd->second->partVec_.push_back(part); + } + + // non-solver; dedx; allow for element-based shifted + nodalGradAlgDriver_.register_face_algorithm( + algType, part, "tdr_nodal_grad", &tdrNp1, &dedxNone, edgeNodalGradient_); +} + +//-------------------------------------------------------------------------- +//-------- register_symmetry_bc -------------------------------------------- +//-------------------------------------------------------------------------- +void +TotalDissipationRateEquationSystem::register_symmetry_bc( + stk::mesh::Part* part, + const stk::topology& /*theTopo*/, + const SymmetryBoundaryConditionData& /* symmetryBCData */) +{ + + // algorithm type + const AlgorithmType algType = SYMMETRY; + + // np1 + ScalarFieldType& tdrNp1 = tdr_->field_of_state(stk::mesh::StateNP1); + VectorFieldType& dedxNone = dedx_->field_of_state(stk::mesh::StateNone); + + // non-solver; dedx; allow for element-based shifted + nodalGradAlgDriver_.register_face_algorithm( + algType, part, "tdr_nodal_grad", &tdrNp1, &dedxNone, edgeNodalGradient_); +} + +//-------------------------------------------------------------------------- +//-------- register_non_conformal_bc --------------------------------------- +//-------------------------------------------------------------------------- +void +TotalDissipationRateEquationSystem::register_non_conformal_bc( + stk::mesh::Part* part, const stk::topology& /*theTopo*/) +{ + + const AlgorithmType algType = NON_CONFORMAL; + + // np1 + ScalarFieldType& tdrNp1 = tdr_->field_of_state(stk::mesh::StateNP1); + VectorFieldType& dedxNone = dedx_->field_of_state(stk::mesh::StateNone); + + // non-solver; contribution to dedx; DG algorithm decides on locations for + // integration points + if (edgeNodalGradient_) { + nodalGradAlgDriver_.register_face_algorithm( + algType, part, "tdr_nodal_grad", &tdrNp1, &dedxNone, edgeNodalGradient_); + } else { + // proceed with DG + nodalGradAlgDriver_ + .register_legacy_algorithm( + algType, part, "tdr_nodal_grad", &tdrNp1, &dedxNone); + } + + // solver; lhs; same for edge and element-based scheme + std::map::iterator itsi = + solverAlgDriver_->solverAlgMap_.find(algType); + if (itsi == solverAlgDriver_->solverAlgMap_.end()) { + AssembleScalarNonConformalSolverAlgorithm* theAlg = + new AssembleScalarNonConformalSolverAlgorithm( + realm_, part, this, tdr_, evisc_); + solverAlgDriver_->solverAlgMap_[algType] = theAlg; + } else { + itsi->second->partVec_.push_back(part); + } +} + +//-------------------------------------------------------------------------- +//-------- register_overset_bc --------------------------------------------- +//-------------------------------------------------------------------------- +void +TotalDissipationRateEquationSystem::register_overset_bc() +{ + create_constraint_algorithm(tdr_); + + equationSystems_.register_overset_field_update(tdr_, 1, 1); +} + +//-------------------------------------------------------------------------- +//-------- initialize ------------------------------------------------------ +//-------------------------------------------------------------------------- +void +TotalDissipationRateEquationSystem::initialize() +{ + solverAlgDriver_->initialize_connectivity(); + linsys_->finalizeLinearSystem(); +} + +//-------------------------------------------------------------------------- +//-------- reinitialize_linear_system -------------------------------------- +//-------------------------------------------------------------------------- +void +TotalDissipationRateEquationSystem::reinitialize_linear_system() +{ + // If this is decoupled overset simulation and the user has requested that the + // linear system be reused, then do nothing + if (decoupledOverset_ && linsys_->config().reuseLinSysIfPossible()) + return; + + // delete linsys + delete linsys_; + + // create new solver + std::string solverName = + realm_.equationSystems_.get_solver_block_name("total_dissipation_rate"); + LinearSolver* solver = realm_.root()->linearSolvers_->reinitialize_solver( + solverName, realm_.name(), EQ_TOT_DISS_RATE); + linsys_ = LinearSystem::create(realm_, 1, this, solver); + + // initialize + solverAlgDriver_->initialize_connectivity(); + linsys_->finalizeLinearSystem(); +} + +//-------------------------------------------------------------------------- +//-------- assemble_nodal_gradient() --------------------------------------- +//-------------------------------------------------------------------------- +void +TotalDissipationRateEquationSystem::assemble_nodal_gradient() +{ + const double timeA = -NaluEnv::self().nalu_time(); + nodalGradAlgDriver_.execute(); + timerMisc_ += (NaluEnv::self().nalu_time() + timeA); +} + +//-------------------------------------------------------------------------- +//-------- compute_effective_flux_coeff() ---------------------------------- +//-------------------------------------------------------------------------- +void +TotalDissipationRateEquationSystem::compute_effective_diff_flux_coeff() +{ + const double timeA = -NaluEnv::self().nalu_time(); + effDiffFluxAlg_->execute(); + timerMisc_ += (NaluEnv::self().nalu_time() + timeA); +} + +//-------------------------------------------------------------------------- +//-------- predict_state() ------------------------------------------------- +//-------------------------------------------------------------------------- +void +TotalDissipationRateEquationSystem::predict_state() +{ + const auto& ngpMesh = realm_.ngp_mesh(); + const auto& fieldMgr = realm_.ngp_field_manager(); + const auto& tdrN = fieldMgr.get_field( + tdr_->field_of_state(stk::mesh::StateN).mesh_meta_data_ordinal()); + auto& tdrNp1 = fieldMgr.get_field( + tdr_->field_of_state(stk::mesh::StateNP1).mesh_meta_data_ordinal()); + + const auto& meta = realm_.meta_data(); + const stk::mesh::Selector sel = + (meta.locally_owned_part() | meta.globally_shared_part() | + meta.aura_part()) & + stk::mesh::selectField(*tdr_); + nalu_ngp::field_copy(ngpMesh, sel, tdrNp1, tdrN); +} + +} // namespace nalu +} // namespace sierra diff --git a/src/TurbKineticEnergyEquationSystem.C b/src/TurbKineticEnergyEquationSystem.C index e72a780239..74036e1476 100644 --- a/src/TurbKineticEnergyEquationSystem.C +++ b/src/TurbKineticEnergyEquationSystem.C @@ -59,7 +59,9 @@ #include #include #include +#include #include +#include // ngp #include @@ -146,7 +148,7 @@ TurbKineticEnergyEquationSystem::TurbKineticEnergyEquationSystem( if (!check_for_valid_turblence_model(turbulenceModel_)) { throw std::runtime_error( "User has requested TurbKinEnergyEqs, however, turbulence model is not " - "KSGS, SST, SST_DES, SST_IDDES, or SST_AMS"); + "KSGS, SST, SST_DES, SST_IDDES, KE, SST_AMS, or KO"); } // create projected nodal gradient equation system @@ -165,6 +167,8 @@ TurbKineticEnergyEquationSystem::check_for_valid_turblence_model( case TurbulenceModel::SST_DES: case TurbulenceModel::SST_AMS: case TurbulenceModel::SST_IDDES: + case TurbulenceModel::KE: + case TurbulenceModel::KO: return true; default: return false; @@ -300,6 +304,11 @@ TurbKineticEnergyEquationSystem::register_interior_algorithm( break; case TurbulenceModel::SST_IDDES: nodeAlg.add_kernel(realm_.meta_data()); + case TurbulenceModel::KE: + nodeAlg.add_kernel(realm_.meta_data()); + break; + case TurbulenceModel::KO: + nodeAlg.add_kernel(realm_.meta_data()); break; default: std::runtime_error("TKEEqSys: Invalid turbulence model"); @@ -347,6 +356,19 @@ TurbKineticEnergyEquationSystem::register_interior_algorithm( realm_, part, visc_, tvisc_, evisc_, sigmaKOne, sigmaKTwo)); break; } + case TurbulenceModel::KO: { + effDiffFluxCoeffAlg_.reset(new EffDiffFluxCoeffAlg( + realm_, part, visc_, tvisc_, evisc_, 1.0, 2.0, realm_.is_turbulent())); + break; + } + case TurbulenceModel::KE: { + const double sigmaK = realm_.get_turb_model_constant(TM_sigmaK); + effDiffFluxCoeffAlg_.reset(new EffDiffFluxCoeffAlg( + realm_, part, visc_, tvisc_, evisc_, 1.0, sigmaK, + realm_.is_turbulent())); + break; + } + default: throw std::runtime_error("Unsupported turbulence model in TurbKe"); } @@ -774,6 +796,8 @@ TurbKineticEnergyEquationSystem::initial_work() (meta.locally_owned_part() | meta.globally_shared_part()) & stk::mesh::selectField(*tke_); + ngpTke.sync_to_device(); + nalu_ngp::run_entity_algorithm( "clip_tke", ngpMesh, stk::topology::NODE_RANK, sel, @@ -802,6 +826,8 @@ TurbKineticEnergyEquationSystem::post_external_data_transfer_work() (meta.locally_owned_part() | meta.globally_shared_part()) & stk::mesh::selectField(*tke_); + ngpTke.sync_to_device(); + nalu_ngp::run_entity_algorithm( "clip_tke", ngpMesh, stk::topology::NODE_RANK, sel, @@ -875,6 +901,8 @@ TurbKineticEnergyEquationSystem::update_and_clip() auto ngpTke = fieldMgr.get_field( tke_->mesh_meta_data_ordinal()); + ngpTke.sync_to_device(); + nalu_ngp::run_entity_par_reduce( "tke_update_and_clip", ngpMesh, stk::topology::NODE_RANK, sel, diff --git a/src/WilcoxKOmegaEquationSystem.C b/src/WilcoxKOmegaEquationSystem.C new file mode 100644 index 0000000000..92c52a3915 --- /dev/null +++ b/src/WilcoxKOmegaEquationSystem.C @@ -0,0 +1,394 @@ +// 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 +#include +#include +#include +#include +#include +#include +#include + +// ngp +#include "FieldTypeDef.h" +#include "ngp_algorithms/GeometryAlgDriver.h" +#include "ngp_algorithms/WallFuncGeometryAlg.h" +#include "ngp_utils/NgpLoopUtils.h" +#include "ngp_utils/NgpFieldUtils.h" +#include "ngp_utils/NgpFieldBLAS.h" + +// stk_util +#include +#include "utils/StkHelpers.h" + +// stk_mesh/base/fem +#include +#include +#include + +#include + +// stk_io +#include + +// basic c++ +#include +#include +#include + +// ngp +#include "ngp_utils/NgpFieldBLAS.h" + +namespace sierra { +namespace nalu { + +//========================================================================== +// Class Definition +//========================================================================== +// WilcoxKOmegaEquationSystem - manage KOmega +//========================================================================== +//-------------------------------------------------------------------------- +//-------- constructor ----------------------------------------------------- +//-------------------------------------------------------------------------- +WilcoxKOmegaEquationSystem::WilcoxKOmegaEquationSystem( + EquationSystems& eqSystems) + : EquationSystem(eqSystems, "WilcoxKOmegaWrap"), + tkeEqSys_(std::make_unique(eqSystems)), + sdrEqSys_(std::make_unique(eqSystems)), + tke_(NULL), + sdr_(NULL), + minDistanceToWall_(NULL), + isInit_(true), + resetAMSAverages_(realm_.solutionOptions_->resetAMSAverages_) +{ + // push back EQ to manager + realm_.push_equation_to_systems(this); +} + +//-------------------------------------------------------------------------- +//-------- destructor ------------------------------------------------------ +//-------------------------------------------------------------------------- +WilcoxKOmegaEquationSystem::~WilcoxKOmegaEquationSystem() {} + +void +WilcoxKOmegaEquationSystem::load(const YAML::Node& node) +{ + EquationSystem::load(node); + + if (realm_.query_for_overset()) { + tkeEqSys_->decoupledOverset_ = decoupledOverset_; + tkeEqSys_->numOversetIters_ = numOversetIters_; + sdrEqSys_->decoupledOverset_ = decoupledOverset_; + sdrEqSys_->numOversetIters_ = numOversetIters_; + } +} + +//-------------------------------------------------------------------------- +//-------- initialize ------------------------------------------------------ +//-------------------------------------------------------------------------- +void +WilcoxKOmegaEquationSystem::initialize() +{ + // let equation systems that are owned some information + tkeEqSys_->convergenceTolerance_ = convergenceTolerance_; + sdrEqSys_->convergenceTolerance_ = convergenceTolerance_; +} + +//-------------------------------------------------------------------------- +//-------- register_nodal_fields ------------------------------------------- +//-------------------------------------------------------------------------- +void +WilcoxKOmegaEquationSystem::register_nodal_fields(stk::mesh::Part* part) +{ + + stk::mesh::MetaData& meta_data = realm_.meta_data(); + const int numStates = realm_.number_of_states(); + + // re-register tke and sdr for convenience + tke_ = &(meta_data.declare_field( + stk::topology::NODE_RANK, "turbulent_ke", numStates)); + stk::mesh::put_field_on_mesh(*tke_, *part, nullptr); + sdr_ = &(meta_data.declare_field( + stk::topology::NODE_RANK, "specific_dissipation_rate", numStates)); + stk::mesh::put_field_on_mesh(*sdr_, *part, nullptr); + + // SST parameters that everyone needs + minDistanceToWall_ = &(meta_data.declare_field( + stk::topology::NODE_RANK, "minimum_distance_to_wall")); + stk::mesh::put_field_on_mesh(*minDistanceToWall_, *part, nullptr); + + // add to restart field + realm_.augment_restart_variable_list("minimum_distance_to_wall"); +} + +//-------------------------------------------------------------------------- +//-------- register_interior_algorithm ------------------------------------- +//-------------------------------------------------------------------------- +void +WilcoxKOmegaEquationSystem::register_interior_algorithm( + stk::mesh::Part* /*part*/) +{ + // nothing to do here... +} + +//-------------------------------------------------------------------------- +//-------- register_wall_bc ------------------------------------------------ +//-------------------------------------------------------------------------- +void +WilcoxKOmegaEquationSystem::register_wall_bc( + stk::mesh::Part* part, + const stk::topology& partTopo, + const WallBoundaryConditionData& wallBCData) +{ + // determine if using RANS for ABL + WallUserData userData = wallBCData.userData_; + bool RANSAblBcApproach = userData.RANSAblBcApproach_; + + // push mesh part + wallBcPart_.push_back(part); + + auto& meta = realm_.meta_data(); + auto& assembledWallArea = meta.declare_field( + stk::topology::NODE_RANK, "assembled_wall_area_wf"); + stk::mesh::put_field_on_mesh(assembledWallArea, *part, nullptr); + auto& assembledWallNormDist = meta.declare_field( + stk::topology::NODE_RANK, "assembled_wall_normal_distance"); + stk::mesh::put_field_on_mesh(assembledWallNormDist, *part, nullptr); + auto& wallNormDistBip = meta.declare_field( + meta.side_rank(), "wall_normal_distance_bip"); + auto* meFC = MasterElementRepo::get_surface_master_element(partTopo); + const int numScsBip = meFC->num_integration_points(); + stk::mesh::put_field_on_mesh(wallNormDistBip, *part, numScsBip, nullptr); +} + +//-------------------------------------------------------------------------- +//-------- solve_and_update ------------------------------------------------ +//-------------------------------------------------------------------------- +void +WilcoxKOmegaEquationSystem::solve_and_update() +{ + // wrap timing + // SST_FIXME: deal with timers; all on misc for SSTEqs double timeA, timeB; + if (isInit_) { + // compute projected nodal gradients + tkeEqSys_->compute_projected_nodal_gradient(); + sdrEqSys_->assemble_nodal_gradient(); + clip_min_distance_to_wall(); + + isInit_ = false; + } else if (realm_.has_mesh_motion()) { + if (realm_.currentNonlinearIteration_ == 1) + clip_min_distance_to_wall(); + } + + // SST effective viscosity for k and omega + tkeEqSys_->compute_effective_diff_flux_coeff(); + sdrEqSys_->compute_effective_diff_flux_coeff(); + + tkeEqSys_->compute_wall_model_parameters(); + sdrEqSys_->compute_wall_model_parameters(); + + // start the iteration loop + for (int k = 0; k < maxIterations_; ++k) { + + NaluEnv::self().naluOutputP0() + << " " << k + 1 << "/" << maxIterations_ << std::setw(15) << std::right + << name_ << std::endl; + + for (int oi = 0; oi < numOversetIters_; ++oi) { + // tke and sdr assemble, load_complete and solve; Jacobi iteration + tkeEqSys_->assemble_and_solve(tkeEqSys_->kTmp_); + sdrEqSys_->assemble_and_solve(sdrEqSys_->wTmp_); + + update_and_clip(); + + if (decoupledOverset_ && realm_.hasOverset_) { + realm_.overset_field_update(tkeEqSys_->tke_, 1, 1); + realm_.overset_field_update(sdrEqSys_->sdr_, 1, 1); + } + } + // compute projected nodal gradients + tkeEqSys_->compute_projected_nodal_gradient(); + sdrEqSys_->assemble_nodal_gradient(); + } +} + +/** Perform sanity checks on TKE/TDR fields + */ +void +WilcoxKOmegaEquationSystem::initial_work() +{ + const auto& meshInfo = realm_.mesh_info(); + const auto& meta = meshInfo.meta(); + const auto& ngpMesh = meshInfo.ngp_mesh(); + const auto& fieldMgr = meshInfo.ngp_field_manager(); + + auto& tkeNp1 = fieldMgr.get_field(tke_->mesh_meta_data_ordinal()); + auto& sdrNp1 = fieldMgr.get_field(sdr_->mesh_meta_data_ordinal()); + const stk::mesh::Selector sel = + (meta.locally_owned_part() | meta.globally_shared_part()) & + stk::mesh::selectField(*sdr_); + clip_ko(ngpMesh, sel, tkeNp1, sdrNp1); +} + +void +WilcoxKOmegaEquationSystem::post_external_data_transfer_work() +{ + const auto& meshInfo = realm_.mesh_info(); + const auto& meta = meshInfo.meta(); + const auto& ngpMesh = meshInfo.ngp_mesh(); + const auto& fieldMgr = meshInfo.ngp_field_manager(); + + auto& tkeNp1 = fieldMgr.get_field(tke_->mesh_meta_data_ordinal()); + auto& sdrNp1 = fieldMgr.get_field(sdr_->mesh_meta_data_ordinal()); + + const stk::mesh::Selector owned_and_shared = + (meta.locally_owned_part() | meta.globally_shared_part()); + auto interior_sel = owned_and_shared & stk::mesh::selectField(*sdr_); + clip_ko(ngpMesh, interior_sel, tkeNp1, sdrNp1); + + auto sdrBCField = + meta.get_field(stk::topology::NODE_RANK, "sdr_bc"); + auto tkeBCField = + meta.get_field(stk::topology::NODE_RANK, "tke_bc"); + if (sdrBCField != nullptr) { + ThrowRequire(tkeBCField); + auto bc_sel = owned_and_shared & stk::mesh::selectField(*sdrBCField); + auto ngpTkeBC = + fieldMgr.get_field(tkeBCField->mesh_meta_data_ordinal()); + auto ngpTdrBC = + fieldMgr.get_field(sdrBCField->mesh_meta_data_ordinal()); + clip_ko(ngpMesh, bc_sel, ngpTkeBC, ngpTdrBC); + } +} + +/** Update solution but ensure that TKE and TDR are greater than zero + */ +void +WilcoxKOmegaEquationSystem::update_and_clip() +{ + using MeshIndex = nalu_ngp::NGPMeshTraits<>::MeshIndex; + + const auto& meshInfo = realm_.mesh_info(); + const auto& meta = meshInfo.meta(); + const auto& ngpMesh = meshInfo.ngp_mesh(); + const auto& fieldMgr = meshInfo.ngp_field_manager(); + + auto& tkeNp1 = fieldMgr.get_field(tke_->mesh_meta_data_ordinal()); + auto& sdrNp1 = fieldMgr.get_field(sdr_->mesh_meta_data_ordinal()); + const auto& kTmp = + fieldMgr.get_field(tkeEqSys_->kTmp_->mesh_meta_data_ordinal()); + const auto& wTmp = + fieldMgr.get_field(sdrEqSys_->wTmp_->mesh_meta_data_ordinal()); + + auto* turbViscosity = meta.get_field( + stk::topology::NODE_RANK, "turbulent_viscosity"); + + const stk::mesh::Selector sel = + (meta.locally_owned_part() | meta.globally_shared_part()) & + stk::mesh::selectField(*turbViscosity); + + // Bring class variables to local scope for lambda capture + const double tkeMinVal = tkeMinValue_; + const double sdrMinVal = sdrMinValue_; + + tkeNp1.sync_to_device(); + sdrNp1.sync_to_device(); + + nalu_ngp::run_entity_algorithm( + "KE::update_and_clip", ngpMesh, stk::topology::NODE_RANK, sel, + KOKKOS_LAMBDA(const MeshIndex& mi) { + const double tkeNew = tkeNp1.get(mi, 0) + kTmp.get(mi, 0); + const double sdrNew = sdrNp1.get(mi, 0) + wTmp.get(mi, 0); + + tkeNp1.get(mi, 0) = (tkeNew < 0.0) ? tkeMinVal : tkeNew; + sdrNp1.get(mi, 0) = stk::math::max(sdrNew, sdrMinVal); + }); + + tkeNp1.modify_on_device(); + sdrNp1.modify_on_device(); +} + +void +WilcoxKOmegaEquationSystem::clip_ko( + const stk::mesh::NgpMesh& ngpMesh, + const stk::mesh::Selector& sel, + stk::mesh::NgpField& tke, + stk::mesh::NgpField& sdr) +{ + tke.sync_to_device(); + sdr.sync_to_device(); + + // Bring class variables to local scope for lambda capture + const double tkeMinVal = tkeMinValue_; + const double sdrMinVal = sdrMinValue_; + + nalu_ngp::run_entity_algorithm( + "KE::clip", ngpMesh, stk::topology::NODE_RANK, sel, + KOKKOS_LAMBDA(const nalu_ngp::NGPMeshTraits<>::MeshIndex& mi) { + const double tkeNew = tke.get(mi, 0); + const double sdrNew = sdr.get(mi, 0); + + tke.get(mi, 0) = (tkeNew < 0.0) ? tkeMinVal : tkeNew; + sdr.get(mi, 0) = stk::math::max(sdrNew, sdrMinVal); + }); + tke.modify_on_device(); + sdr.modify_on_device(); +} + +//-------------------------------------------------------------------------- +//-------- clip_min_distance_to_wall --------------------------------------- +//-------------------------------------------------------------------------- +void +WilcoxKOmegaEquationSystem::clip_min_distance_to_wall() +{ + using MeshIndex = nalu_ngp::NGPMeshTraits<>::MeshIndex; + const auto& meshInfo = realm_.mesh_info(); + const auto& ngpMesh = meshInfo.ngp_mesh(); + const auto& meta = meshInfo.meta(); + const auto& fieldMgr = meshInfo.ngp_field_manager(); + + auto& ndtw = + fieldMgr.get_field(minDistanceToWall_->mesh_meta_data_ordinal()); + const auto& wallNormDist = + nalu_ngp::get_ngp_field(meshInfo, "assembled_wall_normal_distance"); + + const stk::mesh::Selector sel = + (meta.locally_owned_part() | meta.globally_shared_part()) & + stk::mesh::selectUnion(wallBcPart_); + + ndtw.sync_to_device(); + + nalu_ngp::run_entity_algorithm( + "SST::clip_ndtw", ngpMesh, stk::topology::NODE_RANK, sel, + KOKKOS_LAMBDA(const MeshIndex& mi) { + const double minD = ndtw.get(mi, 0); + + ndtw.get(mi, 0) = stk::math::max(minD, wallNormDist.get(mi, 0)); + }); + ndtw.modify_on_device(); + + stk::mesh::parallel_max(realm_.bulk_data(), {minDistanceToWall_}); + if (realm_.hasPeriodic_) { + realm_.periodic_field_max(minDistanceToWall_, 1); + } +} + +void +WilcoxKOmegaEquationSystem::post_iter_work() +{ + // nothing to do here ... +} + +} // namespace nalu +} // namespace sierra diff --git a/src/ngp_algorithms/CMakeLists.txt b/src/ngp_algorithms/CMakeLists.txt index 3a618b909e..bee20484b8 100644 --- a/src/ngp_algorithms/CMakeLists.txt +++ b/src/ngp_algorithms/CMakeLists.txt @@ -18,6 +18,8 @@ target_sources(nalu PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/MetricTensorElemAlg.C ${CMAKE_CURRENT_SOURCE_DIR}/TurbViscKsgsAlg.C ${CMAKE_CURRENT_SOURCE_DIR}/TurbViscSSTAlg.C + ${CMAKE_CURRENT_SOURCE_DIR}/TurbViscKEAlg.C + ${CMAKE_CURRENT_SOURCE_DIR}/TurbViscKOAlg.C ${CMAKE_CURRENT_SOURCE_DIR}/WallFuncGeometryAlg.C ${CMAKE_CURRENT_SOURCE_DIR}/ABLWallFrictionVelAlg.C ${CMAKE_CURRENT_SOURCE_DIR}/ABLWallFluxesAlg.C diff --git a/src/ngp_algorithms/TurbViscKEAlg.C b/src/ngp_algorithms/TurbViscKEAlg.C new file mode 100644 index 0000000000..69b537bed9 --- /dev/null +++ b/src/ngp_algorithms/TurbViscKEAlg.C @@ -0,0 +1,75 @@ +// 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 "ngp_algorithms/TurbViscKEAlg.h" +#include "ngp_utils/NgpLoopUtils.h" +#include "ngp_utils/NgpTypes.h" +#include "ngp_utils/NgpFieldManager.h" +#include "Realm.h" +#include "utils/StkHelpers.h" +#include "stk_mesh/base/MetaData.hpp" +#include "stk_mesh/base/NgpMesh.hpp" + +namespace sierra { +namespace nalu { + +TurbViscKEAlg::TurbViscKEAlg( + Realm& realm, + stk::mesh::Part* part, + ScalarFieldType* tvisc, + const bool useAverages) + : Algorithm(realm, part), + tviscField_(tvisc), + density_(get_field_ordinal(realm.meta_data(), "density")), + tke_(get_field_ordinal(realm.meta_data(), "turbulent_ke")), + tdr_(get_field_ordinal(realm.meta_data(), "total_dissipation_rate")), + tvisc_(tvisc->mesh_meta_data_ordinal()), + dplus_(get_field_ordinal(realm.meta_data(), "dplus_wall_function")), + cMu_(realm.get_turb_model_constant(TM_cMu)), + fMuExp_(realm.get_turb_model_constant(TM_fMuExp)) +{ +} + +void +TurbViscKEAlg::execute() +{ + using Traits = nalu_ngp::NGPMeshTraits; + + const auto& meta = realm_.meta_data(); + + stk::mesh::Selector sel = + (meta.locally_owned_part() | meta.globally_shared_part()) & + stk::mesh::selectField(*tviscField_); + + const auto& meshInfo = realm_.mesh_info(); + const auto ngpMesh = meshInfo.ngp_mesh(); + const auto& fieldMgr = meshInfo.ngp_field_manager(); + const auto density = fieldMgr.get_field(density_); + const auto tke = fieldMgr.get_field(tke_); + const auto tdr = fieldMgr.get_field(tdr_); + const auto dplus = fieldMgr.get_field(dplus_); + auto tvisc = fieldMgr.get_field(tvisc_); + + const DblType cMu = cMu_; + const DblType fMuExp = fMuExp_; + + nalu_ngp::run_entity_algorithm( + "TurbViscKEAlg", ngpMesh, stk::topology::NODE_RANK, sel, + KOKKOS_LAMBDA(const Traits::MeshIndex& meshIdx) { + const DblType fMu = 1.0 - stk::math::exp(fMuExp * dplus.get(meshIdx, 0)); + + tvisc.get(meshIdx, 0) = cMu * fMu * density.get(meshIdx, 0) * + tke.get(meshIdx, 0) * tke.get(meshIdx, 0) / + stk::math::max(tdr.get(meshIdx, 0), 1.0e-16); + }); + tvisc.modify_on_device(); +} + +} // namespace nalu +} // namespace sierra diff --git a/src/ngp_algorithms/TurbViscKOAlg.C b/src/ngp_algorithms/TurbViscKOAlg.C new file mode 100644 index 0000000000..4b6613be9d --- /dev/null +++ b/src/ngp_algorithms/TurbViscKOAlg.C @@ -0,0 +1,99 @@ +// 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 "ngp_algorithms/TurbViscKOAlg.h" +#include "ngp_utils/NgpLoopUtils.h" +#include "ngp_utils/NgpTypes.h" +#include "ngp_utils/NgpFieldManager.h" +#include "Realm.h" +#include "utils/StkHelpers.h" +#include "stk_mesh/base/MetaData.hpp" +#include "stk_mesh/base/NgpMesh.hpp" + +namespace sierra { +namespace nalu { + +TurbViscKOAlg::TurbViscKOAlg( + Realm& realm, + stk::mesh::Part* part, + ScalarFieldType* tvisc, + const bool useAverages) + : Algorithm(realm, part), + tviscField_(tvisc), + density_(get_field_ordinal(realm.meta_data(), "density")), + viscosity_(get_field_ordinal(realm.meta_data(), "viscosity")), + tke_(get_field_ordinal(realm.meta_data(), "turbulent_ke")), + sdr_(get_field_ordinal(realm.meta_data(), "specific_dissipation_rate")), + minDistance_( + get_field_ordinal(realm.meta_data(), "minimum_distance_to_wall")), + dudx_(get_field_ordinal( + realm.meta_data(), (useAverages) ? "average_dudx" : "dudx")), + tvisc_(tvisc->mesh_meta_data_ordinal()), + aOne_(realm.get_turb_model_constant(TM_aOne)), + betaStar_(realm.get_turb_model_constant(TM_betaStar)) +{ +} + +void +TurbViscKOAlg::execute() +{ + using Traits = nalu_ngp::NGPMeshTraits; + + const auto& meta = realm_.meta_data(); + + stk::mesh::Selector sel = + (meta.locally_owned_part() | meta.globally_shared_part()) & + stk::mesh::selectField(*tviscField_); + + const auto& meshInfo = realm_.mesh_info(); + const auto ngpMesh = meshInfo.ngp_mesh(); + const auto& fieldMgr = meshInfo.ngp_field_manager(); + const auto density = fieldMgr.get_field(density_); + const auto visc = fieldMgr.get_field(viscosity_); + const auto tke = fieldMgr.get_field(tke_); + const auto sdr = fieldMgr.get_field(sdr_); + const auto dudx = fieldMgr.get_field(dudx_); + auto tvisc = fieldMgr.get_field(tvisc_); + + const int nDim = meta.spatial_dimension(); + + nalu_ngp::run_entity_algorithm( + "TurbViscKOAlg", ngpMesh, stk::topology::NODE_RANK, sel, + KOKKOS_LAMBDA(const Traits::MeshIndex& meshIdx) { + DblType sijMag = 0.0; + for (int i = 0; i < nDim; ++i) { + const int offSet = nDim * i; + for (int j = 0; j < nDim; ++j) { + const DblType rateOfStrain = 0.5 * (dudx.get(meshIdx, offSet + j) + + dudx.get(meshIdx, nDim * j + i)); + sijMag += rateOfStrain * rateOfStrain; + } + } + sijMag = stk::math::sqrt(2.0 * sijMag); + + const DblType ReTurb = density.get(meshIdx, 0) * tke.get(meshIdx, 0) / + sdr.get(meshIdx, 0) / visc.get(meshIdx, 0); + + const DblType rK = 6.0; + const DblType a0Star = 0.072 / 3.0; + const DblType aStar = (a0Star + ReTurb / rK) / (1.0 + ReTurb / rK); + + const DblType omegaHat = stk::math::max( + sdr.get(meshIdx, 0), + (7.0 / 8.0) * (sijMag / stk::math::sqrt(0.09 / aStar))); + + tvisc.get(meshIdx, 0) = aStar * density.get(meshIdx, 0) * + tke.get(meshIdx, 0) / + stk::math::max(sdr.get(meshIdx, 0), 1.e-8); + }); + tvisc.modify_on_device(); +} + +} // namespace nalu +} // namespace sierra diff --git a/src/node_kernels/CMakeLists.txt b/src/node_kernels/CMakeLists.txt index 4e10cfbd25..cb93c27266 100644 --- a/src/node_kernels/CMakeLists.txt +++ b/src/node_kernels/CMakeLists.txt @@ -14,13 +14,17 @@ target_sources(nalu PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/ScalarGclNodeKernel.C ${CMAKE_CURRENT_SOURCE_DIR}/ScalarMassBDFNodeKernel.C ${CMAKE_CURRENT_SOURCE_DIR}/SDRSSTAMSNodeKernel.C + ${CMAKE_CURRENT_SOURCE_DIR}/TDRKENodeKernel.C ${CMAKE_CURRENT_SOURCE_DIR}/TKESSTAMSNodeKernel.C ${CMAKE_CURRENT_SOURCE_DIR}/TKERodiNodeKernel.C ${CMAKE_CURRENT_SOURCE_DIR}/TKEKsgsNodeKernel.C ${CMAKE_CURRENT_SOURCE_DIR}/TKESSTNodeKernel.C + ${CMAKE_CURRENT_SOURCE_DIR}/TKEKENodeKernel.C + ${CMAKE_CURRENT_SOURCE_DIR}/TKEKONodeKernel.C ${CMAKE_CURRENT_SOURCE_DIR}/TKESSTDESNodeKernel.C ${CMAKE_CURRENT_SOURCE_DIR}/TKESSTIDDESNodeKernel.C ${CMAKE_CURRENT_SOURCE_DIR}/SDRSSTNodeKernel.C + ${CMAKE_CURRENT_SOURCE_DIR}/SDRKONodeKernel.C ${CMAKE_CURRENT_SOURCE_DIR}/SDRSSTDESNodeKernel.C ${CMAKE_CURRENT_SOURCE_DIR}/WallDistNodeKernel.C ${CMAKE_CURRENT_SOURCE_DIR}/BLTGammaM2015NodeKernel.C diff --git a/src/node_kernels/SDRKONodeKernel.C b/src/node_kernels/SDRKONodeKernel.C new file mode 100644 index 0000000000..465612fc4e --- /dev/null +++ b/src/node_kernels/SDRKONodeKernel.C @@ -0,0 +1,142 @@ +// 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 "node_kernels/SDRKONodeKernel.h" +#include "Realm.h" +#include "SolutionOptions.h" +#include "SimdInterface.h" +#include "utils/StkHelpers.h" + +#include "stk_mesh/base/MetaData.hpp" +#include "stk_mesh/base/Types.hpp" + +namespace sierra { +namespace nalu { + +SDRKONodeKernel::SDRKONodeKernel(const stk::mesh::MetaData& meta) + : NGPNodeKernel(), + tkeID_(get_field_ordinal(meta, "turbulent_ke")), + sdrID_(get_field_ordinal(meta, "specific_dissipation_rate")), + densityID_(get_field_ordinal(meta, "density")), + tviscID_(get_field_ordinal(meta, "turbulent_viscosity")), + viscID_(get_field_ordinal(meta, "viscosity")), + dudxID_(get_field_ordinal(meta, "dudx")), + dkdxID_(get_field_ordinal(meta, "dkdx")), + dwdxID_(get_field_ordinal(meta, "dwdx")), + dualNodalVolumeID_(get_field_ordinal(meta, "dual_nodal_volume")), + nDim_(meta.spatial_dimension()) +{ +} + +void +SDRKONodeKernel::setup(Realm& realm) +{ + const auto& fieldMgr = realm.ngp_field_manager(); + + tke_ = fieldMgr.get_field(tkeID_); + sdr_ = fieldMgr.get_field(sdrID_); + density_ = fieldMgr.get_field(densityID_); + tvisc_ = fieldMgr.get_field(tviscID_); + visc_ = fieldMgr.get_field(viscID_); + dudx_ = fieldMgr.get_field(dudxID_); + dkdx_ = fieldMgr.get_field(dkdxID_); + dwdx_ = fieldMgr.get_field(dwdxID_); + dualNodalVolume_ = fieldMgr.get_field(dualNodalVolumeID_); + + // Update turbulence model constants + betaStar_ = realm.get_turb_model_constant(TM_betaStar); + tkeProdLimitRatio_ = realm.get_turb_model_constant(TM_tkeProdLimitRatio); + sigmaWTwo_ = realm.get_turb_model_constant(TM_sigmaWTwo); + betaOne_ = realm.get_turb_model_constant(TM_betaOne); + betaTwo_ = realm.get_turb_model_constant(TM_betaTwo); + gammaOne_ = realm.get_turb_model_constant(TM_gammaOne); + gammaTwo_ = realm.get_turb_model_constant(TM_gammaTwo); +} + +void +SDRKONodeKernel::execute( + NodeKernelTraits::LhsType& lhs, + NodeKernelTraits::RhsType& rhs, + const stk::mesh::FastMeshIndex& node) +{ + using DblType = NodeKernelTraits::DblType; + + const DblType tke = tke_.get(node, 0); + const DblType sdr = sdr_.get(node, 0); + const DblType density = density_.get(node, 0); + const DblType tvisc = tvisc_.get(node, 0); + const DblType visc = visc_.get(node, 0); + const DblType dVol = dualNodalVolume_.get(node, 0); + + DblType Pk = 0.0; + DblType crossDiff = 0.0; + for (int i = 0; i < nDim_; ++i) { + crossDiff += dkdx_.get(node, i) * dwdx_.get(node, i); + const int offset = nDim_ * i; + for (int j = 0; j < nDim_; ++j) { + const auto dudxij = dudx_.get(node, offset + j); + Pk += dudxij * (dudxij + dudx_.get(node, j * nDim_ + i)); + } + } + Pk *= tvisc; + + DblType chi_numer = 0.0; + for (int i = 0; i < nDim_; ++i) { + for (int j = 0; j < nDim_; ++j) { + for (int k = 0; k < nDim_; ++k) { + const auto rot_ij = 0.5 * (dudx_.get(node, i * nDim_ + j) - + dudx_.get(node, j * nDim_ + i)); + const auto rot_jk = 0.5 * (dudx_.get(node, j * nDim_ + k) - + dudx_.get(node, k * nDim_ + j)); + const auto str_ki = 0.5 * (dudx_.get(node, k * nDim_ + i) + + dudx_.get(node, i * nDim_ + k)); + chi_numer += rot_ij * rot_jk * str_ki; + } + } + } + + // JAM: Changes for SWH LowRe + const NodeKernelTraits::DblType alpha0_star = 0.072 / 3.0; + const NodeKernelTraits::DblType alpha_inf = 0.52; + const NodeKernelTraits::DblType alpha0 = 1.0 / 9.0; + const NodeKernelTraits::DblType Rk = 6.0; + const NodeKernelTraits::DblType R_beta = 8.0; + const NodeKernelTraits::DblType Rw = 2.95; + const NodeKernelTraits::DblType ReT = density * tke / sdr / visc; + const DblType Rbeta = 8.0; + const DblType betaStarLowRe = + betaStar_ * (4.0 / 15.0 + stk::math::pow(ReT / Rbeta, 4.0)) / + (1.0 + stk::math::pow(ReT / Rbeta, 4.0)); + DblType Dk = betaStarLowRe * density * sdr * tke; + + const DblType chi_omega = stk::math::abs( + chi_numer / stk::math::pow(0.09 * stk::math::max(sdr, 1.e-8), 3.0)); + const DblType beta = + 0.072 * (1.0 + 70.0 * chi_omega) / (1.0 + 80.0 * chi_omega); + + // JAM: Added for SWH LowRe + const NodeKernelTraits::DblType alpha_star = + (alpha0_star + ReT / Rk) / (1.0 + ReT / Rk); + const NodeKernelTraits::DblType alpha = + (alpha_inf / alpha_star) * ((alpha0 + ReT / Rw) / (1.0 + ReT / Rw)); + + // Pw includes 1/tvisc scaling; tvisc may be zero at a dirichlet low Re + // approach (clip) + // JAM: Changes for SWH LowRe, check densities... + const NodeKernelTraits::DblType Pw = + alpha * Pk * sdr / stk::math::max(tke, 1.e-12); + // Production term with appropriate clipping of tvisc + const DblType Dw = beta * density * sdr * sdr; + + rhs(0) += (Pw - Dw) * dVol; + lhs(0, 0) += 2.0 * beta * density * sdr * dVol; +} + +} // namespace nalu +} // namespace sierra diff --git a/src/node_kernels/SDRSSTDESNodeKernel.C b/src/node_kernels/SDRSSTDESNodeKernel.C index 136084dfba..76c8fa38dc 100644 --- a/src/node_kernels/SDRSSTDESNodeKernel.C +++ b/src/node_kernels/SDRSSTDESNodeKernel.C @@ -52,9 +52,6 @@ SDRSSTDESNodeKernel::setup(Realm& realm) fOneBlend_ = fieldMgr.get_field(fOneBlendID_); cellLengthScale_ = fieldMgr.get_field(cellLengthScaleID_); - const std::string dofName = "specific_dissipation_rate"; - relaxFac_ = realm.solutionOptions_->get_relaxation_factor(dofName); - // Update turbulence model constants betaStar_ = realm.get_turb_model_constant(TM_betaStar); tkeProdLimitRatio_ = realm.get_turb_model_constant(TM_tkeProdLimitRatio); diff --git a/src/node_kernels/SDRSSTNodeKernel.C b/src/node_kernels/SDRSSTNodeKernel.C index 974e6c2f8f..fe06ea632d 100644 --- a/src/node_kernels/SDRSSTNodeKernel.C +++ b/src/node_kernels/SDRSSTNodeKernel.C @@ -49,9 +49,6 @@ SDRSSTNodeKernel::setup(Realm& realm) dualNodalVolume_ = fieldMgr.get_field(dualNodalVolumeID_); fOneBlend_ = fieldMgr.get_field(fOneBlendID_); - const std::string dofName = "specific_dissipation_rate"; - relaxFac_ = realm.solutionOptions_->get_relaxation_factor(dofName); - // Update turbulence model constants betaStar_ = realm.get_turb_model_constant(TM_betaStar); tkeProdLimitRatio_ = realm.get_turb_model_constant(TM_tkeProdLimitRatio); diff --git a/src/node_kernels/TDRKENodeKernel.C b/src/node_kernels/TDRKENodeKernel.C new file mode 100644 index 0000000000..c48bd424b3 --- /dev/null +++ b/src/node_kernels/TDRKENodeKernel.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 "node_kernels/TDRKENodeKernel.h" +#include "Realm.h" +#include "SolutionOptions.h" +#include "SimdInterface.h" +#include "utils/StkHelpers.h" + +#include "stk_mesh/base/MetaData.hpp" +#include "stk_mesh/base/Types.hpp" + +namespace sierra { +namespace nalu { + +TDRKENodeKernel::TDRKENodeKernel(const stk::mesh::MetaData& meta) + : NGPNodeKernel(), + tkeID_(get_field_ordinal(meta, "turbulent_ke")), + tdrID_(get_field_ordinal(meta, "total_dissipation_rate")), + densityID_(get_field_ordinal(meta, "density")), + tviscID_(get_field_ordinal(meta, "turbulent_viscosity")), + viscID_(get_field_ordinal(meta, "viscosity")), + wallDistID_(get_field_ordinal(meta, "minimum_distance_to_wall")), + dplusID_(get_field_ordinal(meta, "dplus_wall_function")), + dudxID_(get_field_ordinal(meta, "dudx")), + dualNodalVolumeID_(get_field_ordinal(meta, "dual_nodal_volume")), + nDim_(meta.spatial_dimension()) +{ +} + +void +TDRKENodeKernel::setup(Realm& realm) +{ + const auto& fieldMgr = realm.ngp_field_manager(); + + tke_ = fieldMgr.get_field(tkeID_); + tdr_ = fieldMgr.get_field(tdrID_); + density_ = fieldMgr.get_field(densityID_); + tvisc_ = fieldMgr.get_field(tviscID_); + visc_ = fieldMgr.get_field(viscID_); + wallDist_ = fieldMgr.get_field(wallDistID_); + dplus_ = fieldMgr.get_field(dplusID_); + dudx_ = fieldMgr.get_field(dudxID_); + dualNodalVolume_ = fieldMgr.get_field(dualNodalVolumeID_); + + // Update turbulence model constants + cEpsOne_ = realm.get_turb_model_constant(TM_cEpsOne); + cEpsTwo_ = realm.get_turb_model_constant(TM_cEpsTwo); + fOne_ = realm.get_turb_model_constant(TM_fOne); +} + +void +TDRKENodeKernel::execute( + NodeKernelTraits::LhsType& lhs, + NodeKernelTraits::RhsType& rhs, + const stk::mesh::FastMeshIndex& node) +{ + using DblType = NodeKernelTraits::DblType; + + const DblType tke = tke_.get(node, 0); + const DblType tdr = tdr_.get(node, 0); + const DblType density = density_.get(node, 0); + const DblType tvisc = tvisc_.get(node, 0); + const DblType visc = visc_.get(node, 0); + const DblType wallDist = wallDist_.get(node, 0); + const DblType dplus = dplus_.get(node, 0); + const DblType dVol = dualNodalVolume_.get(node, 0); + + DblType Pk = 0.0; + for (int i = 0; i < nDim_; ++i) { + const int offset = nDim_ * i; + for (int j = 0; j < nDim_; ++j) { + const auto dudxij = dudx_.get(node, offset + j); + Pk += dudxij * (dudxij + dudx_.get(node, j * nDim_ + i)); + } + } + Pk *= tvisc; + + const double Re_t = density * tke * tke / visc / stk::math::max(tdr, 1.0e-16); + const double fTwo = 1.0 - 0.4 / 1.8 * stk::math::exp(-Re_t * Re_t / 36.0); + + // Production term with appropriate clipping of tvisc + const DblType Pe = cEpsOne_ * fOne_ * Pk * tdr / stk::math::max(tke, 1.0e-16); + const DblType DeFac = + cEpsTwo_ * fTwo * density * tdr / stk::math::max(tke, 1.0e-16); + const DblType De = DeFac * tdr; + const DblType LeFac = + 2.0 * visc * stk::math::exp(-0.5 * dplus) / wallDist / wallDist; + const DblType Le = -LeFac * tdr; + + rhs(0) += (Pe - De + Le) * dVol; + lhs(0, 0) += (2.0 * DeFac + LeFac) * dVol; +} + +} // namespace nalu +} // namespace sierra diff --git a/src/node_kernels/TKEKENodeKernel.C b/src/node_kernels/TKEKENodeKernel.C new file mode 100644 index 0000000000..510f737681 --- /dev/null +++ b/src/node_kernels/TKEKENodeKernel.C @@ -0,0 +1,93 @@ +// 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 "node_kernels/TKEKENodeKernel.h" +#include "Realm.h" +#include "SolutionOptions.h" +#include "SimdInterface.h" +#include "utils/StkHelpers.h" + +#include "stk_mesh/base/MetaData.hpp" +#include "stk_mesh/base/Types.hpp" + +namespace sierra { +namespace nalu { + +TKEKENodeKernel::TKEKENodeKernel(const stk::mesh::MetaData& meta) + : NGPNodeKernel(), + tkeID_(get_field_ordinal(meta, "turbulent_ke")), + tdrID_(get_field_ordinal(meta, "total_dissipation_rate")), + densityID_(get_field_ordinal(meta, "density")), + tviscID_(get_field_ordinal(meta, "turbulent_viscosity")), + viscID_(get_field_ordinal(meta, "viscosity")), + dudxID_(get_field_ordinal(meta, "dudx")), + wallDistID_(get_field_ordinal(meta, "minimum_distance_to_wall")), + dualNodalVolumeID_(get_field_ordinal(meta, "dual_nodal_volume")), + nDim_(meta.spatial_dimension()) +{ +} + +void +TKEKENodeKernel::setup(Realm& realm) +{ + const auto& fieldMgr = realm.ngp_field_manager(); + + tke_ = fieldMgr.get_field(tkeID_); + tdr_ = fieldMgr.get_field(tdrID_); + density_ = fieldMgr.get_field(densityID_); + visc_ = fieldMgr.get_field(viscID_); + tvisc_ = fieldMgr.get_field(tviscID_); + dudx_ = fieldMgr.get_field(dudxID_); + wallDist_ = fieldMgr.get_field(wallDistID_); + dualNodalVolume_ = fieldMgr.get_field(dualNodalVolumeID_); + + // Update turbulence model constants + betaStar_ = realm.get_turb_model_constant(TM_betaStar); + tkeProdLimitRatio_ = realm.get_turb_model_constant(TM_tkeProdLimitRatio); +} + +void +TKEKENodeKernel::execute( + NodeKernelTraits::LhsType& lhs, + NodeKernelTraits::RhsType& rhs, + const stk::mesh::FastMeshIndex& node) +{ + using DblType = NodeKernelTraits::DblType; + + // See https://turbmodels.larc.nasa.gov/sst.html for details + + const DblType tke = tke_.get(node, 0); + const DblType tdr = tdr_.get(node, 0); + const DblType density = density_.get(node, 0); + const DblType tvisc = tvisc_.get(node, 0); + const DblType visc = visc_.get(node, 0); + const DblType wallDist = wallDist_.get(node, 0); + const DblType dVol = dualNodalVolume_.get(node, 0); + + DblType Pk = 0.0; + for (int i = 0; i < nDim_; ++i) { + const int offset = nDim_ * i; + for (int j = 0; j < nDim_; ++j) { + const auto dudxij = dudx_.get(node, offset + j); + Pk += dudxij * (dudxij + dudx_.get(node, j * nDim_ + i)); + } + } + Pk *= tvisc; + + DblType Dk = density * tdr; + + const DblType lFac = 2.0 * visc / wallDist / wallDist; + DblType Lk = -lFac * tke; + + rhs(0) += (Pk - Dk + Lk) * dVol; + lhs(0, 0) += lFac * dVol; +} + +} // namespace nalu +} // namespace sierra diff --git a/src/node_kernels/TKEKONodeKernel.C b/src/node_kernels/TKEKONodeKernel.C new file mode 100644 index 0000000000..f214d8e7f3 --- /dev/null +++ b/src/node_kernels/TKEKONodeKernel.C @@ -0,0 +1,91 @@ +// 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 "node_kernels/TKEKONodeKernel.h" +#include "Realm.h" +#include "SolutionOptions.h" +#include "SimdInterface.h" +#include "utils/StkHelpers.h" + +#include "stk_mesh/base/MetaData.hpp" +#include "stk_mesh/base/Types.hpp" + +namespace sierra { +namespace nalu { + +TKEKONodeKernel::TKEKONodeKernel(const stk::mesh::MetaData& meta) + : NGPNodeKernel(), + tkeID_(get_field_ordinal(meta, "turbulent_ke")), + sdrID_(get_field_ordinal(meta, "specific_dissipation_rate")), + densityID_(get_field_ordinal(meta, "density")), + tviscID_(get_field_ordinal(meta, "turbulent_viscosity")), + viscID_(get_field_ordinal(meta, "viscosity")), + dudxID_(get_field_ordinal(meta, "dudx")), + dualNodalVolumeID_(get_field_ordinal(meta, "dual_nodal_volume")), + nDim_(meta.spatial_dimension()) +{ +} + +void +TKEKONodeKernel::setup(Realm& realm) +{ + const auto& fieldMgr = realm.ngp_field_manager(); + + tke_ = fieldMgr.get_field(tkeID_); + sdr_ = fieldMgr.get_field(sdrID_); + density_ = fieldMgr.get_field(densityID_); + tvisc_ = fieldMgr.get_field(tviscID_); + visc_ = fieldMgr.get_field(viscID_); + dudx_ = fieldMgr.get_field(dudxID_); + dualNodalVolume_ = fieldMgr.get_field(dualNodalVolumeID_); + + // Update turbulence model constants + betaStar_ = realm.get_turb_model_constant(TM_betaStar); + tkeProdLimitRatio_ = realm.get_turb_model_constant(TM_tkeProdLimitRatio); +} + +void +TKEKONodeKernel::execute( + NodeKernelTraits::LhsType& lhs, + NodeKernelTraits::RhsType& rhs, + const stk::mesh::FastMeshIndex& node) +{ + using DblType = NodeKernelTraits::DblType; + + const DblType tke = tke_.get(node, 0); + const DblType sdr = sdr_.get(node, 0); + const DblType density = density_.get(node, 0); + const DblType tvisc = tvisc_.get(node, 0); + const DblType visc = visc_.get(node, 0); + const DblType dVol = dualNodalVolume_.get(node, 0); + + DblType Pk = 0.0; + for (int i = 0; i < nDim_; ++i) { + const int offset = nDim_ * i; + for (int j = 0; j < nDim_; ++j) { + const auto dudxij = dudx_.get(node, offset + j); + Pk += dudxij * (dudxij + dudx_.get(node, j * nDim_ + i)); + } + } + Pk *= tvisc; + + // JAM: Changes for SWH LowRe + const DblType ReT = density * tke / sdr / visc; + const DblType Rbeta = 8.0; + const DblType betaStarLowRe = + betaStar_ * (4.0 / 15.0 + stk::math::pow(ReT / Rbeta, 4.0)) / + (1.0 + stk::math::pow(ReT / Rbeta, 4.0)); + DblType Dk = betaStarLowRe * density * sdr * tke; + + rhs(0) += (Pk - Dk) * dVol; + lhs(0, 0) += betaStarLowRe * density * sdr * dVol; +} + +} // namespace nalu +} // namespace sierra diff --git a/src/node_kernels/TKEKsgsNodeKernel.C b/src/node_kernels/TKEKsgsNodeKernel.C index fe2c252912..b8e044685a 100644 --- a/src/node_kernels/TKEKsgsNodeKernel.C +++ b/src/node_kernels/TKEKsgsNodeKernel.C @@ -42,9 +42,6 @@ TKEKsgsNodeKernel::setup(Realm& realm) dudx_ = fieldMgr.get_field(dudxID_); dualNodalVolume_ = fieldMgr.get_field(dualNodalVolumeID_); - const std::string dofName = "turbulent_ke"; - relaxFac_ = realm.solutionOptions_->get_relaxation_factor(dofName); - // Update turbulence model constants cEps_ = realm.get_turb_model_constant(TM_cEps); tkeProdLimitRatio_ = realm.get_turb_model_constant(TM_tkeProdLimitRatio); diff --git a/src/node_kernels/TKESSTDESNodeKernel.C b/src/node_kernels/TKESSTDESNodeKernel.C index 491d84842e..0023ffc9e3 100644 --- a/src/node_kernels/TKESSTDESNodeKernel.C +++ b/src/node_kernels/TKESSTDESNodeKernel.C @@ -48,9 +48,6 @@ TKESSTDESNodeKernel::setup(Realm& realm) maxLenScale_ = fieldMgr.get_field(maxLenScaleID_); fOneBlend_ = fieldMgr.get_field(fOneBlendID_); - const std::string dofName = "turbulent_ke"; - relaxFac_ = realm.solutionOptions_->get_relaxation_factor(dofName); - // Update turbulence model constants betaStar_ = realm.get_turb_model_constant(TM_betaStar); tkeProdLimitRatio_ = realm.get_turb_model_constant(TM_tkeProdLimitRatio); diff --git a/src/node_kernels/TKESSTIDDESNodeKernel.C b/src/node_kernels/TKESSTIDDESNodeKernel.C index 80d1fa5508..2a183c1a46 100644 --- a/src/node_kernels/TKESSTIDDESNodeKernel.C +++ b/src/node_kernels/TKESSTIDDESNodeKernel.C @@ -49,9 +49,6 @@ TKESSTIDDESNodeKernel::setup(Realm& realm) fOneBlend_ = fieldMgr.get_field(fOneBlendID_); ransIndicator_ = fieldMgr.get_field(ransIndicatorID_); - const std::string dofName = "turbulent_ke"; - relaxFac_ = realm.solutionOptions_->get_relaxation_factor(dofName); - // Update turbulence model constants betaStar_ = realm.get_turb_model_constant(TM_betaStar); tkeProdLimitRatio_ = realm.get_turb_model_constant(TM_tkeProdLimitRatio); diff --git a/src/node_kernels/TKESSTNodeKernel.C b/src/node_kernels/TKESSTNodeKernel.C index 7b79719431..4437287482 100644 --- a/src/node_kernels/TKESSTNodeKernel.C +++ b/src/node_kernels/TKESSTNodeKernel.C @@ -43,9 +43,6 @@ TKESSTNodeKernel::setup(Realm& realm) dudx_ = fieldMgr.get_field(dudxID_); dualNodalVolume_ = fieldMgr.get_field(dualNodalVolumeID_); - const std::string dofName = "turbulent_ke"; - relaxFac_ = realm.solutionOptions_->get_relaxation_factor(dofName); - // Update turbulence model constants betaStar_ = realm.get_turb_model_constant(TM_betaStar); tkeProdLimitRatio_ = realm.get_turb_model_constant(TM_tkeProdLimitRatio); diff --git a/unit_tests/kernels/UnitTestKernelUtils.C b/unit_tests/kernels/UnitTestKernelUtils.C index c35f2533c4..97bfd16bda 100644 --- a/unit_tests/kernels/UnitTestKernelUtils.C +++ b/unit_tests/kernels/UnitTestKernelUtils.C @@ -168,6 +168,18 @@ struct TrigFieldFunction std::sin(a * pi * z)); } + void tdr(const double* coords, double* qField) const + { + double x = coords[0]; + double y = coords[1]; + double z = coords[2]; + + qField[0] = 2*tdrnot + tdrnot * ( + std::cos(a * pi * x) * + std::sin(a * pi * y) * + std::sin(a * pi * z)); + } + void dwdx(const double* coords, double* qField) const { const double x = coords[0]; @@ -235,6 +247,12 @@ struct TrigFieldFunction qField[0] = 10*x+10; } + void dplus_wall_function(const double* coords, double* qField) const + { + double x = coords[0]; + qField[0] = 2*x+2; + } + void dhdx(const double* /*coords*/, double* qField) const { qField[0] = 30.0; @@ -273,6 +291,9 @@ private: /// Factor for sdr field static constexpr double sdrnot{1.0}; + /// Factor for tdr field + static constexpr double tdrnot{1.0}; + /// Factor for tvisc field static constexpr double tviscnot{1.0}; @@ -318,6 +339,8 @@ void init_trigonometric_field( funcPtr = &TrigFieldFunction::dkdx; else if (fieldName == "specific_dissipation_rate") funcPtr = &TrigFieldFunction::sdr; + else if (fieldName == "total_dissipation_rate") + funcPtr = &TrigFieldFunction::tdr; else if (fieldName == "dwdx") funcPtr = &TrigFieldFunction::dwdx; else if (fieldName == "dhdx") @@ -330,6 +353,8 @@ void init_trigonometric_field( funcPtr = &TrigFieldFunction::sst_f_one_blending; else if (fieldName == "minimum_distance_to_wall") funcPtr = &TrigFieldFunction::minimum_distance_to_wall; + else if (fieldName == "dplus_wall_function") + funcPtr = &TrigFieldFunction::dplus_wall_function; else funcPtr = nullptr; @@ -461,6 +486,14 @@ void sdr_test_function( init_trigonometric_field(bulk, coordinates, sdr); } +void tdr_test_function( + const stk::mesh::BulkData& bulk, + const VectorFieldType& coordinates, + ScalarFieldType& tdr) +{ + init_trigonometric_field(bulk, coordinates, tdr); +} + void dwdx_test_function( const stk::mesh::BulkData& bulk, const VectorFieldType& coordinates, @@ -500,6 +533,14 @@ void minimum_distance_to_wall_test_function( init_trigonometric_field(bulk, coordinates, minimum_distance_to_wall); } +void dplus_test_function( + const stk::mesh::BulkData& bulk, + const VectorFieldType& coordinates, + ScalarFieldType& dplus) +{ + init_trigonometric_field(bulk, coordinates, dplus); +} + void property_from_mixture_fraction_test_function( const stk::mesh::BulkData& bulk, const ScalarFieldType& mixFraction, diff --git a/unit_tests/kernels/UnitTestKernelUtils.h b/unit_tests/kernels/UnitTestKernelUtils.h index 28affbd18f..76effa7689 100644 --- a/unit_tests/kernels/UnitTestKernelUtils.h +++ b/unit_tests/kernels/UnitTestKernelUtils.h @@ -85,6 +85,11 @@ void sdr_test_function( const VectorFieldType& coordinates, ScalarFieldType& sdr); +void tdr_test_function( + const stk::mesh::BulkData& bulk, + const VectorFieldType& coordinates, + ScalarFieldType& tdr); + void dwdx_test_function( const stk::mesh::BulkData& bulk, const VectorFieldType& coordinates, @@ -110,6 +115,11 @@ void minimum_distance_to_wall_test_function( const VectorFieldType& coordinates, ScalarFieldType& minimum_distance_to_wall); +void dplus_test_function( + const stk::mesh::BulkData& bulk, + const VectorFieldType& coordinates, + ScalarFieldType& dplus); + void property_from_mixture_fraction_test_function( const stk::mesh::BulkData&, const ScalarFieldType& mixFraction, @@ -771,6 +781,131 @@ class SSTKernelHex8Mesh : public LowMachKernelHex8Mesh ScalarFieldType* pecletFactor_ {nullptr}; }; +/** Test Fixture for the KE Kernels + * + */ +class KEKernelHex8Mesh : public LowMachKernelHex8Mesh +{ +public: + KEKernelHex8Mesh() + : LowMachKernelHex8Mesh(), + tke_(&meta_->declare_field( + stk::topology::NODE_RANK, "turbulent_ke")), + tdr_(&meta_->declare_field( + stk::topology::NODE_RANK, "total_dissipation_rate")), + visc_(&meta_->declare_field( + stk::topology::NODE_RANK, "viscosity")), + tvisc_(&meta_->declare_field( + stk::topology::NODE_RANK, "turbulent_viscosity")), + minDistance_(&meta_->declare_field( + stk::topology::NODE_RANK, "minimum_distance_to_wall")), + dplus_(&meta_->declare_field( + stk::topology::NODE_RANK, "dplus_wall_function")), + dudx_(&meta_->declare_field( + stk::topology::NODE_RANK, "dudx")) + { + stk::mesh::put_field_on_mesh(*tke_, meta_->universal_part(), 1, nullptr); + stk::mesh::put_field_on_mesh(*tdr_, meta_->universal_part(), 1, nullptr); + stk::mesh::put_field_on_mesh(*visc_, meta_->universal_part(), 1, nullptr); + stk::mesh::put_field_on_mesh(*tvisc_, meta_->universal_part(), 1, nullptr); + stk::mesh::put_field_on_mesh(*minDistance_, meta_->universal_part(), 1, nullptr); + stk::mesh::put_field_on_mesh(*dplus_, meta_->universal_part(), 1, nullptr); + stk::mesh::put_field_on_mesh( + *dudx_, meta_->universal_part(), spatialDim_ * spatialDim_, nullptr); + } + + virtual ~KEKernelHex8Mesh() {} + + virtual void fill_mesh_and_init_fields( + const bool doPerturb = false, const bool generateSidesets = false) override + { + LowMachKernelHex8Mesh::fill_mesh_and_init_fields(doPerturb, generateSidesets); + stk::mesh::field_fill(0.2, *visc_); + stk::mesh::field_fill(0.3, *tvisc_); + unit_test_kernel_utils::density_test_function( + *bulk_, *coordinates_, *density_); + unit_test_kernel_utils::tke_test_function(*bulk_, *coordinates_, *tke_); + unit_test_kernel_utils::tdr_test_function(*bulk_, *coordinates_, *tdr_); + unit_test_kernel_utils::minimum_distance_to_wall_test_function(*bulk_, *coordinates_, *minDistance_); + unit_test_kernel_utils::dplus_test_function(*bulk_, *coordinates_, *dplus_); + unit_test_kernel_utils::dudx_test_function(*bulk_, *coordinates_, *dudx_); + } + + ScalarFieldType* tke_{nullptr}; + ScalarFieldType* tdr_{nullptr}; + ScalarFieldType* visc_{nullptr}; + ScalarFieldType* tvisc_{nullptr}; + ScalarFieldType* minDistance_{nullptr}; + ScalarFieldType* dplus_{nullptr}; + GenericFieldType* dudx_{nullptr}; +}; + +/** Test Fixture for the KO Kernels + * + */ +class KOKernelHex8Mesh : public LowMachKernelHex8Mesh +{ +public: + KOKernelHex8Mesh() + : LowMachKernelHex8Mesh(), + tke_(&meta_->declare_field( + stk::topology::NODE_RANK, "turbulent_ke")), + sdr_(&meta_->declare_field( + stk::topology::NODE_RANK, "specific_dissipation_rate")), + visc_(&meta_->declare_field( + stk::topology::NODE_RANK, "viscosity")), + tvisc_(&meta_->declare_field( + stk::topology::NODE_RANK, "turbulent_viscosity")), + minDistance_(&meta_->declare_field( + stk::topology::NODE_RANK, "minimum_distance_to_wall")), + dudx_(&meta_->declare_field( + stk::topology::NODE_RANK, "dudx")), + dkdx_(&meta_->declare_field( + stk::topology::NODE_RANK, "dkdx")), + dwdx_(&meta_->declare_field( + stk::topology::NODE_RANK, "dwdx")) + { + stk::mesh::put_field_on_mesh(*tke_, meta_->universal_part(), 1, nullptr); + stk::mesh::put_field_on_mesh(*sdr_, meta_->universal_part(), 1, nullptr); + stk::mesh::put_field_on_mesh(*visc_, meta_->universal_part(), 1, nullptr); + stk::mesh::put_field_on_mesh(*tvisc_, meta_->universal_part(), 1, nullptr); + stk::mesh::put_field_on_mesh(*minDistance_, meta_->universal_part(), 1, nullptr); + stk::mesh::put_field_on_mesh( + *dudx_, meta_->universal_part(), spatialDim_ * spatialDim_, nullptr); + stk::mesh::put_field_on_mesh( + *dkdx_, meta_->universal_part(), spatialDim_, nullptr); + stk::mesh::put_field_on_mesh( + *dwdx_, meta_->universal_part(), spatialDim_, nullptr); + } + + virtual ~KOKernelHex8Mesh() {} + + virtual void fill_mesh_and_init_fields( + const bool doPerturb = false, const bool generateSidesets = false) override + { + LowMachKernelHex8Mesh::fill_mesh_and_init_fields(doPerturb, generateSidesets); + stk::mesh::field_fill(0.2, *visc_); + stk::mesh::field_fill(0.3, *tvisc_); + unit_test_kernel_utils::density_test_function( + *bulk_, *coordinates_, *density_); + unit_test_kernel_utils::tke_test_function(*bulk_, *coordinates_, *tke_); + unit_test_kernel_utils::sdr_test_function(*bulk_, *coordinates_, *sdr_); + unit_test_kernel_utils::minimum_distance_to_wall_test_function(*bulk_, *coordinates_, *minDistance_); + unit_test_kernel_utils::dudx_test_function(*bulk_, *coordinates_, *dudx_); + stk::mesh::field_fill(0.0, *dkdx_); + stk::mesh::field_fill(0.0, *dwdx_); + } + + ScalarFieldType* tke_{nullptr}; + ScalarFieldType* sdr_{nullptr}; + ScalarFieldType* visc_{nullptr}; + ScalarFieldType* tvisc_{nullptr}; + ScalarFieldType* minDistance_{nullptr}; + GenericFieldType* dudx_{nullptr}; + VectorFieldType* dkdx_{nullptr}; + VectorFieldType* dwdx_{nullptr}; +}; + /** Test Fixture for the Turbulence Kernels * */ diff --git a/unit_tests/ngp_algorithms/CMakeLists.txt b/unit_tests/ngp_algorithms/CMakeLists.txt index a837273da7..dc15568992 100644 --- a/unit_tests/ngp_algorithms/CMakeLists.txt +++ b/unit_tests/ngp_algorithms/CMakeLists.txt @@ -8,6 +8,8 @@ target_sources(${utest_ex_name} PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestMdotAlg.C ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestTurbViscKsgsAlg.C ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestTurbViscSSTAlg.C + ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestTurbViscKEAlg.C + ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestTurbViscKOAlg.C ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestGeometryAlg.C ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestSDRWallAlg.C ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestNodalGradPOpenBoundary.C diff --git a/unit_tests/ngp_algorithms/UnitTestTurbViscKEAlg.C b/unit_tests/ngp_algorithms/UnitTestTurbViscKEAlg.C new file mode 100644 index 0000000000..d0faf4ee9e --- /dev/null +++ b/unit_tests/ngp_algorithms/UnitTestTurbViscKEAlg.C @@ -0,0 +1,63 @@ +// 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 "kernels/UnitTestKernelUtils.h" +#include "UnitTestHelperObjects.h" + +#include "ngp_algorithms/TurbViscKEAlg.h" + +TEST_F(KEKernelHex8Mesh, NGP_turb_visc_ke_alg) +{ + // Only execute for 1 processor runs + if (bulk_->parallel_size() > 1) return; + + KEKernelHex8Mesh::fill_mesh_and_init_fields(); + + // Initialize turbulence parameters in solution options + solnOpts_.initialize_turbulence_constants(); + + unit_test_utils::HelperObjects helperObjs( + bulk_, stk::topology::HEX_8, 1, partVec_[0]); + + sierra::nalu::TurbViscKEAlg TurbViscKEAlg( + helperObjs.realm, partVec_[0], tvisc_); + + TurbViscKEAlg.execute(); + + const auto& fieldMgr = helperObjs.realm.mesh_info().ngp_field_manager(); + auto ngpTvisc = fieldMgr.get_field(tvisc_->mesh_meta_data_ordinal()); + ngpTvisc.modify_on_device(); + ngpTvisc.sync_to_host(); + + { + std::vector expectedValues = { + 0.0040927529208101276, + 0.004756620887334956, + 0.0047455106720311118, + 0.0042834431659046351, + 0.0024056598081291371, + 0.0027958716083218245, + 0.0016322044534763486, + 0.0017904108086812805, + }; + + const double tol = 1.0e-15; + + stk::mesh::Selector sel = meta_->universal_part(); + const auto& bkts = bulk_->get_buckets(stk::topology::NODE_RANK, sel); + + int ii = 0; + for (const auto* b: bkts) + for (const auto node: *b) { + const double* tvisc = stk::mesh::field_data(*tvisc_, node); + EXPECT_NEAR(tvisc[0], expectedValues[ii++], tol); + } + } +} diff --git a/unit_tests/ngp_algorithms/UnitTestTurbViscKOAlg.C b/unit_tests/ngp_algorithms/UnitTestTurbViscKOAlg.C new file mode 100644 index 0000000000..2107f905f3 --- /dev/null +++ b/unit_tests/ngp_algorithms/UnitTestTurbViscKOAlg.C @@ -0,0 +1,63 @@ +// 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 "kernels/UnitTestKernelUtils.h" +#include "UnitTestHelperObjects.h" + +#include "ngp_algorithms/TurbViscKOAlg.h" + +TEST_F(KOKernelHex8Mesh, NGP_turb_visc_ko_alg) +{ + // Only execute for 1 processor runs + if (bulk_->parallel_size() > 1) return; + + KOKernelHex8Mesh::fill_mesh_and_init_fields(); + + // Initialize turbulence parameters in solution options + solnOpts_.initialize_turbulence_constants(); + + unit_test_utils::HelperObjects helperObjs( + bulk_, stk::topology::HEX_8, 1, partVec_[0]); + + sierra::nalu::TurbViscKOAlg TurbViscKOAlg( + helperObjs.realm, partVec_[0], tvisc_); + + TurbViscKOAlg.execute(); + + const auto& fieldMgr = helperObjs.realm.mesh_info().ngp_field_manager(); + auto ngpTvisc = fieldMgr.get_field(tvisc_->mesh_meta_data_ordinal()); + ngpTvisc.modify_on_device(); + ngpTvisc.sync_to_host(); + + { + std::vector expectedValues = { + 0.46763636363636363, + 0.20271993944117145, + 0.34820558301166943, + 0.11992191196521775, + 0.20271993944117145, + 0.083672109203660194, + 0.074293947075371264, + 0.031038745209797804, + }; + + const double tol = 1.0e-15; + + stk::mesh::Selector sel = meta_->universal_part(); + const auto& bkts = bulk_->get_buckets(stk::topology::NODE_RANK, sel); + + int ii = 0; + for (const auto* b: bkts) + for (const auto node: *b) { + const double* tvisc = stk::mesh::field_data(*tvisc_, node); + EXPECT_NEAR(tvisc[0], expectedValues[ii++], tol); + } + } +} diff --git a/unit_tests/node_kernels/CMakeLists.txt b/unit_tests/node_kernels/CMakeLists.txt index 961e97cc5e..fc75379603 100644 --- a/unit_tests/node_kernels/CMakeLists.txt +++ b/unit_tests/node_kernels/CMakeLists.txt @@ -15,6 +15,8 @@ target_sources(${utest_ex_name} PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestAMSNodeKernel.C ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestTKERodiNodeKernel.C ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestSSTNodeKernel.C + ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestKENodeKernel.C + ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestKONodeKernel.C ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestKsgsNodeKernel.C ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestWallDistNode.C ) diff --git a/unit_tests/node_kernels/UnitTestKENodeKernel.C b/unit_tests/node_kernels/UnitTestKENodeKernel.C new file mode 100644 index 0000000000..1a0a2b5a79 --- /dev/null +++ b/unit_tests/node_kernels/UnitTestKENodeKernel.C @@ -0,0 +1,124 @@ +// 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 "kernels/UnitTestKernelUtils.h" +#include "UnitTestUtils.h" +#include "UnitTestHelperObjects.h" + +#include "node_kernels/TKEKENodeKernel.h" +#include "node_kernels/TDRKENodeKernel.h" + +namespace { +namespace hex8_golds { +namespace tke_ke { +static constexpr double rhs[8] = { + -0.251, -0.14719631307311828, + -0.14835082157030577, -0.029604928511851453, + -0.14794631307311829, -0.086622875703131569, + -0.1158765303693139, -0.0037418778749914866, +}; + +static constexpr double lhs[8][8] = { + {0.0005, 0, 0, 0, 0, 0, 0, 0, }, + {0, 0.000125, 0, 0, 0, 0, 0, 0, }, + {0, 0, 0.0005, 0, 0, 0, 0, 0, }, + {0, 0, 0, 0.000125, 0, 0, 0, 0, }, + {0, 0, 0, 0, 0.0005, 0, 0, 0, }, + {0, 0, 0, 0, 0, 0.000125, 0, 0, }, + {0, 0, 0, 0, 0, 0, 0.0005, 0, }, + {0, 0, 0, 0, 0, 0, 0, 0.000125, }, +}; +} // tke_ke + +namespace tdr_ke { +static constexpr double rhs[8] = { + -0.44415022703895984, -0.24202453432573701, + -0.18769279143193071, -0.050569760373841061, + -0.2423585799460993, -0.13070571477449358, + -0.19014590536955761, -0.011309382975531565, +}; + +static constexpr double lhs[8][8] = { + {0.44396628731837412, 0, 0, 0, 0, 0, 0, 0, }, + {0, 0.24200761741533242, 0, 0, 0, 0, 0, 0, }, + {0, 0, 0.18750885171134499, 0, 0, 0, 0, 0, }, + {0, 0, 0, 0.11280579803843291, 0, 0, 0, 0, }, + {0, 0, 0, 0, 0.24217464022551358, 0, 0, 0, }, + {0, 0, 0, 0, 0, 0.13068879786408899, 0, 0, }, + {0, 0, 0, 0, 0, 0, 0.1430786684579321, 0, }, + {0, 0, 0, 0, 0, 0, 0, 0.0770742118579643,}, +}; +} // tdr_ke +} // hex8_golds +} + +TEST_F(KEKernelHex8Mesh, NGP_tke_ke_node) +{ + // Only execute for 1 processor runs + if (bulk_->parallel_size() > 1) return; + + fill_mesh_and_init_fields(); + + // Setup solution options + solnOpts_.meshMotion_ = false; + solnOpts_.externalMeshDeformation_ = false; + solnOpts_.initialize_turbulence_constants(); + + unit_test_utils::NodeHelperObjects helperObjs( + bulk_, stk::topology::HEX_8, 1, partVec_[0]); + + helperObjs.nodeAlg->add_kernel(*meta_); + + helperObjs.execute(); + + Kokkos::deep_copy(helperObjs.linsys->hostNumSumIntoCalls_, helperObjs.linsys->numSumIntoCalls_); + EXPECT_EQ(helperObjs.linsys->lhs_.extent(0), 8u); + EXPECT_EQ(helperObjs.linsys->lhs_.extent(1), 8u); + EXPECT_EQ(helperObjs.linsys->rhs_.extent(0), 8u); + EXPECT_EQ(helperObjs.linsys->hostNumSumIntoCalls_(0), 8u); + + namespace hex8_golds = hex8_golds::tke_ke; + unit_test_kernel_utils::expect_all_near( + helperObjs.linsys->rhs_, hex8_golds::rhs, 1.0e-12); + unit_test_kernel_utils::expect_all_near<8>( + helperObjs.linsys->lhs_, hex8_golds::lhs, 1.0e-12); +} + +TEST_F(KEKernelHex8Mesh, NGP_tdr_ke_node) +{ + // Only execute for 1 processor runs + if (bulk_->parallel_size() > 1) return; + + fill_mesh_and_init_fields(); + + // Setup solution options + solnOpts_.meshMotion_ = false; + solnOpts_.externalMeshDeformation_ = false; + solnOpts_.initialize_turbulence_constants(); + + unit_test_utils::NodeHelperObjects helperObjs( + bulk_, stk::topology::HEX_8, 1, partVec_[0]); + + helperObjs.nodeAlg->add_kernel(*meta_); + + helperObjs.execute(); + + Kokkos::deep_copy(helperObjs.linsys->hostNumSumIntoCalls_, helperObjs.linsys->numSumIntoCalls_); + EXPECT_EQ(helperObjs.linsys->lhs_.extent(0), 8u); + EXPECT_EQ(helperObjs.linsys->lhs_.extent(1), 8u); + EXPECT_EQ(helperObjs.linsys->rhs_.extent(0), 8u); + EXPECT_EQ(helperObjs.linsys->hostNumSumIntoCalls_(0), 8u); + + namespace hex8_golds = hex8_golds::tdr_ke; + unit_test_kernel_utils::expect_all_near( + helperObjs.linsys->rhs_, hex8_golds::rhs, 1.0e-12); + unit_test_kernel_utils::expect_all_near<8>( + helperObjs.linsys->lhs_, hex8_golds::lhs, 1.0e-12); +} diff --git a/unit_tests/node_kernels/UnitTestKONodeKernel.C b/unit_tests/node_kernels/UnitTestKONodeKernel.C new file mode 100644 index 0000000000..f4b2e615e4 --- /dev/null +++ b/unit_tests/node_kernels/UnitTestKONodeKernel.C @@ -0,0 +1,124 @@ +// 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 "kernels/UnitTestKernelUtils.h" +#include "UnitTestUtils.h" +#include "UnitTestHelperObjects.h" + +#include "node_kernels/TKEKONodeKernel.h" +#include "node_kernels/SDRKONodeKernel.h" + +namespace { +namespace hex8_golds { +namespace tke_ko { +static constexpr double rhs[8] = { + -0.016368777801313281, -0.00740039033812396, + -0.011709654168113211, 0.051874097062954767, + -0.00740039033812396, -0.00417063108566584, + -0.0068417453675951607, 0.053763691195545707, +}; + +static constexpr double lhs[8][8] = { +{0.0081843889006566403, 0, 0, 0, 0, 0, 0, 0, }, + {0, 0.00370019516906198, 0, 0, 0, 0, 0, 0, }, + {0, 0, 0.0041685949894791578, 0, 0, 0, 0, 0, }, + {0, 0, 0, 0.0021018912401700412, 0, 0, 0, 0, }, + {0, 0, 0, 0, 0.00370019516906198, 0, 0, 0, }, + {0, 0, 0, 0, 0, 0.00208531554283292, 0, 0, }, + {0, 0, 0, 0, 0, 0, 0.0027637516740426134, 0, }, + {0, 0, 0, 0, 0, 0, 0, 0.0014536892633176806, }, +}; +} // tke_ko + +namespace sdr_ko { +static constexpr double rhs[8] = { + -0.035999999999999997, -0.021160269082529031, + -0.021160269082529031, 0.029003272915464909, + -0.021160269082529031, -0.012437694101250944, + -0.021910289694610025, 0.053914090927629096, +}; + +static constexpr double lhs[8][8] = { + {0.035999999999999997, 0, 0, 0, 0, 0, 0, 0, }, + {0, 0.021160269082529031, 0, 0, 0, 0, 0, 0, }, + {0, 0, 0.021160269082529031, 0, 0, 0, 0, 0, }, + {0, 0, 0, 0.012437694101250944, 0, 0, 0, 0, }, + {0, 0, 0, 0, 0.021160269082529031, 0, 0, 0, }, + {0, 0, 0, 0, 0, 0.012437694101250944, 0, 0, }, + {0, 0, 0, 0, 0, 0, 0.016507982338594577, 0, }, + {0, 0, 0, 0, 0, 0, 0, 0.0087169431652403921,}, +}; +} // sdr_ko +} // hex8_golds +} + +TEST_F(KOKernelHex8Mesh, NGP_tke_ko_node) +{ + // Only execute for 1 processor runs + if (bulk_->parallel_size() > 1) return; + + fill_mesh_and_init_fields(); + + // Setup solution options + solnOpts_.meshMotion_ = false; + solnOpts_.externalMeshDeformation_ = false; + solnOpts_.initialize_turbulence_constants(); + + unit_test_utils::NodeHelperObjects helperObjs( + bulk_, stk::topology::HEX_8, 1, partVec_[0]); + + helperObjs.nodeAlg->add_kernel(*meta_); + + helperObjs.execute(); + + Kokkos::deep_copy(helperObjs.linsys->hostNumSumIntoCalls_, helperObjs.linsys->numSumIntoCalls_); + EXPECT_EQ(helperObjs.linsys->lhs_.extent(0), 8u); + EXPECT_EQ(helperObjs.linsys->lhs_.extent(1), 8u); + EXPECT_EQ(helperObjs.linsys->rhs_.extent(0), 8u); + EXPECT_EQ(helperObjs.linsys->hostNumSumIntoCalls_(0), 8u); + + namespace hex8_golds = hex8_golds::tke_ko; + unit_test_kernel_utils::expect_all_near( + helperObjs.linsys->rhs_, hex8_golds::rhs, 1.0e-12); + unit_test_kernel_utils::expect_all_near<8>( + helperObjs.linsys->lhs_, hex8_golds::lhs, 1.0e-12); +} + +TEST_F(KOKernelHex8Mesh, NGP_sdr_ko_node) +{ + // Only execute for 1 processor runs + if (bulk_->parallel_size() > 1) return; + + fill_mesh_and_init_fields(); + + // Setup solution options + solnOpts_.meshMotion_ = false; + solnOpts_.externalMeshDeformation_ = false; + solnOpts_.initialize_turbulence_constants(); + + unit_test_utils::NodeHelperObjects helperObjs( + bulk_, stk::topology::HEX_8, 1, partVec_[0]); + + helperObjs.nodeAlg->add_kernel(*meta_); + + helperObjs.execute(); + + Kokkos::deep_copy(helperObjs.linsys->hostNumSumIntoCalls_, helperObjs.linsys->numSumIntoCalls_); + EXPECT_EQ(helperObjs.linsys->lhs_.extent(0), 8u); + EXPECT_EQ(helperObjs.linsys->lhs_.extent(1), 8u); + EXPECT_EQ(helperObjs.linsys->rhs_.extent(0), 8u); + EXPECT_EQ(helperObjs.linsys->hostNumSumIntoCalls_(0), 8u); + + namespace hex8_golds = hex8_golds::sdr_ko; + unit_test_kernel_utils::expect_all_near( + helperObjs.linsys->rhs_, hex8_golds::rhs, 1.0e-12); + unit_test_kernel_utils::expect_all_near<8>( + helperObjs.linsys->lhs_, hex8_golds::lhs, 1.0e-12); +}