From afeede7d0cfe3a94173585c076740d4bc1136d66 Mon Sep 17 00:00:00 2001 From: Harish Date: Sun, 5 Jan 2025 16:49:39 -0700 Subject: [PATCH] Compatible with mesoscale zone --- .../icns/source_terms/CMakeLists.txt | 1 + .../icns/source_terms/DragForcing.cpp | 14 +- .../icns/source_terms/WindSpongeForcing.H | 45 ++++++ .../icns/source_terms/WindSpongeForcing.cpp | 153 ++++++++++++++++++ .../source_terms/TempSpongeForcing.H | 15 +- .../source_terms/TempSpongeForcing.cpp | 127 +++++++++++++-- .../tke/source_terms/KransAxell.H | 20 ++- .../tke/source_terms/KransAxell.cpp | 152 ++++++++++++++--- amr-wind/physics/TerrainDrag.cpp | 3 +- amr-wind/turbulence/RANS/KLAxell.H | 2 +- amr-wind/turbulence/RANS/KLAxell.cpp | 4 +- 11 files changed, 486 insertions(+), 50 deletions(-) create mode 100644 amr-wind/equation_systems/icns/source_terms/WindSpongeForcing.H create mode 100644 amr-wind/equation_systems/icns/source_terms/WindSpongeForcing.cpp diff --git a/amr-wind/equation_systems/icns/source_terms/CMakeLists.txt b/amr-wind/equation_systems/icns/source_terms/CMakeLists.txt index 1d9c83c943..7a03eecae5 100644 --- a/amr-wind/equation_systems/icns/source_terms/CMakeLists.txt +++ b/amr-wind/equation_systems/icns/source_terms/CMakeLists.txt @@ -16,4 +16,5 @@ target_sources(${amr_wind_lib_name} PRIVATE NonLinearSGSTerm.cpp DragForcing.cpp ForestForcing.cpp + WindSpongeForcing.cpp ) diff --git a/amr-wind/equation_systems/icns/source_terms/DragForcing.cpp b/amr-wind/equation_systems/icns/source_terms/DragForcing.cpp index 8939b70519..aa10761d30 100644 --- a/amr-wind/equation_systems/icns/source_terms/DragForcing.cpp +++ b/amr-wind/equation_systems/icns/source_terms/DragForcing.cpp @@ -118,28 +118,28 @@ void DragForcing::operator()( yi_end = sponge_north * std::max(yi_end, 0.0); ystart_damping = sponge_strength * yi_start * yi_start; yend_damping = sponge_strength * yi_end * yi_end; - + const amrex::Real ux1 = vel(i, j, k, 0); + const amrex::Real uy1 = vel(i, j, k, 1); + const amrex::Real uz1 = vel(i, j, k, 2); const auto idx = interp::bisection_search(device_vel_ht, device_vel_ht + vsize, z); const amrex::Real spongeVelX = (vsize > 0) ? interp::linear_impl( device_vel_ht, device_vel_vals, z, idx, 3, 0) - : 0.0; + : ux1; const amrex::Real spongeVelY = (vsize > 0) ? interp::linear_impl( device_vel_ht, device_vel_vals, z, idx, 3, 1) - : 0.0; + : uy1; const amrex::Real spongeVelZ = (vsize > 0) ? interp::linear_impl( device_vel_ht, device_vel_vals, z, idx, 3, 2) - : 0.0; + : uz1; amrex::Real Dxz = 0.0; amrex::Real Dyz = 0.0; amrex::Real bc_forcing_x = 0; amrex::Real bc_forcing_y = 0; - const amrex::Real ux1 = vel(i, j, k, 0); - const amrex::Real uy1 = vel(i, j, k, 1); - const amrex::Real uz1 = vel(i, j, k, 2); + const amrex::Real m = std::sqrt(ux1 * ux1 + uy1 * uy1 + uz1 * uz1); if (drag(i, j, k) == 1 && (!is_laminar)) { const amrex::Real ux2 = vel(i, j, k + 1, 0); diff --git a/amr-wind/equation_systems/icns/source_terms/WindSpongeForcing.H b/amr-wind/equation_systems/icns/source_terms/WindSpongeForcing.H new file mode 100644 index 0000000000..24ad0bfd13 --- /dev/null +++ b/amr-wind/equation_systems/icns/source_terms/WindSpongeForcing.H @@ -0,0 +1,45 @@ +#ifndef WINDSPONGEFORCING_H +#define WINDSPONGEFORCING_H + +#include "amr-wind/equation_systems/icns/MomentumSource.H" +#include "amr-wind/core/SimTime.H" +#include "amr-wind/CFDSim.H" + +namespace amr_wind::pde::icns { + +class WindSpongeForcing : public MomentumSource::Register +{ +public: + static std::string identifier() { return "WindSpongeForcing"; } + + explicit WindSpongeForcing(const CFDSim& sim); + + ~WindSpongeForcing() 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 SimTime& m_time; + const amrex::AmrCore& m_mesh; + std::string m_1d_rans; + const Field& m_velocity; + amrex::Vector m_wind_heights; + amrex::Vector m_u_values; + amrex::Vector m_v_values; + amrex::Vector m_w_values; + amrex::Gpu::DeviceVector m_wind_heights_d; + amrex::Gpu::DeviceVector m_u_values_d; + amrex::Gpu::DeviceVector m_v_values_d; + amrex::Gpu::DeviceVector m_w_values_d; + amrex::Real m_meso_start{600}; + amrex::Real m_meso_timescale{30}; + const CFDSim& m_sim; +}; + +} // namespace amr_wind::pde::icns +#endif \ No newline at end of file diff --git a/amr-wind/equation_systems/icns/source_terms/WindSpongeForcing.cpp b/amr-wind/equation_systems/icns/source_terms/WindSpongeForcing.cpp new file mode 100644 index 0000000000..4af83ba9e6 --- /dev/null +++ b/amr-wind/equation_systems/icns/source_terms/WindSpongeForcing.cpp @@ -0,0 +1,153 @@ + +#include "amr-wind/equation_systems/icns/source_terms/WindSpongeForcing.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::icns { + +WindSpongeForcing::WindSpongeForcing(const CFDSim& sim) + : m_time(sim.time()) + , m_mesh(sim.mesh()) + , m_velocity(sim.repo().get_field("velocity")) + , m_sim(sim) +{ + amrex::ParmParse pp_abl("ABL"); + 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, value5; + while (ransfile >> value1 >> value2 >> value3 >> value4 >> value5) { + m_wind_heights.push_back(value1); + m_u_values.push_back(value2); + m_v_values.push_back(value3); + m_w_values.push_back(value4); + } + } else { + amrex::Abort("Cannot find 1-D RANS profile file " + m_1d_rans); + } + pp_abl.query("meso_sponge_start", m_meso_start); + pp_abl.query("meso_timescale", m_meso_timescale); + int num_wind_values = static_cast(m_wind_heights.size()); + m_wind_heights_d.resize(num_wind_values); + m_u_values_d.resize(num_wind_values); + m_v_values_d.resize(num_wind_values); + m_w_values_d.resize(num_wind_values); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, m_wind_heights.begin(), m_wind_heights.end(), + m_wind_heights_d.begin()); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, m_u_values.begin(), m_u_values.end(), + m_u_values_d.begin()); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, m_v_values.begin(), m_v_values.end(), + m_v_values_d.begin()); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, m_w_values.begin(), m_w_values.end(), + m_w_values_d.begin()); +} + +WindSpongeForcing::~WindSpongeForcing() = default; + +void WindSpongeForcing::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& velocity = + m_velocity.state(field_impl::dof_state(fstate))(lev).const_array(mfi); + const amrex::Real sponge_start = m_meso_start; + const amrex::Real meso_timescale = m_time.delta_t(); + const auto vsize = m_wind_heights_d.size(); + const auto* wind_heights_d = m_wind_heights_d.data(); + const auto* u_values_d = m_u_values_d.data(); + const auto* v_values_d = m_v_values_d.data(); + const auto* w_values_d = m_w_values_d.data(); + const bool has_terrain = this->m_sim.repo().field_exists("terrain_height"); + if (has_terrain) { + auto* const m_terrain_height = + &this->m_sim.repo().get_field("terrain_height"); + const auto& terrain_height = (*m_terrain_height)(lev).const_array(mfi); + amrex::ParallelFor( + bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + const amrex::Real z = std::max( + prob_lo[2] + (k + 0.5) * dx[2] - terrain_height(i, j, k), + 0.5 * dx[2]); + const amrex::Real zi = std::max( + (z - sponge_start) / (prob_hi[2] - sponge_start), 0.0); + amrex::Real ref_windx = velocity(i, j, k, 0); + amrex::Real ref_windy = velocity(i, j, k, 1); + amrex::Real ref_windz = velocity(i, j, k, 2); + if (zi > 0) { + ref_windx = (vsize > 0) + ? interp::linear( + wind_heights_d, + wind_heights_d + vsize, u_values_d, z) + : velocity(i, j, k, 0); + ref_windy = (vsize > 0) + ? interp::linear( + wind_heights_d, + wind_heights_d + vsize, v_values_d, z) + : velocity(i, j, k, 1); + ref_windz = (vsize > 0) + ? interp::linear( + wind_heights_d, + wind_heights_d + vsize, w_values_d, z) + : velocity(i, j, k, 2); + } + src_term(i, j, k, 0) -= + 1.0 / meso_timescale * (velocity(i, j, k, 0) - ref_windx); + src_term(i, j, k, 1) -= + 1.0 / meso_timescale * (velocity(i, j, k, 1) - ref_windy); + src_term(i, j, k, 2) -= + 1.0 / meso_timescale * (velocity(i, j, k, 2) - ref_windz); + }); + } else { + 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_windx = velocity(i, j, k, 0); + amrex::Real ref_windy = velocity(i, j, k, 1); + amrex::Real ref_windz = velocity(i, j, k, 2); + if (zi > 0) { + ref_windx = (vsize > 0) + ? interp::linear( + wind_heights_d, + wind_heights_d + vsize, u_values_d, z) + : velocity(i, j, k, 0); + ref_windy = (vsize > 0) + ? interp::linear( + wind_heights_d, + wind_heights_d + vsize, v_values_d, z) + : velocity(i, j, k, 1); + ref_windz = (vsize > 0) + ? interp::linear( + wind_heights_d, + wind_heights_d + vsize, w_values_d, z) + : velocity(i, j, k, 2); + } + src_term(i, j, k, 0) -= + 1.0 / meso_timescale * (velocity(i, j, k, 0) - ref_windx); + src_term(i, j, k, 1) -= + 1.0 / meso_timescale * (velocity(i, j, k, 1) - ref_windy); + src_term(i, j, k, 2) -= + 1.0 / meso_timescale * (velocity(i, j, k, 2) - ref_windz); + }); + } +} + +} // namespace amr_wind::pde::icns diff --git a/amr-wind/equation_systems/temperature/source_terms/TempSpongeForcing.H b/amr-wind/equation_systems/temperature/source_terms/TempSpongeForcing.H index 494d6b7bac..093e20bd41 100644 --- a/amr-wind/equation_systems/temperature/source_terms/TempSpongeForcing.H +++ b/amr-wind/equation_systems/temperature/source_terms/TempSpongeForcing.H @@ -30,7 +30,20 @@ private: 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}; + amrex::Real m_meso_start{600}; + amrex::Real m_meso_timescale{30}; + bool m_horizontal_sponge{false}; + amrex::Real m_sponge_strength{1.0}; + amrex::Real m_sponge_density{1.0}; + amrex::Real m_sponge_distance_west{-1000}; + amrex::Real m_sponge_distance_east{1000}; + amrex::Real m_sponge_distance_south{-1000}; + amrex::Real m_sponge_distance_north{1000}; + int m_sponge_west{0}; + int m_sponge_east{1}; + int m_sponge_south{0}; + int m_sponge_north{1}; + const CFDSim& m_sim; }; } // namespace amr_wind::pde::temperature diff --git a/amr-wind/equation_systems/temperature/source_terms/TempSpongeForcing.cpp b/amr-wind/equation_systems/temperature/source_terms/TempSpongeForcing.cpp index 525b6dfef9..2ba7ef1541 100644 --- a/amr-wind/equation_systems/temperature/source_terms/TempSpongeForcing.cpp +++ b/amr-wind/equation_systems/temperature/source_terms/TempSpongeForcing.cpp @@ -10,13 +10,17 @@ namespace amr_wind::pde::temperature { TempSpongeForcing::TempSpongeForcing(const CFDSim& sim) - : m_mesh(sim.mesh()), m_temperature(sim.repo().get_field("temperature")) + : m_mesh(sim.mesh()) + , m_temperature(sim.repo().get_field("temperature")) + , m_sim(sim) { amrex::ParmParse pp_abl("ABL"); //! Temperature variation as a function of height - pp_abl.query("meso_sponge_start", m_sponge_start); + pp_abl.query("meso_sponge_start", m_meso_start); + pp_abl.query("meso_timescale", m_meso_timescale); pp_abl.getarr("temperature_heights", m_theta_heights); pp_abl.getarr("temperature_values", m_theta_values); + pp_abl.query("horizontal_sponge_temp", m_horizontal_sponge); 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); @@ -27,6 +31,17 @@ TempSpongeForcing::TempSpongeForcing(const CFDSim& sim) amrex::Gpu::copy( amrex::Gpu::hostToDevice, m_theta_values.begin(), m_theta_values.end(), m_theta_values_d.begin()); + amrex::ParmParse pp("DragForcing"); + pp.query("sponge_strength", m_sponge_strength); + pp.query("sponge_density", m_sponge_density); + pp.query("sponge_distance_west", m_sponge_distance_west); + pp.query("sponge_distance_east", m_sponge_distance_east); + pp.query("sponge_distance_south", m_sponge_distance_south); + pp.query("sponge_distance_north", m_sponge_distance_north); + pp.query("sponge_west", m_sponge_west); + pp.query("sponge_east", m_sponge_east); + pp.query("sponge_south", m_sponge_south); + pp.query("sponge_north", m_sponge_north); } TempSpongeForcing::~TempSpongeForcing() = default; @@ -45,24 +60,104 @@ void TempSpongeForcing::operator()( const auto& temperature = m_temperature.state(field_impl::dof_state(fstate))(lev).const_array( mfi); - const amrex::Real sponge_start = m_sponge_start; + const amrex::Real sponge_start = m_meso_start; + const amrex::Real meso_timescale = m_meso_timescale; 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); + const bool has_terrain = this->m_sim.repo().field_exists("terrain_height"); + if (has_terrain) { + auto* const m_terrain_height = + &this->m_sim.repo().get_field("terrain_height"); + const auto& terrain_height = (*m_terrain_height)(lev).const_array(mfi); + amrex::ParallelFor( + bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + const amrex::Real z = std::max( + prob_lo[2] + (k + 0.5) * dx[2] - terrain_height(i, j, 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) -= + 1.0 / meso_timescale * (temperature(i, j, k) - ref_temp); + }); + if (m_horizontal_sponge) { + const amrex::Real sponge_strength = m_sponge_strength; + const amrex::Real sponge_density = m_sponge_density; + const amrex::Real start_east = prob_hi[0] - m_sponge_distance_east; + const amrex::Real start_west = prob_lo[0] - m_sponge_distance_west; + const amrex::Real start_north = + prob_hi[1] - m_sponge_distance_north; + const amrex::Real start_south = + prob_lo[1] - m_sponge_distance_south; + const int sponge_east = m_sponge_east; + const int sponge_west = m_sponge_west; + const int sponge_south = m_sponge_south; + const int sponge_north = m_sponge_north; + amrex::ParallelFor( + bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + const amrex::Real x = prob_lo[0] + (i + 0.5) * dx[0]; + const amrex::Real y = prob_lo[1] + (j + 0.5) * dx[1]; + const amrex::Real z = prob_lo[2] + (k + 0.5) * dx[2]; + amrex::Real xstart_damping = 0; + amrex::Real ystart_damping = 0; + amrex::Real xend_damping = 0; + amrex::Real yend_damping = 0; + amrex::Real xi_end = + (x - start_east) / (prob_hi[0] - start_east); + amrex::Real xi_start = + (start_west - x) / (start_west - prob_lo[0]); + xi_start = sponge_west * std::max(xi_start, 0.0); + xi_end = sponge_east * std::max(xi_end, 0.0); + xstart_damping = + sponge_west * sponge_strength * xi_start * xi_start; + xend_damping = + sponge_east * sponge_strength * xi_end * xi_end; + amrex::Real yi_end = + (y - start_north) / (prob_hi[1] - start_north); + amrex::Real yi_start = + (start_south - y) / (start_south - prob_lo[1]); + yi_start = sponge_south * std::max(yi_start, 0.0); + yi_end = sponge_north * std::max(yi_end, 0.0); + ystart_damping = sponge_strength * yi_start * yi_start; + yend_damping = sponge_strength * yi_end * yi_end; + const amrex::Real 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) -= + (xstart_damping + xend_damping + ystart_damping + + yend_damping) * + (temperature(i, j, k) - sponge_density * ref_temp); + }); } - src_term(i, j, k, 0) -= zi * zi * (temperature(i, j, k) - ref_temp); - }); + } else { + 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) -= + 1.0 / meso_timescale * (temperature(i, j, k) - ref_temp); + }); + } } } // namespace amr_wind::pde::temperature diff --git a/amr-wind/equation_systems/tke/source_terms/KransAxell.H b/amr-wind/equation_systems/tke/source_terms/KransAxell.H index a8e0f0d80a..80c00608ff 100644 --- a/amr-wind/equation_systems/tke/source_terms/KransAxell.H +++ b/amr-wind/equation_systems/tke/source_terms/KransAxell.H @@ -38,19 +38,33 @@ private: amrex::Real m_heat_flux{0.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::Real m_meso_start{600}; 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; - + std::string m_1d_rans; + amrex::Vector m_wind_heights; + amrex::Vector m_tke_values; + amrex::Gpu::DeviceVector m_wind_heights_d; + amrex::Gpu::DeviceVector m_tke_values_d; //! Transport model const transport::TransportModel& m_transport; //! Reference temperature std::unique_ptr m_ref_theta; + bool m_horizontal_sponge{false}; + amrex::Real m_sponge_strength{1.0}; + amrex::Real m_sponge_density{1.0}; + amrex::Real m_sponge_distance_west{-1000}; + amrex::Real m_sponge_distance_east{1000}; + amrex::Real m_sponge_distance_south{-1000}; + amrex::Real m_sponge_distance_north{1000}; + int m_sponge_west{0}; + int m_sponge_east{1}; + int m_sponge_south{0}; + int m_sponge_north{1}; }; } // 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 b763ed2757..a2c5e4e6b2 100644 --- a/amr-wind/equation_systems/tke/source_terms/KransAxell.cpp +++ b/amr-wind/equation_systems/tke/source_terms/KransAxell.cpp @@ -25,12 +25,47 @@ KransAxell::KransAxell(const CFDSim& sim) pp.query("kappa", m_kappa); pp.query("surface_roughness_z0", m_z0); pp.query("surface_temp_flux", m_heat_flux); - pp.query("meso_sponge_start", m_sponge_start); + pp.query("meso_sponge_start", m_meso_start); + pp.query("rans_1dprofile_file", m_1d_rans); + pp.query("horizontal_sponge_tke", m_horizontal_sponge); + 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, value5; + while (ransfile >> value1 >> value2 >> value3 >> value4 >> value5) { + m_wind_heights.push_back(value1); + m_tke_values.push_back(value5); + } + int num_wind_values = static_cast(m_wind_heights.size()); + m_wind_heights_d.resize(num_wind_values); + m_tke_values_d.resize(num_wind_values); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, m_wind_heights.begin(), + m_wind_heights.end(), m_wind_heights_d.begin()); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, m_tke_values.begin(), m_tke_values.end(), + m_tke_values_d.begin()); + } else { + amrex::Abort("Cannot find 1-D RANS profile file " + m_1d_rans); + } { amrex::ParmParse pp_incflow("incflo"); pp_incflow.queryarr("gravity", m_gravity); } m_ref_theta = m_transport.ref_theta(); + amrex::ParmParse pp_drag("DragForcing"); + pp_drag.query("sponge_strength", m_sponge_strength); + pp_drag.query("sponge_density", m_sponge_density); + pp_drag.query("sponge_distance_west", m_sponge_distance_west); + pp_drag.query("sponge_distance_east", m_sponge_distance_east); + pp_drag.query("sponge_distance_south", m_sponge_distance_south); + pp_drag.query("sponge_distance_north", m_sponge_distance_north); + pp_drag.query("sponge_west", m_sponge_west); + pp_drag.query("sponge_east", m_sponge_east); + pp_drag.query("sponge_south", m_sponge_south); + pp_drag.query("sponge_north", m_sponge_north); } KransAxell::~KransAxell() = default; @@ -57,21 +92,25 @@ void KransAxell::operator()( const auto& dt = m_time.delta_t(); 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 bool has_terrain = + this->m_sim.repo().int_field_exists("terrain_blank"); + const amrex::Real sponge_start = m_meso_start; + const auto vsize = m_wind_heights_d.size(); + const auto* wind_heights_d = m_wind_heights_d.data(); + const auto* tke_values_d = m_tke_values_d.data(); 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); - 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 ux = vel(i, j, k + 1, 0); + const amrex::Real uy = vel(i, j, k + 1, 1); const amrex::Real m = std::sqrt(ux * ux + uy * uy); - const amrex::Real ustar = m * kappa / std::log(z / z0); + const amrex::Real ustar = m * kappa / std::log(3 * 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( @@ -79,37 +118,46 @@ void KransAxell::operator()( bcforcing = (ustar * ustar / (Cmu * Cmu) + rans_b - tke_arr(i, j, k)) / dt; } + amrex::Real ref_tke = tke_arr(i, j, k); const amrex::Real zi = std::max((z - sponge_start) / (probhi[2] - sponge_start), 0.0); + if (zi > 0) { + ref_tke = (vsize > 0) ? interp::linear( + wind_heights_d, wind_heights_d + vsize, + tke_values_d, z) + : tke_arr(i, j, k, 0); + } const amrex::Real sponge_forcing = - zi * zi * (tke_arr(i, j, k) - ref_tke); + 1.0 / dt * (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; + src_term(i, j, k) += + shear_prod_arr(i, j, k) + buoy_prod_arr(i, j, k) - + dissip_arr(i, j, k) - + (1 - static_cast(has_terrain)) * (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"); + auto* const m_terrain_height = + &this->m_sim.repo().get_field("terrain_height"); const auto* m_terrain_vf = &this->m_sim.repo().get_field("terrain_vf"); const auto& blank_arr = (*m_terrain_blank)(lev).const_array(mfi); const auto& drag_arr = (*m_terrain_drag)(lev).const_array(mfi); + const auto& terrain_height = (*m_terrain_height)(lev).const_array(mfi); const auto& vf_arr = (*m_terrain_vf)(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 ux = vel(i, j, k + 1, 0); + amrex::Real uy = vel(i, j, k + 1, 1); + 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 ustar = m * kappa / std::log(3.0 * 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( @@ -118,16 +166,84 @@ void KransAxell::operator()( terrainforcing = (ustar * ustar / (Cmu * Cmu) + rans_b - tke_arr(i, j, k)) / dt; + ux = vel(i, j, k, 0); + uy = vel(i, j, k, 1); 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); - + z = std::max( + problo[2] + (k + 0.5) * dx[2] - terrain_height(i, j, k), + 0.5 * dx[2]); + const amrex::Real zi = std::max( + (z - sponge_start) / (probhi[2] - sponge_start), 0.0); + amrex::Real ref_tke = tke_arr(i, j, k); + if (zi > 0) { + ref_tke = (vsize > 0) + ? interp::linear( + wind_heights_d, wind_heights_d + vsize, + tke_values_d, z) + : tke_arr(i, j, k, 0); + } + const amrex::Real sponge_forcing = + 1.0 / dt * (tke_arr(i, j, k) - ref_tke); src_term(i, j, k) += drag_arr(i, j, k) * terrainforcing + - blank_arr(i, j, k) * vf_arr(i, j, k) * dragforcing; + blank_arr(i, j, k) * dragforcing - + static_cast(has_terrain) * sponge_forcing; + ; }); + if (m_horizontal_sponge) { + const amrex::Real sponge_strength = m_sponge_strength; + const amrex::Real sponge_density = m_sponge_density; + const amrex::Real start_east = probhi[0] - m_sponge_distance_east; + const amrex::Real start_west = problo[0] - m_sponge_distance_west; + const amrex::Real start_north = probhi[1] - m_sponge_distance_north; + const amrex::Real start_south = problo[1] - m_sponge_distance_south; + const int sponge_east = m_sponge_east; + const int sponge_west = m_sponge_west; + const int sponge_south = m_sponge_south; + const int sponge_north = m_sponge_north; + amrex::ParallelFor( + bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + const amrex::Real x = problo[0] + (i + 0.5) * dx[0]; + const amrex::Real y = problo[1] + (j + 0.5) * dx[1]; + const amrex::Real z = problo[2] + (k + 0.5) * dx[2]; + amrex::Real xstart_damping = 0; + amrex::Real ystart_damping = 0; + amrex::Real xend_damping = 0; + amrex::Real yend_damping = 0; + amrex::Real xi_end = + (x - start_east) / (probhi[0] - start_east); + amrex::Real xi_start = + (start_west - x) / (start_west - problo[0]); + xi_start = sponge_west * std::max(xi_start, 0.0); + xi_end = sponge_east * std::max(xi_end, 0.0); + xstart_damping = + sponge_west * sponge_strength * xi_start * xi_start; + xend_damping = + sponge_east * sponge_strength * xi_end * xi_end; + amrex::Real yi_end = + (y - start_north) / (probhi[1] - start_north); + amrex::Real yi_start = + (start_south - y) / (start_south - problo[1]); + yi_start = sponge_south * std::max(yi_start, 0.0); + yi_end = sponge_north * std::max(yi_end, 0.0); + ystart_damping = sponge_strength * yi_start * yi_start; + yend_damping = sponge_strength * yi_end * yi_end; + const amrex::Real ref_tke = + (vsize > 0) + ? interp::linear( + wind_heights_d, wind_heights_d + vsize, + tke_values_d, z) + : tke_arr(i, j, k, 0); + src_term(i, j, k, 0) -= + (xstart_damping + xend_damping + ystart_damping + + yend_damping) * + (tke_arr(i, j, k) - sponge_density * ref_tke); + }); + } } } diff --git a/amr-wind/physics/TerrainDrag.cpp b/amr-wind/physics/TerrainDrag.cpp index 6b846774b4..6cd353258f 100644 --- a/amr-wind/physics/TerrainDrag.cpp +++ b/amr-wind/physics/TerrainDrag.cpp @@ -124,8 +124,7 @@ void TerrainDrag::initialize_fields(int level, const amrex::Geometry& geom) static_cast((z <= terrainHt) && (z > prob_lo[2])); levelvf[nbx](i, j, k, 0) = static_cast(levelBlanking[nbx](i, j, k, 0)); - levelheight[nbx](i, j, k, 0) = - std::max(std::abs(z - terrainHt), 0.5 * dx[2]); + levelheight[nbx](i, j, k, 0) = terrainHt; amrex::Real roughz0 = 0.1; if (xrough_size > 0) { diff --git a/amr-wind/turbulence/RANS/KLAxell.H b/amr-wind/turbulence/RANS/KLAxell.H index af19ae1246..bda221f822 100644 --- a/amr-wind/turbulence/RANS/KLAxell.H +++ b/amr-wind/turbulence/RANS/KLAxell.H @@ -59,7 +59,7 @@ private: //! Gravity vector (m/s^2) amrex::Vector m_gravity{0.0, 0.0, -9.81}; amrex::Real m_surf_flux{0}; - amrex::Real m_lengthscale_switch{800}; + amrex::Real m_meso_sponge_start{2000}; }; } // namespace amr_wind::turbulence diff --git a/amr-wind/turbulence/RANS/KLAxell.cpp b/amr-wind/turbulence/RANS/KLAxell.cpp index 44a8f79460..7ab535c57c 100644 --- a/amr-wind/turbulence/RANS/KLAxell.cpp +++ b/amr-wind/turbulence/RANS/KLAxell.cpp @@ -32,7 +32,7 @@ KLAxell::KLAxell(CFDSim& sim) { amrex::ParmParse pp("ABL"); pp.get("surface_temp_flux", m_surf_flux); - pp.query("length_scale_switch", m_lengthscale_switch); + pp.query("meso_sponge_start", m_meso_sponge_start); } { @@ -98,7 +98,7 @@ void KLAxell::update_turbulent_viscosity( 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; + const amrex::Real lengthscale_switch = m_meso_sponge_start; for (int lev = 0; lev < nlevels; ++lev) { const auto& geom = geom_vec[lev]; const auto& problo = repo.mesh().Geom(lev).ProbLoArray();