diff --git a/amr-wind/equation_systems/temperature/source_terms/CMakeLists.txt b/amr-wind/equation_systems/temperature/source_terms/CMakeLists.txt index 9fa2e2057f..e271ef35e6 100644 --- a/amr-wind/equation_systems/temperature/source_terms/CMakeLists.txt +++ b/amr-wind/equation_systems/temperature/source_terms/CMakeLists.txt @@ -2,4 +2,5 @@ target_sources(${amr_wind_lib_name} PRIVATE ABLMesoForcingTemp.cpp BodyForce.cpp HurricaneTempForcing.cpp DragTempForcing.cpp +TempSpongeForcing.cpp ) diff --git a/amr-wind/equation_systems/temperature/source_terms/TempSpongeForcing.H b/amr-wind/equation_systems/temperature/source_terms/TempSpongeForcing.H new file mode 100644 index 0000000000..494d6b7bac --- /dev/null +++ b/amr-wind/equation_systems/temperature/source_terms/TempSpongeForcing.H @@ -0,0 +1,37 @@ +#ifndef TEMPSPONGEFORCING_H +#define TEMPSPONGEFORCING_H + +#include "amr-wind/equation_systems/temperature/TemperatureSource.H" +#include "amr-wind/core/SimTime.H" +#include "amr-wind/CFDSim.H" + +namespace amr_wind::pde::temperature { + +class TempSpongeForcing : public TemperatureSource::Register +{ +public: + static std::string identifier() { return "TempSpongeForcing"; } + + explicit TempSpongeForcing(const CFDSim& sim); + + ~TempSpongeForcing() override; + + void operator()( + const int lev, + const amrex::MFIter& mfi, + const amrex::Box& bx, + const FieldState /*fstate*/, + const amrex::Array4& src_term) const override; + +private: + const amrex::AmrCore& m_mesh; + const Field& m_temperature; + amrex::Vector m_theta_heights; + amrex::Vector m_theta_values; + amrex::Gpu::DeviceVector m_theta_heights_d; + amrex::Gpu::DeviceVector m_theta_values_d; + amrex::Real m_sponge_start{600}; +}; + +} // namespace amr_wind::pde::temperature +#endif \ No newline at end of file diff --git a/amr-wind/equation_systems/temperature/source_terms/TempSpongeForcing.cpp b/amr-wind/equation_systems/temperature/source_terms/TempSpongeForcing.cpp new file mode 100644 index 0000000000..525b6dfef9 --- /dev/null +++ b/amr-wind/equation_systems/temperature/source_terms/TempSpongeForcing.cpp @@ -0,0 +1,68 @@ + +#include "amr-wind/equation_systems/temperature/source_terms/TempSpongeForcing.H" +#include "amr-wind/utilities/IOManager.H" +#include "amr-wind/utilities/linear_interpolation.H" + +#include "AMReX_ParmParse.H" +#include "AMReX_Gpu.H" +#include "AMReX_Random.H" + +namespace amr_wind::pde::temperature { + +TempSpongeForcing::TempSpongeForcing(const CFDSim& sim) + : m_mesh(sim.mesh()), m_temperature(sim.repo().get_field("temperature")) +{ + amrex::ParmParse pp_abl("ABL"); + //! Temperature variation as a function of height + pp_abl.query("meso_sponge_start", m_sponge_start); + pp_abl.getarr("temperature_heights", m_theta_heights); + pp_abl.getarr("temperature_values", m_theta_values); + AMREX_ALWAYS_ASSERT(m_theta_heights.size() == m_theta_values.size()); + const int num_theta_values = static_cast(m_theta_heights.size()); + m_theta_heights_d.resize(num_theta_values); + m_theta_values_d.resize(num_theta_values); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, m_theta_heights.begin(), + m_theta_heights.end(), m_theta_heights_d.begin()); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, m_theta_values.begin(), m_theta_values.end(), + m_theta_values_d.begin()); +} + +TempSpongeForcing::~TempSpongeForcing() = default; + +void TempSpongeForcing::operator()( + const int lev, + const amrex::MFIter& mfi, + const amrex::Box& bx, + const FieldState fstate, + const amrex::Array4& src_term) const +{ + const auto& geom = m_mesh.Geom(lev); + const auto& dx = geom.CellSizeArray(); + const auto& prob_lo = geom.ProbLoArray(); + const auto& prob_hi = geom.ProbHiArray(); + const auto& temperature = + m_temperature.state(field_impl::dof_state(fstate))(lev).const_array( + mfi); + const amrex::Real sponge_start = m_sponge_start; + const auto vsize = m_theta_heights_d.size(); + const auto* theta_heights_d = m_theta_heights_d.data(); + const auto* theta_values_d = m_theta_values_d.data(); + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + const amrex::Real z = prob_lo[2] + (k + 0.5) * dx[2]; + const amrex::Real zi = + std::max((z - sponge_start) / (prob_hi[2] - sponge_start), 0.0); + amrex::Real ref_temp = temperature(i, j, k); + if (zi > 0) { + ref_temp = (vsize > 0) + ? interp::linear( + theta_heights_d, theta_heights_d + vsize, + theta_values_d, z) + : temperature(i, j, k); + } + src_term(i, j, k, 0) -= zi * zi * (temperature(i, j, k) - ref_temp); + }); +} + +} // namespace amr_wind::pde::temperature diff --git a/amr-wind/equation_systems/tke/source_terms/CMakeLists.txt b/amr-wind/equation_systems/tke/source_terms/CMakeLists.txt index beed97be50..c9a97bc713 100644 --- a/amr-wind/equation_systems/tke/source_terms/CMakeLists.txt +++ b/amr-wind/equation_systems/tke/source_terms/CMakeLists.txt @@ -1,4 +1,5 @@ target_sources(${amr_wind_lib_name} PRIVATE KsgsM84Src.cpp KwSSTSrc.cpp + KransAxell.cpp ) diff --git a/amr-wind/equation_systems/tke/source_terms/KransAxell.H b/amr-wind/equation_systems/tke/source_terms/KransAxell.H new file mode 100644 index 0000000000..dc4bbe2066 --- /dev/null +++ b/amr-wind/equation_systems/tke/source_terms/KransAxell.H @@ -0,0 +1,52 @@ +#ifndef KRANSAXELL_H +#define KRANSAXELL_H + +#include "amr-wind/equation_systems/tke/TKESource.H" + +namespace amr_wind::pde::tke { + +/** TKE source term based on Axell 2011 paper + * Axell, L. B., & Liungman, O. (2001). A one-equation turbulence model for + * geophysical applications: comparison with data and the k− ε model. + * Environmental Fluid Mechanics, 1, 71-106. + * \ingroup tke_src turb_model we_abl + */ +class KransAxell : public TKESource::Register +{ +public: + static std::string identifier() { return "KransAxell"; } + + explicit KransAxell(const CFDSim& /*sim*/); + + ~KransAxell() override; + + void operator()( + const int lev, + const amrex::MFIter& mfi, + const amrex::Box& bx, + const FieldState fstate, + const amrex::Array4& src_term) const override; + +private: + Field& m_turb_lscale; + Field& m_shear_prod; + Field& m_buoy_prod; + Field& m_dissip; + 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}; + amrex::Real m_ref_tke{1e-10}; + amrex::Vector m_gravity{0.0, 0.0, -9.81}; + const SimTime& m_time; + const CFDSim& m_sim; + const amrex::AmrCore& m_mesh; + const Field& m_velocity; +}; + +} // namespace amr_wind::pde::tke + +#endif /* KRANSAXELL_H */ diff --git a/amr-wind/equation_systems/tke/source_terms/KransAxell.cpp b/amr-wind/equation_systems/tke/source_terms/KransAxell.cpp new file mode 100644 index 0000000000..66808678af --- /dev/null +++ b/amr-wind/equation_systems/tke/source_terms/KransAxell.cpp @@ -0,0 +1,126 @@ +#include + +#include "amr-wind/equation_systems/tke/source_terms/KransAxell.H" +#include "amr-wind/CFDSim.H" +#include "amr-wind/turbulence/TurbulenceModel.H" +#include "amr-wind/utilities/linear_interpolation.H" +namespace amr_wind::pde::tke { + +KransAxell::KransAxell(const CFDSim& sim) + : m_turb_lscale(sim.repo().get_field("turb_lscale")) + , m_shear_prod(sim.repo().get_field("shear_prod")) + , m_buoy_prod(sim.repo().get_field("buoy_prod")) + , m_dissip(sim.repo().get_field("dissipation")) + , m_tke(sim.repo().get_field("tke")) + , m_time(sim.time()) + , m_sim(sim) + , m_mesh(sim.mesh()) + , m_velocity(sim.repo().get_field("velocity")) +{ + AMREX_ALWAYS_ASSERT(sim.turbulence_model().model_name() == "KLAxell"); + auto coeffs = sim.turbulence_model().model_coeffs(); + amrex::ParmParse pp("ABL"); + 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); + } +} + +KransAxell::~KransAxell() = default; + +void KransAxell::operator()( + const int lev, + const amrex::MFIter& mfi, + const amrex::Box& bx, + const FieldState fstate, + const amrex::Array4& src_term) const +{ + const auto& vel = + m_velocity.state(field_impl::dof_state(fstate))(lev).const_array(mfi); + const auto& tlscale_arr = (this->m_turb_lscale)(lev).array(mfi); + const auto& shear_prod_arr = (this->m_shear_prod)(lev).array(mfi); + 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& 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 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; + 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); + const amrex::Real uy = vel(i, j, k, 1); + const amrex::Real z = problo[2] + (k + 0.5) * dx[2]; + 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 rans_b = std::pow( + std::max(heat_flux, 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; + } + const amrex::Real zi = + std::max((z - sponge_start) / (probhi[2] - sponge_start), 0.0); + const amrex::Real sponge_forcing = + zi * zi * (tke_arr(i, j, k) - ref_tke); + dissip_arr(i, j, k) = std::pow(Cmu, 3) * + std::pow(tke_arr(i, j, k), 1.5) / + (tlscale_arr(i, j, k) + tiny); + src_term(i, j, k) += shear_prod_arr(i, j, k) + buoy_prod_arr(i, j, k) - + dissip_arr(i, j, k) - sponge_forcing + bcforcing; + }); + // Add terrain components + const bool has_terrain = + this->m_sim.repo().int_field_exists("terrain_blank"); + if (has_terrain) { + const auto* const m_terrain_blank = + &this->m_sim.repo().get_int_field("terrain_blank"); + const auto* const m_terrain_drag = + &this->m_sim.repo().get_int_field("terrain_drag"); + const auto& blank_arr = (*m_terrain_blank)(lev).const_array(mfi); + const auto& drag_arr = (*m_terrain_drag)(lev).const_array(mfi); + amrex::ParallelFor( + bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + amrex::Real terrainforcing = 0; + amrex::Real dragforcing = 0; + const amrex::Real ux = vel(i, j, k, 0); + const amrex::Real uy = vel(i, j, k, 1); + 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 rans_b = std::pow( + std::max(heat_flux, 0.0) * kappa * z / std::pow(Cmu, 3), + (2.0 / 3.0)); + terrainforcing = + (ustar * ustar / (Cmu * Cmu) + rans_b - tke_arr(i, j, k)) / + dt; + const amrex::Real uz = vel(i, j, k, 2); + m = std::sqrt(ux * ux + uy * uy + uz * uz); + const amrex::Real Cd = + std::min(10 / (dx[2] * m + tiny), 100 / dx[2]); + dragforcing = -Cd * m * tke_arr(i, j, k, 0); + + src_term(i, j, k) += drag_arr(i, j, k) * terrainforcing + + blank_arr(i, j, k) * dragforcing; + }); + } +} + +} // namespace amr_wind::pde::tke diff --git a/amr-wind/turbulence/RANS/CMakeLists.txt b/amr-wind/turbulence/RANS/CMakeLists.txt index 137db454c3..342f7d4fc9 100644 --- a/amr-wind/turbulence/RANS/CMakeLists.txt +++ b/amr-wind/turbulence/RANS/CMakeLists.txt @@ -1,4 +1,5 @@ target_sources(${amr_wind_lib_name} PRIVATE KOmegaSST.cpp KOmegaSSTIDDES.cpp + KLAxell.cpp ) diff --git a/amr-wind/turbulence/RANS/KLAxell.H b/amr-wind/turbulence/RANS/KLAxell.H new file mode 100644 index 0000000000..1f6175f0b3 --- /dev/null +++ b/amr-wind/turbulence/RANS/KLAxell.H @@ -0,0 +1,69 @@ +#ifndef KLAXELL_H +#define KLAXELL_H + +#include +#include "amr-wind/turbulence/TurbModelBase.H" + +namespace amr_wind::turbulence { + +/** Base class for 1-Equation RANS TKE turbulence model + * \ingroup turb_model + */ +template +class KLAxell : public TurbModelBase +{ +public: + static std::string identifier() + { + return "KLAxell-" + Transport::identifier(); + } + + explicit KLAxell(CFDSim& sim); + + std::string model_name() const override { return "KLAxell"; } + + //! Update the turbulent viscosity field + void update_turbulent_viscosity( + const FieldState fstate, const DiffusionType /*unused*/) override; + + //! Do any post advance work + void post_advance_work() override; + + //! Update the effective thermal diffusivity field + void update_alphaeff(Field& alphaeff) override; + + //! Update the effective scalar diffusivity field + void update_scalar_diff(Field& deff, const std::string& name) override; + + //! Parse turbulence model coefficients + void parse_model_coeffs() override; + + //! Return turbulence model coefficients + TurbulenceModel::CoeffsDictType model_coeffs() const override; + +private: + Field& m_vel; + Field& m_turb_lscale; + Field& m_shear_prod; + Field& m_buoy_prod; + Field& m_dissip; + Field& m_rho; + Field* m_tke{nullptr}; + //! Turbulence constant + amrex::Real m_Cmu{0.556}; + amrex::Real m_Cmu_prime{0.556}; + amrex::Real m_Cb_stable{0.25}; + amrex::Real m_Cb_unstable{0.35}; + amrex::Real m_prandtl{1.0}; + 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}; +}; + +} // namespace amr_wind::turbulence + +#endif /* KLAXELL_H */ diff --git a/amr-wind/turbulence/RANS/KLAxell.cpp b/amr-wind/turbulence/RANS/KLAxell.cpp new file mode 100644 index 0000000000..edf25a8ed2 --- /dev/null +++ b/amr-wind/turbulence/RANS/KLAxell.cpp @@ -0,0 +1,353 @@ +#include "amr-wind/turbulence/RANS/KLAxell.H" +#include "amr-wind/equation_systems/PDEBase.H" +#include "amr-wind/turbulence/TurbModelDefs.H" +#include "amr-wind/fvm/gradient.H" +#include "amr-wind/fvm/strainrate.H" +#include "amr-wind/turbulence/turb_utils.H" +#include "amr-wind/equation_systems/tke/TKE.H" + +#include "AMReX_ParmParse.H" + +namespace amr_wind { +namespace turbulence { + +template +KLAxell::KLAxell(CFDSim& sim) + : TurbModelBase(sim) + , m_vel(sim.repo().get_field("velocity")) + , m_turb_lscale(sim.repo().declare_field("turb_lscale", 1)) + , m_shear_prod(sim.repo().declare_field("shear_prod", 1)) + , m_buoy_prod(sim.repo().declare_field("buoy_prod", 1)) + , m_dissip(sim.repo().declare_field("dissipation", 1)) + , m_rho(sim.repo().get_field("density")) + , m_temperature(sim.repo().get_field("temperature")) +{ + auto& tke_eqn = + sim.pde_manager().register_transport_pde(pde::TKE::pde_name()); + m_tke = &(tke_eqn.fields().field); + auto& phy_mgr = this->m_sim.physics_manager(); + if (!phy_mgr.contains("ABL")) { + amrex::Abort("KLAxell model only works with ABL physics"); + } + { + amrex::ParmParse pp("ABL"); + pp.get("reference_temperature", m_ref_theta); + pp.get("surface_temp_flux", m_surf_flux); + pp.query("length_scale_switch", m_lengthscale_switch); + } + + { + amrex::ParmParse pp("incflo"); + pp.queryarr("gravity", m_gravity); + } + + // TKE source term to be added to PDE + turb_utils::inject_turbulence_src_terms( + pde::TKE::pde_name(), {"KransAxell"}); +} + +template +void KLAxell::parse_model_coeffs() +{ + const std::string coeffs_dict = this->model_name() + "_coeffs"; + amrex::ParmParse pp(coeffs_dict); + pp.query("Cmu", this->m_Cmu); + pp.query("Cmu_prime", this->m_Cmu_prime); + pp.query("Cb_stable", this->m_Cb_stable); + pp.query("Cb_unstable", this->m_Cb_unstable); + pp.query("prandtl", this->m_prandtl); +} + +template +TurbulenceModel::CoeffsDictType KLAxell::model_coeffs() const +{ + return TurbulenceModel::CoeffsDictType{ + {"Cmu", this->m_Cmu}, + {"Cmu_prime", this->m_Cmu_prime}, + {"Cb_stable", this->m_Cb_stable}, + {"Cb_unstable", this->m_Cb_unstable}, + {"prandtl", this->m_prandtl}}; +} + +template +void KLAxell::update_turbulent_viscosity( + const FieldState fstate, const DiffusionType /*unused*/) +{ + BL_PROFILE( + "amr-wind::" + this->identifier() + "::update_turbulent_viscosity"); + + auto gradT = (this->m_sim.repo()).create_scratch_field(3, 0); + fvm::gradient(*gradT, m_temperature.state(fstate)); + + const auto& vel = this->m_vel.state(fstate); + fvm::strainrate(this->m_shear_prod, vel); + + const amrex::GpuArray gravity{ + m_gravity[0], m_gravity[1], m_gravity[2]}; + const amrex::Real beta = 1.0 / m_ref_theta; + const amrex::Real Cmu = m_Cmu; + const amrex::Real Cb_stable = m_Cb_stable; + const amrex::Real Cb_unstable = m_Cb_unstable; + auto& mu_turb = this->mu_turb(); + const auto& den = this->m_rho.state(fstate); + const auto& repo = mu_turb.repo(); + const auto& geom_vec = repo.mesh().Geom(); + const int nlevels = repo.num_active_levels(); + const amrex::Real Rtc = -1.0; + const amrex::Real Rtmin = -3.0; + const amrex::Real lambda = 30.0; + const amrex::Real kappa = 0.41; + const amrex::Real surf_flux = m_surf_flux; + const auto tiny = std::numeric_limits::epsilon(); + const amrex::Real lengthscale_switch = m_lengthscale_switch; + for (int lev = 0; lev < nlevels; ++lev) { + const auto& geom = geom_vec[lev]; + const auto& problo = repo.mesh().Geom(lev).ProbLoArray(); + const amrex::Real dz = geom.CellSize()[2]; + for (amrex::MFIter mfi(mu_turb(lev)); mfi.isValid(); ++mfi) { + const auto& bx = mfi.tilebox(); + const auto& mu_arr = mu_turb(lev).array(mfi); + const auto& rho_arr = den(lev).const_array(mfi); + const auto& gradT_arr = (*gradT)(lev).array(mfi); + const auto& tlscale_arr = (this->m_turb_lscale)(lev).array(mfi); + 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); + //! Add terrain components + const bool has_terrain = + this->m_sim.repo().int_field_exists("terrain_blank"); + if (has_terrain) { + const auto* m_terrain_height = + &this->m_sim.repo().get_field("terrain_height"); + const auto* m_terrain_blank = + &this->m_sim.repo().get_int_field("terrain_blank"); + const auto& ht_arr = (*m_terrain_height)(lev).const_array(mfi); + const auto& blank_arr = + (*m_terrain_blank)(lev).const_array(mfi); + amrex::ParallelFor( + bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + amrex::Real stratification = + -(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; + const amrex::Real z = std::max( + problo[2] + (k + 0.5) * dz - ht_arr(i, j, k), + 0.5 * dz); + const amrex::Real lscale_s = + (lambda * kappa * z) / (lambda + kappa * z); + const amrex::Real lscale_b = + Cb_stable * std::sqrt( + tke_arr(i, j, k) / + std::max(stratification, tiny)); + amrex::Real epsilon = std::pow(Cmu, 3) * + std::pow(tke_arr(i, j, k), 1.5) / + (tlscale_arr(i, j, k) + tiny); + amrex::Real Rt = + std::pow(tke_arr(i, j, k) / epsilon, 2) * + stratification; + Rt = (Rt > Rtc) + ? Rt + : std::max( + Rt, Rt - std::pow(Rt - Rtc, 2) / + (Rt + Rtmin - 2 * Rtc)); + tlscale_arr(i, j, k) = + (stratification > 0) + ? std::sqrt( + std::pow(lscale_s * lscale_b, 2) / + (std::pow(lscale_s, 2) + + std::pow(lscale_b, 2))) + : lscale_s * + std::sqrt( + 1.0 - + std::pow(Cmu, 6.0) * + std::pow(Cb_unstable, -2.0) * Rt); + tlscale_arr(i, j, k) = + (stratification > 0) + ? std::min( + tlscale_arr(i, j, k), + std::sqrt( + Cmu * tke_arr(i, j, k) / + stratification)) + : tlscale_arr(i, j, k); + tlscale_arr(i, j, k) = (std::abs(surf_flux) < 1e-5 && + z <= lengthscale_switch) + ? lscale_s + : tlscale_arr(i, j, k); + Rt = (std::abs(surf_flux) < 1e-5 && + z <= lengthscale_switch) + ? 0.0 + : Rt; + const amrex::Real Cmu_Rt = + (Cmu + 0.108 * Rt) / + (1 + 0.308 * Rt + 0.00837 * std::pow(Rt, 2)); + mu_arr(i, j, k) = rho_arr(i, j, k) * Cmu_Rt * + tlscale_arr(i, j, k) * + std::sqrt(tke_arr(i, j, k)) * + (1 - blank_arr(i, j, k)); + const amrex::Real Cmu_prime_Rt = Cmu / (1 + 0.277 * Rt); + const amrex::Real muPrime = + rho_arr(i, j, k) * Cmu_prime_Rt * + tlscale_arr(i, j, k) * std::sqrt(tke_arr(i, j, k)) * + (1 - blank_arr(i, j, k)); + buoy_prod_arr(i, j, k) = -muPrime * stratification; + shear_prod_arr(i, j, k) *= + shear_prod_arr(i, j, k) * mu_arr(i, j, k); + }); + } else { + amrex::ParallelFor( + bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + amrex::Real stratification = + -(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; + const amrex::Real z = problo[2] + (k + 0.5) * dz; + const amrex::Real lscale_s = + (lambda * kappa * z) / (lambda + kappa * z); + const amrex::Real lscale_b = + Cb_stable * std::sqrt( + tke_arr(i, j, k) / + std::max(stratification, tiny)); + amrex::Real epsilon = std::pow(Cmu, 3) * + std::pow(tke_arr(i, j, k), 1.5) / + (tlscale_arr(i, j, k) + tiny); + amrex::Real Rt = + std::pow(tke_arr(i, j, k) / epsilon, 2) * + stratification; + Rt = (Rt > Rtc) + ? Rt + : std::max( + Rt, Rt - std::pow(Rt - Rtc, 2) / + (Rt + Rtmin - 2 * Rtc)); + tlscale_arr(i, j, k) = + (stratification > 0) + ? std::sqrt( + std::pow(lscale_s * lscale_b, 2) / + (std::pow(lscale_s, 2) + + std::pow(lscale_b, 2))) + : lscale_s * + std::sqrt( + 1.0 - + std::pow(Cmu, 6.0) * + std::pow(Cb_unstable, -2.0) * Rt); + tlscale_arr(i, j, k) = + (stratification > 0) + ? std::min( + tlscale_arr(i, j, k), + std::sqrt( + Cmu * tke_arr(i, j, k) / + stratification)) + : tlscale_arr(i, j, k); + tlscale_arr(i, j, k) = (std::abs(surf_flux) < 1e-5 && + z <= lengthscale_switch) + ? lscale_s + : tlscale_arr(i, j, k); + Rt = (std::abs(surf_flux) < 1e-5 && + z <= lengthscale_switch) + ? 0.0 + : Rt; + const amrex::Real Cmu_Rt = + (Cmu + 0.108 * Rt) / + (1 + 0.308 * Rt + 0.00837 * std::pow(Rt, 2)); + mu_arr(i, j, k) = rho_arr(i, j, k) * Cmu_Rt * + tlscale_arr(i, j, k) * + std::sqrt(tke_arr(i, j, k)); + const amrex::Real Cmu_prime_Rt = Cmu / (1 + 0.277 * Rt); + const amrex::Real muPrime = + rho_arr(i, j, k) * Cmu_prime_Rt * + tlscale_arr(i, j, k) * std::sqrt(tke_arr(i, j, k)); + buoy_prod_arr(i, j, k) = -muPrime * stratification; + shear_prod_arr(i, j, k) *= + shear_prod_arr(i, j, k) * mu_arr(i, j, k); + }); + } + } + } + + mu_turb.fillpatch(this->m_sim.time().current_time()); +} + +template +void KLAxell::update_alphaeff(Field& alphaeff) +{ + + BL_PROFILE("amr-wind::" + this->identifier() + "::update_alphaeff"); + auto lam_alpha = (this->m_transport).alpha(); + auto& mu_turb = this->m_mu_turb; + auto& repo = mu_turb.repo(); + auto gradT = (this->m_sim.repo()).create_scratch_field(3, 0); + 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 amrex::Real Cmu = m_Cmu; + const int nlevels = repo.num_active_levels(); + for (int lev = 0; lev < nlevels; ++lev) { + for (amrex::MFIter mfi(mu_turb(lev)); mfi.isValid(); ++mfi) { + const auto& bx = mfi.tilebox(); + const auto& muturb_arr = mu_turb(lev).array(mfi); + const auto& alphaeff_arr = alphaeff(lev).array(mfi); + const auto& lam_diff_arr = (*lam_alpha)(lev).array(mfi); + 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 amrex::Real Rtc = -1.0; + const amrex::Real Rtmin = -3.0; + amrex::ParallelFor( + bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + amrex::Real stratification = + -(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; + amrex::Real epsilon = std::pow(Cmu, 3) * + std::pow(tke_arr(i, j, k), 1.5) / + tlscale_arr(i, j, k); + amrex::Real Rt = std::pow(tke_arr(i, j, k) / epsilon, 2) * + stratification; + Rt = (Rt > Rtc) ? Rt + : std::max( + Rt, Rt - std::pow(Rt - Rtc, 2) / + (Rt + Rtmin - 2 * Rtc)); + const amrex::Real prandtlRt = + (1 + 0.193 * Rt) / (1 + 0.0302 * Rt); + alphaeff_arr(i, j, k) = + lam_diff_arr(i, j, k) + muturb_arr(i, j, k) / prandtlRt; + }); + } + } + + alphaeff.fillpatch(this->m_sim.time().current_time()); +} + +template +void KLAxell::update_scalar_diff( + Field& deff, const std::string& name) +{ + + BL_PROFILE("amr-wind::" + this->identifier() + "::update_scalar_diff"); + + if (name == pde::TKE::var_name()) { + auto& mu_turb = this->mu_turb(); + deff.setVal(0.0); + field_ops::saxpy( + deff, 2.0, mu_turb, 0, 0, deff.num_comp(), deff.num_grow()); + } else { + amrex::Abort( + "KLAxell:update_scalar_diff not implemented for field " + name); + } +} + +template +void KLAxell::post_advance_work() +{ + + BL_PROFILE("amr-wind::" + this->identifier() + "::post_advance_work"); +} + +} // namespace turbulence + +INSTANTIATE_TURBULENCE_MODEL(KLAxell); + +} // namespace amr_wind diff --git a/amr-wind/wind_energy/ABLFieldInit.H b/amr-wind/wind_energy/ABLFieldInit.H index 0fe98676b6..284e688456 100644 --- a/amr-wind/wind_energy/ABLFieldInit.H +++ b/amr-wind/wind_energy/ABLFieldInit.H @@ -65,18 +65,25 @@ private: amrex::Vector m_theta_heights; amrex::Vector m_theta_values; + //! Adding option for wind heights similar to temperature heights + //! Speed-up RANS calculation using 1-D profile for flat surface + bool m_initial_wind_profile{false}; + //! File name for 1-D data file + std::string m_1d_rans; + amrex::Vector m_wind_heights; amrex::Vector m_u_values; amrex::Vector m_v_values; - + amrex::Vector m_tke_values; ///@} // Device copies of the above arrays amrex::Gpu::DeviceVector m_thht_d; amrex::Gpu::DeviceVector m_thvv_d; + amrex::Gpu::DeviceVector m_windht_d; amrex::Gpu::DeviceVector m_prof_u_d; amrex::Gpu::DeviceVector m_prof_v_d; - + amrex::Gpu::DeviceVector m_prof_tke_d; //! Initial density field amrex::Real m_rho; diff --git a/amr-wind/wind_energy/ABLFieldInit.cpp b/amr-wind/wind_energy/ABLFieldInit.cpp index bbddad29b0..53341c205a 100644 --- a/amr-wind/wind_energy/ABLFieldInit.cpp +++ b/amr-wind/wind_energy/ABLFieldInit.cpp @@ -7,7 +7,7 @@ #include "amr-wind/utilities/trig_ops.H" #include "AMReX_Gpu.H" #include "AMReX_ParmParse.H" - +#include "amr-wind/utilities/linear_interpolation.H" namespace amr_wind { ABLFieldInit::ABLFieldInit() @@ -23,7 +23,24 @@ ABLFieldInit::ABLFieldInit() void ABLFieldInit::initialize_from_inputfile() { amrex::ParmParse pp_abl("ABL"); - + //! Check for wind profile + pp_abl.query("initial_wind_profile", m_initial_wind_profile); + if (m_initial_wind_profile) { + pp_abl.query("rans_1dprofile_file", m_1d_rans); + if (!m_1d_rans.empty()) { + std::ifstream ransfile(m_1d_rans, std::ios::in); + if (!ransfile.good()) { + amrex::Abort("Cannot find 1-D RANS profile file " + m_1d_rans); + } + amrex::Real value1, value2, value3, value4; + while (ransfile >> value1 >> value2 >> value3 >> value4) { + m_wind_heights.push_back(value1); + m_u_values.push_back(value2); + m_v_values.push_back(value3); + m_tke_values.push_back(value4); + } + } + } // Temperature variation as a function of height pp_abl.getarr("temperature_heights", m_theta_heights); pp_abl.getarr("temperature_values", m_theta_values); @@ -61,7 +78,8 @@ void ABLFieldInit::initialize_from_inputfile() pp_incflo.get("density", m_rho); } else { pp_mphase.get("density_fluid2", m_rho); - // Note: density field will later be overwritten by MultiPhase post_init + // Note: density field will later be overwritten by MultiPhase + // post_init } amrex::ParmParse pp_forcing("ABLForcing"); @@ -84,7 +102,25 @@ void ABLFieldInit::initialize_from_inputfile() m_thht_d.resize(num_theta_values); m_thvv_d.resize(num_theta_values); - + if (m_initial_wind_profile) { + int num_wind_values = static_cast(m_wind_heights.size()); + m_windht_d.resize(num_wind_values); + m_prof_u_d.resize(num_wind_values); + m_prof_v_d.resize(num_wind_values); + m_prof_tke_d.resize(num_wind_values); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, m_wind_heights.begin(), + m_wind_heights.end(), m_windht_d.begin()); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, m_u_values.begin(), m_u_values.end(), + m_prof_u_d.begin()); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, m_v_values.begin(), m_v_values.end(), + m_prof_v_d.begin()); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, m_tke_values.begin(), m_tke_values.end(), + m_prof_tke_d.begin()); + } amrex::Gpu::copy( amrex::Gpu::hostToDevice, m_theta_heights.begin(), m_theta_heights.end(), m_thht_d.begin()); @@ -158,6 +194,7 @@ void ABLFieldInit::operator()( const amrex::Real rho_init = m_rho; const int ntvals = static_cast(m_theta_heights.size()); + const int nwvals = static_cast(m_wind_heights.size()); const amrex::Real* th = m_thht_d.data(); const amrex::Real* tv = m_thvv_d.data(); @@ -193,13 +230,39 @@ void ABLFieldInit::operator()( } } + temperature(i, j, k, 0) += theta; + velocity(i, j, k, 0) += umean_prof; + velocity(i, j, k, 1) += vmean_prof; + }); + } else if (m_initial_wind_profile) { + //! RANS 1-D profile + const amrex::Real* windh = m_windht_d.data(); + const amrex::Real* uu = m_prof_u_d.data(); + const amrex::Real* vv = m_prof_v_d.data(); + + amrex::ParallelFor( + vbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + const amrex::Real z = problo[2] + (k + 0.5) * dx[2]; + + density(i, j, k) = rho_init; + const amrex::Real theta = + (ntvals > 0) ? interp::linear(th, th + ntvals, tv, z) + : tv[0]; + const amrex::Real umean_prof = + (nwvals > 0) ? interp::linear(windh, windh + nwvals, uu, z) + : uu[0]; + const amrex::Real vmean_prof = + (nwvals > 0) ? interp::linear(windh, windh + nwvals, vv, z) + : vv[0]; + temperature(i, j, k, 0) += theta; velocity(i, j, k, 0) += umean_prof; velocity(i, j, k, 1) += vmean_prof; }); } else { /* - * Set uniform/linear wind profile with specified temperature profile + * Set uniform/linear wind profile with specified temperature + * profile */ const bool linear_profile = m_linear_profile; @@ -254,8 +317,9 @@ void ABLFieldInit::operator()( }); } - // velocity perturbations may be added on top of the simple wind profiles - // specified in the input file or the general profiles from a netcdf input + // velocity perturbations may be added on top of the simple wind + // profiles specified in the input file or the general profiles from + // a netcdf input if (m_perturb_vel) { const amrex::Real aval = m_Uperiods * 2.0 * pi / (probhi[1] - problo[1]); @@ -324,7 +388,7 @@ void ABLFieldInit::init_tke( { tke_fab.setVal(m_tke_init); - if (!m_tke_init_profile) { + if (!m_tke_init_profile && !m_initial_wind_profile) { return; } @@ -340,19 +404,35 @@ void ABLFieldInit::init_tke( ++mfi) { const auto& bx = mfi.tilebox(); const auto& tke = tke_fab.array(mfi); - - // Profile definition from Beare et al. (2006) - amrex::ParallelFor( - bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - const amrex::Real z = problo[2] + (k + 0.5) * dx[2]; - - if (z < tke_cutoff_height) { - tke(i, j, k) = tke_init_factor * - std::pow(1. - z / tke_cutoff_height, 3); - } else { - tke(i, j, k) = 1.e-20; - } - }); + const auto tiny = std::numeric_limits::epsilon(); + if (m_initial_wind_profile) { + const int nwvals = static_cast(m_wind_heights.size()); + const amrex::Real* windh = m_windht_d.data(); + const amrex::Real* tke_data = m_prof_tke_d.data(); + amrex::ParallelFor( + bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + const amrex::Real z = problo[2] + (k + 0.5) * dx[2]; + const amrex::Real tke_prof = + (nwvals > 0) + ? interp::linear(windh, windh + nwvals, tke_data, z) + : tiny; + + tke(i, j, k) = tke_prof; + }); + } else { + // Profile definition from Beare et al. (2006) + amrex::ParallelFor( + bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + const amrex::Real z = problo[2] + (k + 0.5) * dx[2]; + + if (z < tke_cutoff_height) { + tke(i, j, k) = tke_init_factor * + std::pow(1. - z / tke_cutoff_height, 3); + } else { + tke(i, j, k) = tiny; + } + }); + } } } diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 7c474209b8..8da6d177b2 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -283,6 +283,7 @@ if(AMR_WIND_ENABLE_NETCDF) add_test_re(abl_sampling_netcdf) add_test_re(abl_kosovic_neutral_ib) add_test_re(nrel_precursor) + add_test_re(abl_wallrans_neutral) endif() if(AMR_WIND_ENABLE_MASA) diff --git a/test/test_files/abl_wallrans_neutral/abl_wallrans_neutral.inp b/test/test_files/abl_wallrans_neutral/abl_wallrans_neutral.inp new file mode 100644 index 0000000000..39f29cb9ea --- /dev/null +++ b/test/test_files/abl_wallrans_neutral/abl_wallrans_neutral.inp @@ -0,0 +1,86 @@ +# Geometry +geometry.prob_lo = 0.0 0.0 0.0 +geometry.prob_hi = 2048 2048 1024 +geometry.is_periodic = 1 1 0 +# Grid +amr.n_cell = 64 64 64 +amr.max_level = 0 +# Simulation control parameters +time.max_step = 10 +time.initial_dt = 0.5 +time.cfl = 0.5 +time.fixed_dt = -1 +time.init_shrink = 0.1 +time.regrid_interval = -1 +time.plot_interval = 86400 +time.checkpoint_interval = 86400 +#io +#io.restart_file = "chk72000" +io.output_default_variables = false +io.outputs = velocity temperature tke mu_turb +# incflo +incflo.physics = ABL +incflo.density = 1.225 +incflo.gravity = 0. 0. -9.81 # Gravitational force (3D) +incflo.velocity = 10.0 0.0 0.0 +incflo.verbose = 0 +incflo.initial_iterations = 8 +incflo.do_initial_proj = true +incflo.constant_density = true +incflo.use_godunov = true +incflo.godunov_type = "weno_z" +incflo.diffusion_type = 2 +# transport equation parameters +transport.model = ConstTransport +transport.viscosity = 0.0 +transport.laminar_prandtl = 0.7 +transport.turbulent_prandtl = 0.333 +# turbulence equation parameters +turbulence.model = KLAxell +TKE.source_terms = KransAxell +# Atmospheric boundary layer +ABL.Uperiods = 12.0 +ABL.Vperiods = 12.0 +ABL.cutoff_height = 50.0 +ABL.deltaU = 0 +ABL.deltaV = 0 +ABL.kappa = .41 +ABL.normal_direction = 2 +ABL.perturb_ref_height = 50.0 +ABL.perturb_velocity = true +ABL.perturb_temperature = false +ABL.reference_temperature = 300. +ABL.stats_output_format = netcdf +ABL.surface_roughness_z0 = 0.1 +ABL.initial_wind_profile = true +ABL.rans_1dprofile_file = "rans_1d.info" +ABL.temperature_heights = 0 800 900 1800 2700 +ABL.temperature_values = 300 300 308 311 314 +ABL.wall_shear_stress_type = local +ABL.surface_temp_flux = 0.0 +# Source +ICNS.source_terms = BoussinesqBuoyancy CoriolisForcing GeostrophicForcing RayleighDamping +BoussinesqBuoyancy.reference_temperature = 300.0 +BoussinesqBuoyancy.thermal_expansion_coeff = 0.003333 +#CoriolisForcing.east_vector = 1.0 0.0 0.0 +#CoriolisForcing.latitude = 45.0 +#CoriolisForcing.north_vector = 0.0 1.0 0.0 +#CoriolisForcing.rotational_time_period = 86400.0 +CoriolisForcing.east_vector = 1.0 0.0 0.0 +CoriolisForcing.north_vector = 0.0 1.0 0.0 +CoriolisForcing.latitude = 90.0 +CoriolisForcing.rotational_time_period = 125663.706143592 +GeostrophicForcing.geostrophic_wind = 10 0.0 0.0 +RayleighDamping.reference_velocity = 10 0.0 0.0 +RayleighDamping.length_sloped_damping = 150 +RayleighDamping.length_complete_damping = 50 +RayleighDamping.time_scale = 20.0 +# BC +zhi.type = "slip_wall" +zhi.temperature_type = "fixed_gradient" +zhi.temperature = 0.0 +zhi.tke_type = "fixed_gradient" +zlo.type = "wall_model" + + + diff --git a/test/test_files/abl_wallrans_neutral/rans_1d.info b/test/test_files/abl_wallrans_neutral/rans_1d.info new file mode 100644 index 0000000000..7148cd774c --- /dev/null +++ b/test/test_files/abl_wallrans_neutral/rans_1d.info @@ -0,0 +1,5 @@ +0 3.6715 0.95701 0.894618 +502.618 15 1.86236e-10 6.6423e-08 +507.853 15 0 1.01178e-14 +952.88 15 0 1e-15 +1000 15 0 1.83246e-08 \ No newline at end of file diff --git a/unit_tests/turbulence/CMakeLists.txt b/unit_tests/turbulence/CMakeLists.txt index dc87c3adf4..0baafb7bdc 100644 --- a/unit_tests/turbulence/CMakeLists.txt +++ b/unit_tests/turbulence/CMakeLists.txt @@ -3,5 +3,6 @@ target_sources(${amr_wind_unit_test_exe_name} test_turbulence_init.cpp test_turbulence_LES.cpp + test_turbulence_RANS.cpp test_turbulence_LES_bc.cpp ) diff --git a/unit_tests/turbulence/test_turbulence_RANS.cpp b/unit_tests/turbulence/test_turbulence_RANS.cpp new file mode 100644 index 0000000000..8afd34310c --- /dev/null +++ b/unit_tests/turbulence/test_turbulence_RANS.cpp @@ -0,0 +1,184 @@ +#include "gtest/gtest.h" +#include "aw_test_utils/MeshTest.H" +#include "amr-wind/turbulence/TurbulenceModel.H" +#include "aw_test_utils/test_utils.H" + +namespace amr_wind_tests { + +namespace { + +void init_strain_field(amr_wind::Field& fld, amrex::Real srate) +{ + const auto& mesh = fld.repo().mesh(); + const int nlevels = fld.repo().num_active_levels(); + amrex::Real offset = + (fld.field_location() == amr_wind::FieldLoc::CELL) ? 0.5 : 0.0; + for (int lev = 0; lev < nlevels; ++lev) { + const auto& dx = mesh.Geom(lev).CellSizeArray(); + const auto& problo = mesh.Geom(lev).ProbLoArray(); + + for (amrex::MFIter mfi(fld(lev)); mfi.isValid(); ++mfi) { + auto bx = mfi.growntilebox(); + const auto& farr = fld(lev).array(mfi); + + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + const amrex::Real x = problo[0] + (i + offset) * dx[0]; + const amrex::Real y = problo[1] + (j + offset) * dx[1]; + const amrex::Real z = problo[2] + (k + offset) * dx[2]; + + farr(i, j, k, 0) = x / sqrt(6.0) * srate; + farr(i, j, k, 1) = y / sqrt(6.0) * srate; + farr(i, j, k, 2) = z / sqrt(6.0) * srate; + }); + } + } +} + +void init_temperature_field(amr_wind::Field& fld, amrex::Real tgrad) +{ + const auto& mesh = fld.repo().mesh(); + const int nlevels = fld.repo().num_active_levels(); + + amrex::Real offset = + (fld.field_location() == amr_wind::FieldLoc::CELL) ? 0.5 : 0.0; + + for (int lev = 0; lev < nlevels; ++lev) { + const auto& dx = mesh.Geom(lev).CellSizeArray(); + const auto& problo = mesh.Geom(lev).ProbLoArray(); + + for (amrex::MFIter mfi(fld(lev)); mfi.isValid(); ++mfi) { + auto bx = mfi.growntilebox(); + const auto& farr = fld(lev).array(mfi); + + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + const amrex::Real z = problo[2] + (k + offset) * dx[2]; + + farr(i, j, k, 0) = z * tgrad; + }); + } + } +} + +} // namespace + +class TurbRANSTest : public MeshTest +{ +protected: + void populate_parameters() override + { + MeshTest::populate_parameters(); + + { + amrex::ParmParse pp("amr"); + amrex::Vector ncell{{10, 10, 64}}; + pp.addarr("n_cell", ncell); + pp.add("blocking_factor", 2); + } + { + amrex::ParmParse pp("geometry"); + amrex::Vector problo{{0.0, 0.0, 0.0}}; + amrex::Vector probhi{{1024.0, 1024.0, 1024.0}}; + pp.addarr("prob_lo", problo); + pp.addarr("prob_hi", probhi); + } + } +}; + +TEST_F(TurbRANSTest, test_1eqKrans_setup_calc) +{ + // Parser inputs for turbulence model + const amrex::Real Tref = 265; + const amrex::Real gravz = 10.0; + const amrex::Real rho0 = 1.2; + { + amrex::ParmParse pp("turbulence"); + pp.add("model", (std::string) "KLAxell"); + } + { + amrex::ParmParse pp("incflo"); + amrex::Vector physics{"ABL"}; + pp.addarr("physics", physics); + pp.add("density", rho0); + amrex::Vector vvec{8.0, 0.0, 0.0}; + pp.addarr("velocity", vvec); + amrex::Vector gvec{0.0, 0.0, -gravz}; + pp.addarr("gravity", gvec); + } + { + amrex::ParmParse pp("ABL"); + pp.add("reference_temperature", Tref); + pp.add("surface_temp_rate", 0.0); + pp.add("initial_wind_profile", true); + amrex::Vector t_hts{0.0, 100.0, 4000.0}; + pp.addarr("temperature_heights", t_hts); + pp.addarr("wind_heights", t_hts); + amrex::Vector t_vals{265.0, 265.0, 265.0}; + pp.addarr("temperature_values", t_vals); + amrex::Vector u_vals{8.0, 8.0, 8.0}; + pp.addarr("u_values", u_vals); + amrex::Vector v_vals{0.0, 0.0, 0.0}; + pp.addarr("v_values", v_vals); + amrex::Vector tke_vals{0.1, 0.1, 0.1}; + pp.addarr("tke_values", tke_vals); + pp.add("surface_temp_flux", 0.0); + } + + // Initialize necessary parts of solver + populate_parameters(); + initialize_mesh(); + auto& pde_mgr = sim().pde_manager(); + pde_mgr.register_icns(); + sim().init_physics(); + + // Create turbulence model + sim().create_turbulence_model(); + // Get turbulence model + auto& tmodel = sim().turbulence_model(); + + // Get coefficients + auto model_dict = tmodel.model_coeffs(); + + // Constants for fields + const amrex::Real srate = 20.0; + const amrex::Real Tgz = 0.0; + const amrex::Real lambda = 30.0; + const amrex::Real kappa = 0.41; + const amrex::Real x3 = 1016; + const amrex::Real lscale_s = (lambda * kappa * x3) / (lambda + kappa * x3); + const amrex::Real tlscale_val = lscale_s; + const amrex::Real tke_val = 0.1; + // Set up velocity field with constant strainrate + auto& vel = sim().repo().get_field("velocity"); + init_strain_field(vel, srate); + // Set up uniform unity density field + auto& dens = sim().repo().get_field("density"); + dens.setVal(rho0); + // Set up temperature field with constant gradient in z + auto& temp = sim().repo().get_field("temperature"); + init_temperature_field(temp, Tgz); + // Give values to tlscale and tke arrays + auto& tlscale = sim().repo().get_field("turb_lscale"); + tlscale.setVal(tlscale_val); + auto& tke = sim().repo().get_field("tke"); + tke.setVal(tke_val); + + // Update turbulent viscosity directly + tmodel.update_turbulent_viscosity( + amr_wind::FieldState::New, DiffusionType::Crank_Nicolson); + const auto& muturb = sim().repo().get_field("mu_turb"); + + // Check values of turbulent viscosity + const auto max_val = utils::field_max(muturb); + const amrex::Real Cmu = 0.556; + const amrex::Real epsilon = + std::pow(Cmu, 3) * std::pow(tke_val, 1.5) / (tlscale_val + 1e-3); + const amrex::Real stratification = 0.0; + const amrex::Real Rt = std::pow(tke_val / epsilon, 2) * stratification; + const amrex::Real Cmu_Rt = + (0.556 + 0.108 * Rt) / (1 + 0.308 * Rt + 0.00837 * std::pow(Rt, 2)); + const amrex::Real tol = 0.12; + const amrex::Real nut_max = rho0 * Cmu_Rt * tlscale_val * sqrt(tke_val); + EXPECT_NEAR(max_val, nut_max, tol); +} + +} // namespace amr_wind_tests