From 299e469f294e5fa02f0fc8d8842985af8a4c0eeb Mon Sep 17 00:00:00 2001 From: Marc Henry de Frahan Date: Mon, 2 Dec 2024 16:47:10 -0700 Subject: [PATCH] big update for ref temp --- .../icns/source_terms/ABLMeanBoussinesq.H | 6 +- .../icns/source_terms/ABLMeanBoussinesq.cpp | 73 +++++++++++++------ .../icns/source_terms/BoussinesqBuoyancy.H | 3 + .../icns/source_terms/BoussinesqBuoyancy.cpp | 4 +- .../source_terms/DragTempForcing.H | 8 +- .../source_terms/DragTempForcing.cpp | 9 ++- .../tke/source_terms/KransAxell.H | 8 +- .../tke/source_terms/KransAxell.cpp | 19 +++-- amr-wind/transport_models/ConstTransport.H | 14 ++++ amr-wind/transport_models/TransportModel.H | 3 + amr-wind/transport_models/TwoPhaseTransport.H | 8 ++ amr-wind/turbulence/LES/AMD.H | 3 - amr-wind/turbulence/LES/OneEqKsgs.H | 3 - amr-wind/turbulence/LES/OneEqKsgs.cpp | 2 +- amr-wind/turbulence/RANS/KLAxell.H | 2 - amr-wind/turbulence/RANS/KLAxell.cpp | 7 +- amr-wind/wind_energy/ABLStats.H | 4 +- amr-wind/wind_energy/ABLStats.cpp | 7 +- amr-wind/wind_energy/ABLWallFunction.cpp | 9 +-- amr-wind/wind_energy/MOData.H | 1 - unit_tests/turbulence/test_turbulence_LES.cpp | 4 +- 21 files changed, 131 insertions(+), 66 deletions(-) diff --git a/amr-wind/equation_systems/icns/source_terms/ABLMeanBoussinesq.H b/amr-wind/equation_systems/icns/source_terms/ABLMeanBoussinesq.H index 258299dd92..cbc3bc4070 100644 --- a/amr-wind/equation_systems/icns/source_terms/ABLMeanBoussinesq.H +++ b/amr-wind/equation_systems/icns/source_terms/ABLMeanBoussinesq.H @@ -43,15 +43,15 @@ private: amrex::Gpu::DeviceVector m_theta_ht; amrex::Gpu::DeviceVector m_theta_vals; - //! Reference temperature (Kelvin) - amrex::Real m_ref_theta{300.0}; - //! Transport model const transport::TransportModel& m_transport; //! Thermal expansion coefficient std::unique_ptr m_beta; + //! Reference temperature + std::unique_ptr m_ref_theta; + int m_axis{2}; bool m_const_profile{false}; diff --git a/amr-wind/equation_systems/icns/source_terms/ABLMeanBoussinesq.cpp b/amr-wind/equation_systems/icns/source_terms/ABLMeanBoussinesq.cpp index 6859b27ebf..aa897066ed 100644 --- a/amr-wind/equation_systems/icns/source_terms/ABLMeanBoussinesq.cpp +++ b/amr-wind/equation_systems/icns/source_terms/ABLMeanBoussinesq.cpp @@ -8,14 +8,6 @@ namespace amr_wind::pde::icns { /** Boussinesq buoyancy source term for ABL simulations - * - * Reads in the following parameters from `ABLMeanBoussinesq` namespace: - * - * - `reference_temperature` (Mandatory) temperature (`T0`) in Kelvin - * - `thermal_expansion_coeff` Optional, default = `1.0 / T0` - * - `gravity` acceleration due to gravity (m/s) - * - `read_temperature_profile` - * - `tprofile_filename` */ ABLMeanBoussinesq::ABLMeanBoussinesq(const CFDSim& sim) : m_mesh(sim.mesh()), m_transport(sim.transport_model()) @@ -24,31 +16,67 @@ ABLMeanBoussinesq::ABLMeanBoussinesq(const CFDSim& sim) const auto& abl = sim.physics_manager().get(); abl.register_mean_boussinesq_term(this); - amrex::ParmParse pp_boussinesq_buoyancy("BoussinesqBuoyancy"); - pp_boussinesq_buoyancy.get("reference_temperature", m_ref_theta); - m_beta = m_transport.beta(); + m_ref_theta = m_transport.ref_theta(); // gravity in `incflo` namespace amrex::ParmParse pp_incflo("incflo"); pp_incflo.queryarr("gravity", m_gravity); - bool read_temp_prof = false; - pp_boussinesq_buoyancy.query("read_temperature_profile", read_temp_prof); - - if ((!pp_boussinesq_buoyancy.contains("read_temperature_profile") && - pp_boussinesq_buoyancy.contains("tprofile_filename")) || - read_temp_prof) { + // Backwards compatibility + amrex::ParmParse pp_boussinesq_buoyancy("BoussinesqBuoyancy"); + amrex::ParmParse pp_abl("ABLMeanBoussinesq"); - m_const_profile = true; + bool read_temp_prof = false; + if (pp_abl.contains("read_temperature_profile")) { + pp_abl.get("read_temperature_profile", read_temp_prof); + if (pp_boussinesq_buoyancy.contains("read_temperature_profile")) { + amrex::Print() + << "WARNING: BoussinesqBuoyancy.read_temperature_profile " + "option has been deprecated in favor of " + "ABLMeanBoussinesq.read_temperature_profile. Ignoring the" + "BoussinesqBuoyancy option in favor of the " + "ABLMeanBoussinesq " + "option." + << std::endl; + } + } else if (pp_boussinesq_buoyancy.contains("read_temperature_profile")) { + amrex::Print() + << "WARNING: BoussinesqBuoyancy.read_temperature_profile option " + "has been deprecated in favor of " + "ABLMeanBoussinesq.read_temperature_profile. Please replace " + "this option." + << std::endl; + pp_boussinesq_buoyancy.get("read_temperature_profile", read_temp_prof); + } - std::string tprofile_filename; + std::string tprofile_filename; + if (pp_abl.contains("temperature_profile_filename")) { + pp_abl.get("temperature_profile_filename", tprofile_filename); + if (pp_boussinesq_buoyancy.contains("tprofile_filename")) { + amrex::Print() << "WARNING: BoussinesqBuoyancy.tprofile_filename " + "option has been deprecated in favor of " + "ABLMeanBoussinesq.temperature_profile_filename. " + "Ignoring the" + "BoussinesqBuoyancy option in favor of the " + "ABLMeanBoussinesq " + "option." + << std::endl; + } + } else if (pp_boussinesq_buoyancy.contains("tprofile_filename")) { + amrex::Print() + << "WARNING: BoussinesqBuoyancy.tprofile_filename option " + "has been deprecated in favor of " + "ABLMeanBoussinesq.temperature_profile_filename. Please replace " + "this option." + << std::endl; pp_boussinesq_buoyancy.get("tprofile_filename", tprofile_filename); + } + if ((read_temp_prof) && (tprofile_filename.empty())) { + m_const_profile = true; read_temperature_profile(tprofile_filename); - } else { - mean_temperature_init(abl.abl_statistics().theta_profile()); } } @@ -64,8 +92,8 @@ void ABLMeanBoussinesq::operator()( { const auto& problo = m_mesh.Geom(lev).ProbLoArray(); const auto& dx = m_mesh.Geom(lev).CellSizeArray(); - const amrex::Real T0 = m_ref_theta; const auto& beta = (*m_beta)(lev).const_array(mfi); + const auto& ref_theta = (*m_ref_theta)(lev).const_array(mfi); const amrex::GpuArray gravity{ m_gravity[0], m_gravity[1], m_gravity[2]}; @@ -79,6 +107,7 @@ void ABLMeanBoussinesq::operator()( amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { amrex::IntVect iv(i, j, k); const amrex::Real ht = problo[idir] + (iv[idir] + 0.5) * dx[idir]; + const amrex::Real T0 = ref_theta(i, j, k); const amrex::Real temp = amr_wind::interp::linear(theights, theights_end, tvals, ht); const amrex::Real fac = beta(i, j, k) * (temp - T0); diff --git a/amr-wind/equation_systems/icns/source_terms/BoussinesqBuoyancy.H b/amr-wind/equation_systems/icns/source_terms/BoussinesqBuoyancy.H index 205a24f478..13651ae0ea 100644 --- a/amr-wind/equation_systems/icns/source_terms/BoussinesqBuoyancy.H +++ b/amr-wind/equation_systems/icns/source_terms/BoussinesqBuoyancy.H @@ -40,6 +40,9 @@ private: //! Thermal expansion coefficient std::unique_ptr m_beta; + + //! Reference temperature + std::unique_ptr m_ref_theta; }; } // namespace amr_wind::pde::icns diff --git a/amr-wind/equation_systems/icns/source_terms/BoussinesqBuoyancy.cpp b/amr-wind/equation_systems/icns/source_terms/BoussinesqBuoyancy.cpp index a7169d17f6..61389b609f 100644 --- a/amr-wind/equation_systems/icns/source_terms/BoussinesqBuoyancy.cpp +++ b/amr-wind/equation_systems/icns/source_terms/BoussinesqBuoyancy.cpp @@ -14,6 +14,7 @@ BoussinesqBuoyancy::BoussinesqBuoyancy(const CFDSim& sim) , m_transport(sim.transport_model()) { m_beta = m_transport.beta(); + m_ref_theta = m_transport.ref_theta(); // gravity in `incflo` namespace amrex::ParmParse pp_incflo("incflo"); @@ -29,7 +30,6 @@ void BoussinesqBuoyancy::operator()( const FieldState fstate, const amrex::Array4& src_term) const { - const amrex::Real T0 = m_transport.reference_temperature(); const amrex::GpuArray gravity{ m_gravity[0], m_gravity[1], m_gravity[2]}; @@ -37,8 +37,10 @@ void BoussinesqBuoyancy::operator()( m_temperature.state(field_impl::phi_state(fstate))(lev).const_array( mfi); const auto& beta = (*m_beta)(lev).const_array(mfi); + const auto& ref_theta = (*m_ref_theta)(lev).const_array(mfi); amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { const amrex::Real T = temp(i, j, k, 0); + const amrex::Real T0 = ref_theta(i, j, k); const amrex::Real fac = beta(i, j, k) * (T0 - T); src_term(i, j, k, 0) += gravity[0] * fac; diff --git a/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.H b/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.H index 8c504d0ffb..6c2c27aafd 100644 --- a/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.H +++ b/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.H @@ -4,6 +4,7 @@ #include "amr-wind/equation_systems/temperature/TemperatureSource.H" #include "amr-wind/core/SimTime.H" #include "amr-wind/CFDSim.H" +#include "amr-wind/transport_models/TransportModel.H" namespace amr_wind::pde::temperature { @@ -29,7 +30,12 @@ private: const Field& m_velocity; const Field& m_temperature; amrex::Real m_drag_coefficient{1.0}; - amrex::Real m_reference_temperature{300.0}; + + //! Transport model + const transport::TransportModel& m_transport; + + //! Reference temperature + std::unique_ptr m_ref_theta; }; } // namespace amr_wind::pde::temperature diff --git a/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.cpp b/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.cpp index 098e806cc0..ae1b8e38ee 100644 --- a/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.cpp +++ b/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.cpp @@ -13,10 +13,11 @@ DragTempForcing::DragTempForcing(const CFDSim& sim) , m_mesh(sim.mesh()) , m_velocity(sim.repo().get_field("velocity")) , m_temperature(sim.repo().get_field("temperature")) + , m_transport(sim.transport_model()) { amrex::ParmParse pp("DragTempForcing"); pp.query("drag_coefficient", m_drag_coefficient); - pp.query("reference_temperature", m_reference_temperature); + m_ref_theta = m_transport.ref_theta(); } DragTempForcing::~DragTempForcing() = default; @@ -44,7 +45,7 @@ void DragTempForcing::operator()( const auto& geom = m_mesh.Geom(lev); const auto& dx = geom.CellSizeArray(); const amrex::Real drag_coefficient = m_drag_coefficient / dx[2]; - const amrex::Real reference_temperature = m_reference_temperature; + const auto& ref_theta = (*m_ref_theta)(lev).const_array(mfi); const auto tiny = std::numeric_limits::epsilon(); const amrex::Real cd_max = 10.0; amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { @@ -54,9 +55,9 @@ void DragTempForcing::operator()( const amrex::Real m = std::sqrt(ux1 * ux1 + uy1 * uy1 + uz1 * uz1); const amrex::Real Cd = std::min(drag_coefficient / (m + tiny), cd_max / dx[2]); + const amrex::Real T0 = ref_theta(i, j, k); src_term(i, j, k, 0) -= - (Cd * (temperature(i, j, k, 0) - reference_temperature) * - blank(i, j, k, 0)); + (Cd * (temperature(i, j, k, 0) - T0) * blank(i, j, k, 0)); }); } diff --git a/amr-wind/equation_systems/tke/source_terms/KransAxell.H b/amr-wind/equation_systems/tke/source_terms/KransAxell.H index dc4bbe2066..a8e0f0d80a 100644 --- a/amr-wind/equation_systems/tke/source_terms/KransAxell.H +++ b/amr-wind/equation_systems/tke/source_terms/KransAxell.H @@ -1,6 +1,7 @@ #ifndef KRANSAXELL_H #define KRANSAXELL_H +#include "amr-wind/transport_models/TransportModel.H" #include "amr-wind/equation_systems/tke/TKESource.H" namespace amr_wind::pde::tke { @@ -35,7 +36,6 @@ private: Field& m_tke; amrex::Real m_Cmu{0.556}; amrex::Real m_heat_flux{0.0}; - amrex::Real m_ref_temp{300.0}; amrex::Real m_z0{0.1}; amrex::Real m_kappa{0.41}; amrex::Real m_sponge_start{600}; @@ -45,6 +45,12 @@ private: const CFDSim& m_sim; const amrex::AmrCore& m_mesh; const Field& m_velocity; + + //! Transport model + const transport::TransportModel& m_transport; + + //! Reference temperature + std::unique_ptr m_ref_theta; }; } // namespace amr_wind::pde::tke diff --git a/amr-wind/equation_systems/tke/source_terms/KransAxell.cpp b/amr-wind/equation_systems/tke/source_terms/KransAxell.cpp index 66808678af..df24297bbe 100644 --- a/amr-wind/equation_systems/tke/source_terms/KransAxell.cpp +++ b/amr-wind/equation_systems/tke/source_terms/KransAxell.cpp @@ -16,6 +16,7 @@ KransAxell::KransAxell(const CFDSim& sim) , m_sim(sim) , m_mesh(sim.mesh()) , m_velocity(sim.repo().get_field("velocity")) + , m_transport(sim.transport_model()) { AMREX_ALWAYS_ASSERT(sim.turbulence_model().model_name() == "KLAxell"); auto coeffs = sim.turbulence_model().model_coeffs(); @@ -23,13 +24,13 @@ KransAxell::KransAxell(const CFDSim& sim) pp.query("Cmu", m_Cmu); pp.query("kappa", m_kappa); pp.query("surface_roughness_z0", m_z0); - pp.query("reference_temperature", m_ref_temp); pp.query("surface_temp_flux", m_heat_flux); pp.query("meso_sponge_start", m_sponge_start); { amrex::ParmParse pp_incflow("incflo"); pp_incflow.queryarr("gravity", m_gravity); } + m_ref_theta = m_transport.ref_theta(); } KransAxell::~KransAxell() = default; @@ -48,20 +49,21 @@ void KransAxell::operator()( const auto& buoy_prod_arr = (this->m_buoy_prod)(lev).array(mfi); const auto& dissip_arr = (this->m_dissip)(lev).array(mfi); const auto& tke_arr = m_tke(lev).array(mfi); + const auto& ref_theta_arr = (*m_ref_theta)(lev).const_array(mfi); const auto& geom = m_mesh.Geom(lev); const auto& problo = m_mesh.Geom(lev).ProbLoArray(); const auto& probhi = m_mesh.Geom(lev).ProbHiArray(); const auto& dx = geom.CellSizeArray(); const auto& dt = m_time.delta_t(); - const amrex::Real ref_temp = m_ref_temp; - const amrex::Real heat_flux = - std::abs(m_gravity[2]) / ref_temp * m_heat_flux; + const amrex::Real heat_flux = m_heat_flux; const amrex::Real Cmu = m_Cmu; const amrex::Real sponge_start = m_sponge_start; const amrex::Real ref_tke = m_ref_tke; const auto tiny = std::numeric_limits::epsilon(); const amrex::Real kappa = m_kappa; const amrex::Real z0 = m_z0; + const amrex::GpuArray gravity{ + m_gravity[0], m_gravity[1], m_gravity[2]}; amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { amrex::Real bcforcing = 0; const amrex::Real ux = vel(i, j, k, 0); @@ -70,9 +72,10 @@ void KransAxell::operator()( if (k == 0) { const amrex::Real m = std::sqrt(ux * ux + uy * uy); const amrex::Real ustar = m * kappa / std::log(z / z0); + const amrex::Real T0 = ref_theta_arr(i, j, k); + const amrex::Real hf = std::abs(gravity[2]) / T0 * heat_flux; const amrex::Real rans_b = std::pow( - std::max(heat_flux, 0.0) * kappa * z / std::pow(Cmu, 3), - (2.0 / 3.0)); + std::max(hf, 0.0) * kappa * z / std::pow(Cmu, 3), (2.0 / 3.0)); bcforcing = (ustar * ustar / (Cmu * Cmu) + rans_b - tke_arr(i, j, k)) / dt; } @@ -105,8 +108,10 @@ void KransAxell::operator()( const amrex::Real z = 0.5 * dx[2]; amrex::Real m = std::sqrt(ux * ux + uy * uy); const amrex::Real ustar = m * kappa / std::log(z / z0); + const amrex::Real T0 = ref_theta_arr(i, j, k); + const amrex::Real hf = std::abs(gravity[2]) / T0 * heat_flux; const amrex::Real rans_b = std::pow( - std::max(heat_flux, 0.0) * kappa * z / std::pow(Cmu, 3), + std::max(hf, 0.0) * kappa * z / std::pow(Cmu, 3), (2.0 / 3.0)); terrainforcing = (ustar * ustar / (Cmu * Cmu) + rans_b - tke_arr(i, j, k)) / diff --git a/amr-wind/transport_models/ConstTransport.H b/amr-wind/transport_models/ConstTransport.H index fbe625bb97..5ac3cdc096 100644 --- a/amr-wind/transport_models/ConstTransport.H +++ b/amr-wind/transport_models/ConstTransport.H @@ -198,6 +198,20 @@ public: return m_reference_temperature; } + //! Return the reference temperature + inline std::unique_ptr ref_theta() const override + { + if (m_reference_temperature < 0.0) { + amrex::Abort("reference temperature was not set"); + } + + auto ref_theta = m_repo.create_scratch_field(1, 1); + for (int lev = 0; lev < m_repo.num_active_levels(); ++lev) { + (*ref_theta)(lev).setVal(m_reference_temperature); + } + return ref_theta; + } + private: //! Reference to the field repository (for creating scratch fields) FieldRepo& m_repo; diff --git a/amr-wind/transport_models/TransportModel.H b/amr-wind/transport_models/TransportModel.H index 60e553fc4f..ef1875939c 100644 --- a/amr-wind/transport_models/TransportModel.H +++ b/amr-wind/transport_models/TransportModel.H @@ -48,6 +48,9 @@ public: //! Reference temperature virtual amrex::Real reference_temperature() const = 0; + + //! Reference temperature + virtual std::unique_ptr ref_theta() const = 0; }; } // namespace amr_wind::transport diff --git a/amr-wind/transport_models/TwoPhaseTransport.H b/amr-wind/transport_models/TwoPhaseTransport.H index 4fada7deae..538cd5a409 100644 --- a/amr-wind/transport_models/TwoPhaseTransport.H +++ b/amr-wind/transport_models/TwoPhaseTransport.H @@ -356,6 +356,14 @@ public: return m_reference_temperature; } + //! Return the reference temperature + inline std::unique_ptr ref_theta() const override + { + amrex::Abort("ref temp not implemented"); + auto ref_theta = m_repo.create_scratch_field(1, 1); + return ref_theta; + } + private: //! Reference to the CFD sim const CFDSim& m_sim; diff --git a/amr-wind/turbulence/LES/AMD.H b/amr-wind/turbulence/LES/AMD.H index 5ca475a7ee..df7e1f6bc3 100644 --- a/amr-wind/turbulence/LES/AMD.H +++ b/amr-wind/turbulence/LES/AMD.H @@ -50,9 +50,6 @@ private: //! discretization) amrex::Real m_C{0.333333333333333}; - //! Reference temperature (Kelvin) - amrex::Real m_ref_theta{300.0}; - //! Wall-normal direction axis int m_normal_dir{2}; diff --git a/amr-wind/turbulence/LES/OneEqKsgs.H b/amr-wind/turbulence/LES/OneEqKsgs.H index 4467470c75..5ce88c053b 100644 --- a/amr-wind/turbulence/LES/OneEqKsgs.H +++ b/amr-wind/turbulence/LES/OneEqKsgs.H @@ -80,9 +80,6 @@ private: //! Gravity vector (m/s^2) amrex::Vector m_gravity{0.0, 0.0, -9.81}; - //! Reference temperature (Kelvin) - amrex::Real m_ref_theta{300.0}; - //! Hybrid RANS-LES with Nalu-wind bool m_hybrid_rl{false}; Field* m_sdr; diff --git a/amr-wind/turbulence/LES/OneEqKsgs.cpp b/amr-wind/turbulence/LES/OneEqKsgs.cpp index ec4cb11e18..d7fafbcb60 100644 --- a/amr-wind/turbulence/LES/OneEqKsgs.cpp +++ b/amr-wind/turbulence/LES/OneEqKsgs.cpp @@ -126,7 +126,7 @@ void OneEqKsgsM84::update_turbulent_viscosity( const auto& tke_arr = (*this->m_tke)(lev).array(mfi); const auto& buoy_prod_arr = (this->m_buoy_prod)(lev).array(mfi); const auto& shear_prod_arr = (this->m_shear_prod)(lev).array(mfi); - const auto& beta_arr = (*beta)(lev).array(mfi); + const auto& beta_arr = (*beta)(lev).const_array(mfi); amrex::ParallelFor( bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { diff --git a/amr-wind/turbulence/RANS/KLAxell.H b/amr-wind/turbulence/RANS/KLAxell.H index 1f6175f0b3..af19ae1246 100644 --- a/amr-wind/turbulence/RANS/KLAxell.H +++ b/amr-wind/turbulence/RANS/KLAxell.H @@ -58,8 +58,6 @@ private: Field& m_temperature; //! Gravity vector (m/s^2) amrex::Vector m_gravity{0.0, 0.0, -9.81}; - //! Reference temperature (Kelvin) - amrex::Real m_ref_theta{300.0}; amrex::Real m_surf_flux{0}; amrex::Real m_lengthscale_switch{800}; }; diff --git a/amr-wind/turbulence/RANS/KLAxell.cpp b/amr-wind/turbulence/RANS/KLAxell.cpp index b48bf7bd48..1d2af84972 100644 --- a/amr-wind/turbulence/RANS/KLAxell.cpp +++ b/amr-wind/turbulence/RANS/KLAxell.cpp @@ -112,7 +112,7 @@ void KLAxell::update_turbulent_viscosity( const auto& tke_arr = (*this->m_tke)(lev).array(mfi); const auto& buoy_prod_arr = (this->m_buoy_prod)(lev).array(mfi); const auto& shear_prod_arr = (this->m_shear_prod)(lev).array(mfi); - const auto& beta_arr = (*beta)(lev).array(mfi); + const auto& beta_arr = (*beta)(lev).const_array(mfi); //! Add terrain components const bool has_terrain = @@ -281,7 +281,7 @@ void KLAxell::update_alphaeff(Field& alphaeff) fvm::gradient(*gradT, m_temperature); const amrex::GpuArray gravity{ m_gravity[0], m_gravity[1], m_gravity[2]}; - const amrex::Real beta = 1.0 / m_ref_theta; + const auto beta = (this->m_transport).beta(); const amrex::Real Cmu = m_Cmu; const int nlevels = repo.num_active_levels(); for (int lev = 0; lev < nlevels; ++lev) { @@ -293,6 +293,7 @@ void KLAxell::update_alphaeff(Field& alphaeff) const auto& tke_arr = (*this->m_tke)(lev).array(mfi); const auto& gradT_arr = (*gradT)(lev).array(mfi); const auto& tlscale_arr = (this->m_turb_lscale)(lev).array(mfi); + const auto& beta_arr = (*beta)(lev).const_array(mfi); const amrex::Real Rtc = -1.0; const amrex::Real Rtmin = -3.0; amrex::ParallelFor( @@ -301,7 +302,7 @@ void KLAxell::update_alphaeff(Field& alphaeff) -(gradT_arr(i, j, k, 0) * gravity[0] + gradT_arr(i, j, k, 1) * gravity[1] + gradT_arr(i, j, k, 2) * gravity[2]) * - beta; + beta_arr(i, j, k); amrex::Real epsilon = std::pow(Cmu, 3) * std::pow(tke_arr(i, j, k), 1.5) / tlscale_arr(i, j, k); diff --git a/amr-wind/wind_energy/ABLStats.H b/amr-wind/wind_energy/ABLStats.H index 633d811edb..1ba8269e3e 100644 --- a/amr-wind/wind_energy/ABLStats.H +++ b/amr-wind/wind_energy/ABLStats.H @@ -3,6 +3,7 @@ #include "amr-wind/wind_energy/ABLStatsBase.H" #include "amr-wind/CFDSim.H" +#include "amr-wind/transport_models/TransportModel.H" #include "amr-wind/utilities/FieldPlaneAveraging.H" #include "amr-wind/utilities/FieldPlaneAveragingFine.H" #include "amr-wind/utilities/SecondMomentAveraging.H" @@ -176,9 +177,6 @@ private: //! Von-Karman constant amrex::Real m_kappa{0.41}; - //! Reference surface temperature - amrex::Real m_ref_theta{300.0}; - //! Variable to store capping inversion height amrex::Real m_zi{0.0}; diff --git a/amr-wind/wind_energy/ABLStats.cpp b/amr-wind/wind_energy/ABLStats.cpp index 95517285cc..ad0aadea50 100644 --- a/amr-wind/wind_energy/ABLStats.cpp +++ b/amr-wind/wind_energy/ABLStats.cpp @@ -55,7 +55,6 @@ void ABLStats::initialize() pp.query("normal_direction", m_normal_dir); AMREX_ASSERT((0 <= m_normal_dir) && (m_normal_dir < AMREX_SPACEDIM)); pp.query("kappa", m_kappa); - pp.get("reference_temperature", m_ref_theta); pp.query("stats_do_energy_budget", m_do_energy_budget); } @@ -337,9 +336,10 @@ void ABLStats::write_ascii() } double wstar = 0.0; + const auto ref_theta = m_sim.transport_model().reference_temperature(); auto Q = m_abl_wall_func.mo().surf_temp_flux; if (Q > 1e-10) { - wstar = std::cbrt(m_gravity * Q * m_zi / m_ref_theta); + wstar = std::cbrt(m_gravity * Q * m_zi / ref_theta); } auto L = m_abl_wall_func.mo().obukhov_len; @@ -525,11 +525,12 @@ void ABLStats::write_netcdf() ncf.var("ustar").put(&ustar, {nt}, {1}); double wstar = 0.0; auto Q = m_abl_wall_func.mo().surf_temp_flux; + const auto ref_theta = m_sim.transport_model().reference_temperature(); ncf.var("Q").put(&Q, {nt}, {1}); auto Tsurf = m_abl_wall_func.mo().surf_temp; ncf.var("Tsurf").put(&Tsurf, {nt}, {1}); if (Q > 1e-10) { - wstar = std::cbrt(m_gravity * Q * m_zi / m_ref_theta); + wstar = std::cbrt(m_gravity * Q * m_zi / ref_theta); } ncf.var("wstar").put(&wstar, {nt}, {1}); double L = m_abl_wall_func.mo().obukhov_len; diff --git a/amr-wind/wind_energy/ABLWallFunction.cpp b/amr-wind/wind_energy/ABLWallFunction.cpp index d4e101eea0..6be8af6074 100644 --- a/amr-wind/wind_energy/ABLWallFunction.cpp +++ b/amr-wind/wind_energy/ABLWallFunction.cpp @@ -1,3 +1,4 @@ +#include "amr-wind/transport_models/TransportModel.H" #include "amr-wind/wind_energy/ABLWallFunction.H" #include "amr-wind/wind_energy/ABL.H" #include "amr-wind/utilities/tensor_ops.H" @@ -42,8 +43,6 @@ ABLWallFunction::ABLWallFunction(const CFDSim& sim) << std::endl; } - pp.get("reference_temperature", m_mo.ref_temp); - if (pp.contains("surface_temp_flux")) { pp.query("surface_temp_flux", m_mo.surf_temp_flux); } else if (pp.contains("surface_temp_rate")) { @@ -54,9 +53,9 @@ ABLWallFunction::ABLWallFunction(const CFDSim& sim) } else { amrex::Print() << "ABLWallFunction: Initial surface temperature not found for " - "ABL. Assuming to be equal to the reference temperature " - << m_mo.ref_temp << std::endl; - m_surf_temp_init = m_mo.ref_temp; + "ABL. Assuming to be equal to the reference temperature" + << std::endl; + m_surf_temp_init = sim.transport_model().reference_temperature(); } if (pp.contains("surface_temp_rate_tstart")) { pp.get("surface_temp_rate_tstart", m_surf_temp_rate_tstart); diff --git a/amr-wind/wind_energy/MOData.H b/amr-wind/wind_energy/MOData.H index a19024f352..d37dab2fbd 100644 --- a/amr-wind/wind_energy/MOData.H +++ b/amr-wind/wind_energy/MOData.H @@ -39,7 +39,6 @@ struct MOData amrex::Real surf_temp_flux{0.0}; ///< Heat flux amrex::Real surf_temp; ///< Instantaneous surface temperature - amrex::Real ref_temp; ///< Reference temperature amrex::Real gamma_m{5.0}; amrex::Real gamma_h{5.0}; diff --git a/unit_tests/turbulence/test_turbulence_LES.cpp b/unit_tests/turbulence/test_turbulence_LES.cpp index 088b9c97fa..cb767e7e26 100644 --- a/unit_tests/turbulence/test_turbulence_LES.cpp +++ b/unit_tests/turbulence/test_turbulence_LES.cpp @@ -276,11 +276,9 @@ TEST_F(TurbLESTest, test_1eqKsgs_setup_calc) initialize_mesh(); auto& pde_mgr = sim().pde_manager(); pde_mgr.register_icns(); + sim().create_transport_model(); sim().init_physics(); - - // Create turbulence model sim().create_turbulence_model(); - // Get turbulence model auto& tmodel = sim().turbulence_model(); // Get coefficients