From 1480a1081fd542021536d0894b410b38db683c32 Mon Sep 17 00:00:00 2001 From: Harish Date: Sun, 5 Jan 2025 13:51:22 -0700 Subject: [PATCH] First push --- .../icns/source_terms/DragForcing.cpp | 12 +++++---- .../source_terms/DragTempForcing.cpp | 5 +++- .../tke/source_terms/KransAxell.cpp | 7 +++-- amr-wind/physics/TerrainDrag.H | 3 +++ amr-wind/physics/TerrainDrag.cpp | 26 +++++++++++++++++++ amr-wind/turbulence/LES/Kosovic.cpp | 9 ++++++- amr-wind/turbulence/RANS/KLAxell.cpp | 13 ++++++---- 7 files changed, 61 insertions(+), 14 deletions(-) diff --git a/amr-wind/equation_systems/icns/source_terms/DragForcing.cpp b/amr-wind/equation_systems/icns/source_terms/DragForcing.cpp index 2a6c6654c1..8939b70519 100644 --- a/amr-wind/equation_systems/icns/source_terms/DragForcing.cpp +++ b/amr-wind/equation_systems/icns/source_terms/DragForcing.cpp @@ -68,6 +68,8 @@ void DragForcing::operator()( const auto& drag = (*m_terrain_drag)(lev).const_array(mfi); auto* const m_terrainz0 = &this->m_sim.repo().get_field("terrainz0"); const auto& terrainz0 = (*m_terrainz0)(lev).const_array(mfi); + const auto* m_terrain_vf = &this->m_sim.repo().get_field("terrain_vf"); + const auto& vf_arr = (*m_terrain_vf)(lev).const_array(mfi); const auto& geom = m_mesh.Geom(lev); const auto& dx = geom.CellSizeArray(); const auto& prob_lo = geom.ProbLoArray(); @@ -161,17 +163,17 @@ void DragForcing::operator()( const amrex::Real CdM = std::min(Cd / (m + tiny), cd_max / scale_factor); src_term(i, j, k, 0) -= - (CdM * m * ux1 * blank(i, j, k) + Dxz * drag(i, j, k) + - bc_forcing_x * drag(i, j, k) + + (CdM * m * ux1 * blank(i, j, k) * vf_arr(i, j, k) + + Dxz * drag(i, j, k) + bc_forcing_x * drag(i, j, k) + (xstart_damping + xend_damping + ystart_damping + yend_damping) * (ux1 - sponge_density * spongeVelX)); src_term(i, j, k, 1) -= - (CdM * m * uy1 * blank(i, j, k) + Dyz * drag(i, j, k) + - bc_forcing_y * drag(i, j, k) + + (CdM * m * uy1 * blank(i, j, k) * vf_arr(i, j, k) + + Dyz * drag(i, j, k) + bc_forcing_y * drag(i, j, k) + (xstart_damping + xend_damping + ystart_damping + yend_damping) * (uy1 - sponge_density * spongeVelY)); src_term(i, j, k, 2) -= - (CdM * m * uz1 * blank(i, j, k) + + (CdM * m * uz1 * blank(i, j, k) * vf_arr(i, j, k) + (xstart_damping + xend_damping + ystart_damping + yend_damping) * (uz1 - sponge_density * spongeVelZ)); }); diff --git a/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.cpp b/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.cpp index f699820d0c..1d839d68cf 100644 --- a/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.cpp +++ b/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.cpp @@ -47,6 +47,8 @@ void DragTempForcing::operator()( auto* const m_terrain_blank = &this->m_sim.repo().get_int_field("terrain_blank"); const auto& blank = (*m_terrain_blank)(lev).const_array(mfi); + const auto* m_terrain_vf = &this->m_sim.repo().get_field("terrain_vf"); + const auto& vf_arr = (*m_terrain_vf)(lev).const_array(mfi); const auto& geom = m_mesh.Geom(lev); const auto& dx = geom.CellSizeArray(); const amrex::Real drag_coefficient = m_drag_coefficient / dx[2]; @@ -62,7 +64,8 @@ void DragTempForcing::operator()( 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) - T0) * blank(i, j, k, 0)); + (Cd * (temperature(i, j, k, 0) - T0) * blank(i, j, k, 0) * + vf_arr(i, j, k, 0)); }); } diff --git a/amr-wind/equation_systems/tke/source_terms/KransAxell.cpp b/amr-wind/equation_systems/tke/source_terms/KransAxell.cpp index df24297bbe..b763ed2757 100644 --- a/amr-wind/equation_systems/tke/source_terms/KransAxell.cpp +++ b/amr-wind/equation_systems/tke/source_terms/KransAxell.cpp @@ -97,8 +97,10 @@ void KransAxell::operator()( &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* 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& 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; @@ -122,8 +124,9 @@ void KransAxell::operator()( 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; + src_term(i, j, k) += + drag_arr(i, j, k) * terrainforcing + + blank_arr(i, j, k) * vf_arr(i, j, k) * dragforcing; }); } } diff --git a/amr-wind/physics/TerrainDrag.H b/amr-wind/physics/TerrainDrag.H index 2009bb0346..1d3bd1a056 100644 --- a/amr-wind/physics/TerrainDrag.H +++ b/amr-wind/physics/TerrainDrag.H @@ -51,6 +51,9 @@ private: //! Roughness fields Field& m_terrainz0; Field& m_terrain_height; + Field& m_terrain_vf; + //! binary, pyramid or EB + std::string m_terrain_cut_model{"pyramid"}; }; } // namespace amr_wind::terraindrag diff --git a/amr-wind/physics/TerrainDrag.cpp b/amr-wind/physics/TerrainDrag.cpp index 24e6d49ad6..71bd0fa203 100644 --- a/amr-wind/physics/TerrainDrag.cpp +++ b/amr-wind/physics/TerrainDrag.cpp @@ -21,20 +21,24 @@ TerrainDrag::TerrainDrag(CFDSim& sim) , m_terrain_drag(sim.repo().declare_int_field("terrain_drag", 1, 1, 1)) , m_terrainz0(sim.repo().declare_field("terrainz0", 1, 1, 1)) , m_terrain_height(sim.repo().declare_field("terrain_height", 1, 1, 1)) + , m_terrain_vf(sim.repo().declare_field("terrain_vf", 1, 1, 1)) { amrex::ParmParse pp(identifier()); pp.query("terrain_file", m_terrain_file); pp.query("roughness_file", m_roughness_file); + pp.query("terrain_cut_model", m_terrain_cut_model); m_sim.io_manager().register_output_int_var("terrain_drag"); m_sim.io_manager().register_output_int_var("terrain_blank"); m_sim.io_manager().register_io_var("terrainz0"); m_sim.io_manager().register_io_var("terrain_height"); + m_sim.io_manager().register_io_var("terrain_vf"); m_terrain_blank.setVal(0.0); m_terrain_drag.setVal(0.0); m_terrainz0.set_default_fillpatch_bc(m_sim.time()); m_terrain_height.set_default_fillpatch_bc(m_sim.time()); + m_terrain_vf.set_default_fillpatch_bc(m_sim.time()); } void TerrainDrag::initialize_fields(int level, const amrex::Geometry& geom) @@ -63,6 +67,7 @@ void TerrainDrag::initialize_fields(int level, const amrex::Geometry& geom) auto& terrainz0 = m_terrainz0(level); auto& terrain_height = m_terrain_height(level); auto& drag = m_terrain_drag(level); + auto& terrain_vf=m_terrain_vf(level); const auto xterrain_size = xterrain.size(); const auto yterrain_size = yterrain.size(); const auto zterrain_size = zterrain.size(); @@ -104,6 +109,8 @@ void TerrainDrag::initialize_fields(int level, const amrex::Geometry& geom) auto levelDrag = drag.arrays(); auto levelz0 = terrainz0.arrays(); auto levelheight = terrain_height.arrays(); + auto levelvf=terrain_vf.arrays(); + std::string terrain_cut_model = m_terrain_cut_model; amrex::ParallelFor( blanking, m_terrain_blank.num_grow(), [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { @@ -115,6 +122,7 @@ void TerrainDrag::initialize_fields(int level, const amrex::Geometry& geom) yterrain_ptr + yterrain_size, zterrain_ptr, x, y); levelBlanking[nbx](i, j, k, 0) = 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]); @@ -138,6 +146,24 @@ void TerrainDrag::initialize_fields(int level, const amrex::Geometry& geom) } }); amrex::Gpu::synchronize(); + if (terrain_cut_model == "pyramid") { + amrex::ParallelFor( + blanking, + [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { + if ((levelDrag[nbx](i, j, k + 1, 0) == 1) && (k > 0)) { + 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]; + const amrex::Real terrainHt = interp::bilinear( + xterrain_ptr, xterrain_ptr + xterrain_size, + yterrain_ptr, yterrain_ptr + yterrain_size, + zterrain_ptr, x, y); + const amrex::Real dz = terrainHt - (z - 0.5 * dx[2]); + levelvf[nbx](i, j, k, 0) = 0.33 * dz / dx[2]; + } + }); + amrex::Gpu::synchronize(); + } } void TerrainDrag::post_regrid_actions() diff --git a/amr-wind/turbulence/LES/Kosovic.cpp b/amr-wind/turbulence/LES/Kosovic.cpp index dfae6d5f0a..66aa820593 100644 --- a/amr-wind/turbulence/LES/Kosovic.cpp +++ b/amr-wind/turbulence/LES/Kosovic.cpp @@ -59,6 +59,8 @@ void Kosovic::update_turbulent_viscosity( const auto* m_terrain_blank = has_terrain ? &this->m_sim.repo().get_int_field("terrain_blank") : nullptr; + const auto* m_terrain_vf = + has_terrain ? &this->m_sim.repo().get_field("terrain_vf") : nullptr; // Populate strainrate into the turbulent viscosity arrays to avoid creating // a temporary buffer fvm::strainrate(mu_turb, vel); @@ -89,6 +91,9 @@ void Kosovic::update_turbulent_viscosity( const auto& blank_arr = has_terrain ? (*m_terrain_blank)(lev).const_array(mfi) : amrex::Array4(); + const auto& vf_arr = has_terrain + ? (*m_terrain_vf)(lev).const_array(mfi) + : amrex::Array4(); amrex::ParallelFor( bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { const amrex::Real rho = rho_arr(i, j, k); @@ -107,7 +112,9 @@ void Kosovic::update_turbulent_viscosity( std::pow(fmu, locSurfaceRANSExp) * ransL) + (1 - locSurfaceFactor) * smag_factor; const amrex::Real blankTerrain = - (has_terrain) ? 1 - blank_arr(i, j, k, 0) : 1.0; + (has_terrain) + ? 1 - blank_arr(i, j, k, 0) * vf_arr(i, j, k, 0) + : 1.0; mu_arr(i, j, k) *= rho * viscosityScale * turnOff * blankTerrain; amrex::Real stressScale = diff --git a/amr-wind/turbulence/RANS/KLAxell.cpp b/amr-wind/turbulence/RANS/KLAxell.cpp index 1d2af84972..44a8f79460 100644 --- a/amr-wind/turbulence/RANS/KLAxell.cpp +++ b/amr-wind/turbulence/RANS/KLAxell.cpp @@ -122,9 +122,12 @@ void KLAxell::update_turbulent_viscosity( &this->m_sim.repo().get_field("terrain_height"); const auto* m_terrain_blank = &this->m_sim.repo().get_int_field("terrain_blank"); + const auto* m_terrain_vf = + &this->m_sim.repo().get_field("terrain_vf"); const auto& ht_arr = (*m_terrain_height)(lev).const_array(mfi); const auto& blank_arr = (*m_terrain_blank)(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 stratification = @@ -182,15 +185,15 @@ void KLAxell::update_turbulent_viscosity( 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)); + 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) * vf_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)); + (1 - blank_arr(i, j, k) * vf_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);