From 27b0cb78df9182872b77b84624833608a184fd46 Mon Sep 17 00:00:00 2001 From: Harish Date: Mon, 23 Jun 2025 05:04:50 -0400 Subject: [PATCH 01/10] First --- .../icns/source_terms/DragForcing.cpp | 79 +++++++++++++++++-- .../source_terms/DragTempForcing.cpp | 50 +++++++++++- .../tke/source_terms/KransAxell.cpp | 7 +- amr-wind/physics/TerrainDrag.H | 2 + amr-wind/physics/TerrainDrag.cpp | 65 +++++++++++++-- 5 files changed, 183 insertions(+), 20 deletions(-) diff --git a/amr-wind/equation_systems/icns/source_terms/DragForcing.cpp b/amr-wind/equation_systems/icns/source_terms/DragForcing.cpp index 80ecaa7f88..f3e03f1998 100644 --- a/amr-wind/equation_systems/icns/source_terms/DragForcing.cpp +++ b/amr-wind/equation_systems/icns/source_terms/DragForcing.cpp @@ -10,6 +10,27 @@ #include "amr-wind/utilities/constants.H" namespace { + +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::GpuArray +compute_target_wind( + const amrex::Real ux, + const amrex::Real uy, + const amrex::Real dx, + const amrex::Real z0, + const amrex::Real kappa) +{ + const amrex::Real wspd = std::sqrt(ux * ux + uy * uy); + const amrex::Real ustar = wspd * kappa / std::log(1.5 * dx / z0); + const amrex::Real wspd_target = ustar / kappa * std::log(0.5 * dx / z0); + const amrex::Real wspd_target_x = + wspd_target * ux / + (amr_wind::constants::EPS + std::sqrt(ux * ux + uy * uy)); + const amrex::Real wspd_target_y = + wspd_target * uy / + (amr_wind::constants::EPS + std::sqrt(ux * ux + uy * uy)); + return {wspd_target_x, wspd_target_y}; +} + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real viscous_drag_calculations( amrex::Real& Dxz, amrex::Real& Dyz, @@ -209,6 +230,7 @@ void DragForcing::operator()( 0.5 * dx[2] / m_monin_obukhov_length, m_beta_m, m_gamma_m) : 0.0; amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + const int cell_drag = (drag(i, j, k) > 0) ? 1 : 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]; @@ -249,7 +271,9 @@ void DragForcing::operator()( amrex::Real Dyz = 0.0; amrex::Real bc_forcing_x = 0; amrex::Real bc_forcing_y = 0; + amrex::Real bc_forcing_z = 0; const amrex::Real m = std::sqrt(ux1 * ux1 + uy1 * uy1 + uz1 * uz1); + const amrex::Real z0 = std::max(terrainz0(i, j, k), z0_min); if (drag(i, j, k) == 1 && (!is_laminar)) { // Check if close enough to interface to use current cell or below int k_off = -1; @@ -273,8 +297,7 @@ void DragForcing::operator()( const amrex::Real uy1r = uy1 - wall_v; const amrex::Real ux2r = vel(i, j, k + 1, 0) - wall_u; const amrex::Real uy2r = vel(i, j, k + 1, 1) - wall_v; - const amrex::Real z0 = std::max(terrainz0(i, j, k), z0_min); - const amrex::Real ustar = viscous_drag_calculations( + amrex::Real ustar = viscous_drag_calculations( Dxz, Dyz, ux1r, uy1r, ux2r, uy2r, z0, dx[2], kappa, non_neutral_neighbour); if (model_form_drag) { @@ -292,6 +315,48 @@ void DragForcing::operator()( // BC forcing pushes nonrelative velocity toward target velocity bc_forcing_x = -(uxTarget - ux1) / dt; bc_forcing_y = -(uyTarget - uy1) / dt; + //! Adding horizonal drag + if(drag(i,j,k)>1){ + //! West + amrex::GpuArray tmp_wind_target = + compute_target_wind( + vel(i - 1, j, k, 2), vel(i - 1, j, k, 1), dx[0], z0, kappa); + bc_forcing_z += + -(tmp_wind_target[0] - uz1) / dt * blank(i - 1, j, k); + bc_forcing_y += + -(tmp_wind_target[1] - uy1) / dt * blank(i - 1, j, k); + //! East + tmp_wind_target = compute_target_wind( + vel(i + 1, j, k, 2), vel(i + 1, j, k, 1), dx[0], z0, kappa); + bc_forcing_z += + -(tmp_wind_target[0] - uz1) / dt * blank(i + 1, j, k); + bc_forcing_y += + -(tmp_wind_target[1] - uy1) / dt * blank(i + 1, j, k); + //! South + tmp_wind_target = compute_target_wind( + vel(i, j - 1, k, 2), vel(i, j - 1, k, 0), dx[1], z0, kappa); + bc_forcing_z += + -(tmp_wind_target[0] - uz1) / dt * blank(i, j - 1, k); + bc_forcing_x += + -(tmp_wind_target[1] - ux1) / dt * blank(i, j - 1, k); + //! North + tmp_wind_target = compute_target_wind( + vel(i, j + 1, k, 2), vel(i, j + 1, k, 0), dx[1], z0, kappa); + bc_forcing_z += + -(tmp_wind_target[0] - uz1) / dt * blank(i, j + 1, k); + bc_forcing_x += + -(tmp_wind_target[1] - ux1) / dt * blank(i, j + 1, k); + const amrex::Real sum_blank_x = + blank(i, j - 1, k) + blank(i, j + 1, k) + blank(i, j, k - 1); + bc_forcing_x /= (sum_blank_x + amr_wind::constants::EPS); + const amrex::Real sum_blank_y = + blank(i - 1, j, k) + blank(i + 1, j, k) + blank(i, j, k - 1); + bc_forcing_y /= (sum_blank_y + amr_wind::constants::EPS); + const amrex::Real sum_blank_z = + blank(i - 1, j, k) + blank(i + 1, j, k) + blank(i, j - 1, k) + + blank(i, j + 1, k); + bc_forcing_z /= (sum_blank_z + amr_wind::constants::EPS); + } } // Target velocity intended for within terrain amrex::Real target_u = 0.; @@ -302,21 +367,21 @@ void DragForcing::operator()( target_v = target_vel_arr(i, j, k, 1); target_w = target_vel_arr(i, j, k, 2); } - const amrex::Real CdM = std::min( Cd / (m + amr_wind::constants::EPS), cd_max / scale_factor); src_term(i, j, k, 0) -= - (CdM * m * (ux1 - target_u) * blank(i, j, k) + Dxz * drag(i, j, k) + - bc_forcing_x * drag(i, j, k) + + (CdM * m * (ux1 - target_u) * blank(i, j, k) + Dxz * cell_drag + + bc_forcing_x * cell_drag + (xstart_damping + xend_damping + ystart_damping + yend_damping) * (ux1 - sponge_density * spongeVelX)); src_term(i, j, k, 1) -= - (CdM * m * (uy1 - target_v) * blank(i, j, k) + Dyz * drag(i, j, k) + - bc_forcing_y * drag(i, j, k) + + (CdM * m * (uy1 - target_v) * blank(i, j, k) + Dyz * cell_drag + + bc_forcing_y * cell_drag + (xstart_damping + xend_damping + ystart_damping + yend_damping) * (uy1 - sponge_density * spongeVelY)); src_term(i, j, k, 2) -= (CdM * m * (uz1 - target_w) * blank(i, j, k) + + bc_forcing_z * cell_drag + (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 b431f8e8f9..5f9425db1c 100644 --- a/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.cpp +++ b/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.cpp @@ -6,7 +6,34 @@ #include "AMReX_ParmParse.H" #include "AMReX_Gpu.H" #include "AMReX_Random.H" +#include "amr-wind/utilities/constants.H" +namespace{ +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real +compute_target_theta( + const amrex::Real ux, + const amrex::Real uy, + const amrex::Real theta, + const amrex::Real monin_obukhov_length, + const amrex::Real gravity_mod, + const amrex::Real dx, + const amrex::Real z0, + const amrex::Real kappa) +{ + const amrex::Real wspd = std::sqrt(ux * ux + uy * uy); + const amrex::Real ustar = wspd * kappa / std::log(1.5 * dx / z0); + const amrex::Real thetastar = + theta * ustar * ustar / + (kappa * gravity_mod * monin_obukhov_length); + const amrex::Real surf_temp = + theta - + thetastar / kappa * (std::log(1.5 * dx / z0)); + const amrex::Real tTarget = + surf_temp + + thetastar / kappa * (std::log(0.5 * dx / z0)); + return tTarget; +} +} namespace amr_wind::pde::temperature { DragTempForcing::DragTempForcing(const CFDSim& sim) @@ -87,6 +114,7 @@ void DragTempForcing::operator()( const amrex::Real cd_max = 10.0; const amrex::Real T0 = m_soil_temperature; amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + const int cell_drag = (drag(i,j,k)> 0) ? 1 : 0; const amrex::Real z0 = std::max(terrainz0(i, j, k), z0_min); const amrex::Real ux1 = vel(i, j, k, 0); const amrex::Real uy1 = vel(i, j, k, 1); @@ -106,13 +134,31 @@ void DragTempForcing::operator()( const amrex::Real tTarget = surf_temp + thetastar / kappa * (std::log(0.5 * dx[2] / z0) - psi_h_cell); - const amrex::Real bc_forcing_t = -(tTarget - theta) / dt; + amrex::Real bc_forcing_t = -(tTarget - theta) / dt; + // //! West + // amrex::Real tmp_temp_target =compute_target_theta(vel(i-1,j,k,2),vel(i-1,j,k,1),theta, + // monin_obukhov_length, gravity_mod, dx[0],z0, kappa); + // bc_forcing_t += -(tmp_temp_target - theta) / dt * blank(i-1,j,k); + // //! East + // tmp_temp_target = compute_target_theta(vel(i+1,j,k,2),vel(i+1,j,k,1),theta, + // monin_obukhov_length, gravity_mod, dx[0],z0, kappa); + // bc_forcing_t += -(tmp_temp_target - theta) / dt * blank(i+1,j,k); + // //! South + // tmp_temp_target = compute_target_theta(vel(i,j-1,k,2),vel(i,j-1,k,0),theta, + // monin_obukhov_length, gravity_mod, dx[1],z0, kappa); + // bc_forcing_t += -(tmp_temp_target - theta) / dt * blank(i,j-1,k); + // //! North + // tmp_temp_target = compute_target_theta(vel(i,j+1,k,2),vel(i,j+1,k,0),theta, + // monin_obukhov_length, gravity_mod, dx[1],z0, kappa); + // bc_forcing_t += -(tmp_temp_target - theta) / dt * blank(i,j+1,k); + // const amrex::Real sum_blank_t = blank(i,j,k-1)+blank(i - 1, j, k) + blank(i + 1, j, k) + blank(i, j - 1, k) + blank(i, j + 1, k); + // bc_forcing_t /= (sum_blank_t + amr_wind::constants::EPS); 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]); src_term(i, j, k, 0) -= (Cd * (theta - T0) * blank(i, j, k, 0) + - bc_forcing_t * drag(i, j, k)); + bc_forcing_t * cell_drag); }); } diff --git a/amr-wind/equation_systems/tke/source_terms/KransAxell.cpp b/amr-wind/equation_systems/tke/source_terms/KransAxell.cpp index e81b9601a6..32a299d126 100644 --- a/amr-wind/equation_systems/tke/source_terms/KransAxell.cpp +++ b/amr-wind/equation_systems/tke/source_terms/KransAxell.cpp @@ -164,9 +164,10 @@ void KransAxell::operator()( const auto& terrainz0 = (*m_terrainz0)(lev).const_array(mfi); amrex::ParallelFor( bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + const int cell_drag = (drag_arr(i, j, k) > 0) ? 1 : 0; const amrex::Real cell_z0 = - drag_arr(i, j, k) * std::max(terrainz0(i, j, k), z0_min) + - (1 - drag_arr(i, j, k)) * z0; + cell_drag * std::max(terrainz0(i, j, k), z0_min) + + (1 - cell_drag) * z0; amrex::Real terrainforcing = 0; amrex::Real dragforcing = 0; amrex::Real ux = vel(i, j, k + 1, 0); @@ -211,7 +212,7 @@ void KransAxell::operator()( 1.0 / dt * (tke_arr(i, j, k) - ref_tke); src_term(i, j, k) = (1 - blank_arr(i, j, k)) * src_term(i, j, k) + - drag_arr(i, j, k) * terrainforcing + + cell_drag * terrainforcing + blank_arr(i, j, k) * dragforcing - static_cast(has_terrain) * (sponge_forcing - bcforcing); diff --git a/amr-wind/physics/TerrainDrag.H b/amr-wind/physics/TerrainDrag.H index 3495eb88eb..86eeb8d411 100644 --- a/amr-wind/physics/TerrainDrag.H +++ b/amr-wind/physics/TerrainDrag.H @@ -5,6 +5,7 @@ #include "amr-wind/core/Field.H" #include "amr-wind/CFDSim.H" + namespace amr_wind::terraindrag { /** Terraindrag Flow physics @@ -56,6 +57,7 @@ private: //! Terrain drag force term IntField& m_terrain_drag; + //! Terrain file std::string m_terrain_file{"terrain.amrwind"}; diff --git a/amr-wind/physics/TerrainDrag.cpp b/amr-wind/physics/TerrainDrag.cpp index dd8096aa15..376695f79d 100644 --- a/amr-wind/physics/TerrainDrag.cpp +++ b/amr-wind/physics/TerrainDrag.cpp @@ -8,7 +8,8 @@ #include "amr-wind/utilities/IOManager.H" #include "amr-wind/utilities/io_utils.H" #include "amr-wind/utilities/linear_interpolation.H" - +#include "amr-wind/fvm/gradient.H" +#include "amr-wind/utilities/constants.H" namespace amr_wind::terraindrag { namespace {} // namespace @@ -140,16 +141,64 @@ void TerrainDrag::initialize_fields(int level, const amrex::Geometry& geom) levelz0[nbx](i, j, k, 0) = roughz0; }); amrex::Gpu::streamSynchronize(); + // amrex::ParallelFor( + // blanking, [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { + // const amrex::GpuArray cell_blanking{ + // levelBlanking[nbx](i, j, k - 1, 0), + // levelBlanking[nbx](i - 1, j, k, 0), + // levelBlanking[nbx](i + 1, j, k, 0), + // levelBlanking[nbx](i, j - 1, k, 0), + // levelBlanking[nbx](i, j + 1, k, 0)}; + // levelDrag[nbx](i, j, k, 0) = 0; + // for (int index = 0; index <= 4; ++index) { + // if ((levelBlanking[nbx](i, j, k, 0) == 0) && + // (cell_blanking[index] == 1) && (levelDrag[nbx](i, j, k, 0) == 0)) { + // levelDrag[nbx](i, j, k, 0) = index + 1; + // } + // } + // }); + //amrex::Gpu::streamSynchronize(); + + // amrex::ParallelFor( + // blanking, + // [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { + // if ((levelBlanking[nbx](i, j, k, 0) == 0) && + // (levelBlanking[nbx](i, j, k - 1, 0) == 1)) { + // levelDrag[nbx](i, j, k, 0) = 1; + // } else { + // levelDrag[nbx](i, j, k, 0) = 0; + // } + // }); + // amrex::Gpu::streamSynchronize(); + // amrex::ParallelFor( + // blanking, [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { + // const amrex::GpuArray cell_blanking{ + // levelBlanking[nbx](i - 1, j, k, 0), + // levelBlanking[nbx](i + 1, j, k, 0), + // levelBlanking[nbx](i, j - 1, k, 0), + // levelBlanking[nbx](i, j + 1, k, 0)}; + // levelDrag[nbx](i, j, k, 0) = 0; + // for (int index = 0; index <= 3; ++index) { + // if ((levelBlanking[nbx](i, j, k, 0) == 0) && + // (cell_blanking[index] == 1) ) { + // levelDrag[nbx](i, j, k, 0) = 1; + // } + // } + // }); + // amrex::Gpu::streamSynchronize(); amrex::ParallelFor( blanking, [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { - if ((levelBlanking[nbx](i, j, k, 0) == 0) && (k > 0) && - (levelBlanking[nbx](i, j, k - 1, 0) == 1)) { - levelDrag[nbx](i, j, k, 0) = 1; - } else { - levelDrag[nbx](i, j, k, 0) = 0; - } + const int blankxp= std::abs(levelBlanking[nbx](i + 1, j, k, 0) - levelBlanking[nbx](i, j, k, 0)); + const int blankxm= std::abs(levelBlanking[nbx](i - 1, j, k, 0) - levelBlanking[nbx](i, j, k, 0)); + const int blankyp= std::abs(levelBlanking[nbx](i, j + 1, k, 0) - levelBlanking[nbx](i, j, k, 0)); + const int blankym= std::abs(levelBlanking[nbx](i, j - 1, k, 0) - levelBlanking[nbx](i, j, k, 0)); + const int blankzp= std::abs(levelBlanking[nbx](i, j, k + 1, 0) - levelBlanking[nbx](i, j, k, 0)); + const int blankzm= std::abs(levelBlanking[nbx](i, j, k - 1, 0) - levelBlanking[nbx](i, j, k, 0)); + levelDrag[nbx](i,j,k,0)= ((k==0 && levelBlanking[nbx](i, j, k, 0) == 1))? + 0: + (blankxp + blankxm + blankyp + blankym + blankzp + blankzm); }); - amrex::Gpu::streamSynchronize(); + } void TerrainDrag::post_init_actions() From 85eff5a1fa7ad0db82b84f0f6046e113ebde989f Mon Sep 17 00:00:00 2001 From: Harish Date: Mon, 23 Jun 2025 05:10:11 -0400 Subject: [PATCH 02/10] Second --- .../icns/source_terms/DragForcing.cpp | 83 ++++++++++--------- .../source_terms/DragTempForcing.cpp | 36 ++++---- amr-wind/physics/TerrainDrag.H | 2 - amr-wind/physics/TerrainDrag.cpp | 47 +++++++---- 4 files changed, 92 insertions(+), 76 deletions(-) diff --git a/amr-wind/equation_systems/icns/source_terms/DragForcing.cpp b/amr-wind/equation_systems/icns/source_terms/DragForcing.cpp index f3e03f1998..f0384058ff 100644 --- a/amr-wind/equation_systems/icns/source_terms/DragForcing.cpp +++ b/amr-wind/equation_systems/icns/source_terms/DragForcing.cpp @@ -316,46 +316,49 @@ void DragForcing::operator()( bc_forcing_x = -(uxTarget - ux1) / dt; bc_forcing_y = -(uyTarget - uy1) / dt; //! Adding horizonal drag - if(drag(i,j,k)>1){ - //! West - amrex::GpuArray tmp_wind_target = - compute_target_wind( - vel(i - 1, j, k, 2), vel(i - 1, j, k, 1), dx[0], z0, kappa); - bc_forcing_z += - -(tmp_wind_target[0] - uz1) / dt * blank(i - 1, j, k); - bc_forcing_y += - -(tmp_wind_target[1] - uy1) / dt * blank(i - 1, j, k); - //! East - tmp_wind_target = compute_target_wind( - vel(i + 1, j, k, 2), vel(i + 1, j, k, 1), dx[0], z0, kappa); - bc_forcing_z += - -(tmp_wind_target[0] - uz1) / dt * blank(i + 1, j, k); - bc_forcing_y += - -(tmp_wind_target[1] - uy1) / dt * blank(i + 1, j, k); - //! South - tmp_wind_target = compute_target_wind( - vel(i, j - 1, k, 2), vel(i, j - 1, k, 0), dx[1], z0, kappa); - bc_forcing_z += - -(tmp_wind_target[0] - uz1) / dt * blank(i, j - 1, k); - bc_forcing_x += - -(tmp_wind_target[1] - ux1) / dt * blank(i, j - 1, k); - //! North - tmp_wind_target = compute_target_wind( - vel(i, j + 1, k, 2), vel(i, j + 1, k, 0), dx[1], z0, kappa); - bc_forcing_z += - -(tmp_wind_target[0] - uz1) / dt * blank(i, j + 1, k); - bc_forcing_x += - -(tmp_wind_target[1] - ux1) / dt * blank(i, j + 1, k); - const amrex::Real sum_blank_x = - blank(i, j - 1, k) + blank(i, j + 1, k) + blank(i, j, k - 1); - bc_forcing_x /= (sum_blank_x + amr_wind::constants::EPS); - const amrex::Real sum_blank_y = - blank(i - 1, j, k) + blank(i + 1, j, k) + blank(i, j, k - 1); - bc_forcing_y /= (sum_blank_y + amr_wind::constants::EPS); - const amrex::Real sum_blank_z = - blank(i - 1, j, k) + blank(i + 1, j, k) + blank(i, j - 1, k) + - blank(i, j + 1, k); - bc_forcing_z /= (sum_blank_z + amr_wind::constants::EPS); + if (drag(i, j, k) > 1) { + //! West + amrex::GpuArray tmp_wind_target = + compute_target_wind( + vel(i - 1, j, k, 2), vel(i - 1, j, k, 1), dx[0], z0, + kappa); + bc_forcing_z += + -(tmp_wind_target[0] - uz1) / dt * blank(i - 1, j, k); + bc_forcing_y += + -(tmp_wind_target[1] - uy1) / dt * blank(i - 1, j, k); + //! East + tmp_wind_target = compute_target_wind( + vel(i + 1, j, k, 2), vel(i + 1, j, k, 1), dx[0], z0, kappa); + bc_forcing_z += + -(tmp_wind_target[0] - uz1) / dt * blank(i + 1, j, k); + bc_forcing_y += + -(tmp_wind_target[1] - uy1) / dt * blank(i + 1, j, k); + //! South + tmp_wind_target = compute_target_wind( + vel(i, j - 1, k, 2), vel(i, j - 1, k, 0), dx[1], z0, kappa); + bc_forcing_z += + -(tmp_wind_target[0] - uz1) / dt * blank(i, j - 1, k); + bc_forcing_x += + -(tmp_wind_target[1] - ux1) / dt * blank(i, j - 1, k); + //! North + tmp_wind_target = compute_target_wind( + vel(i, j + 1, k, 2), vel(i, j + 1, k, 0), dx[1], z0, kappa); + bc_forcing_z += + -(tmp_wind_target[0] - uz1) / dt * blank(i, j + 1, k); + bc_forcing_x += + -(tmp_wind_target[1] - ux1) / dt * blank(i, j + 1, k); + const amrex::Real sum_blank_x = blank(i, j - 1, k) + + blank(i, j + 1, k) + + blank(i, j, k - 1); + bc_forcing_x /= (sum_blank_x + amr_wind::constants::EPS); + const amrex::Real sum_blank_y = blank(i - 1, j, k) + + blank(i + 1, j, k) + + blank(i, j, k - 1); + bc_forcing_y /= (sum_blank_y + amr_wind::constants::EPS); + const amrex::Real sum_blank_z = + blank(i - 1, j, k) + blank(i + 1, j, k) + + blank(i, j - 1, k) + blank(i, j + 1, k); + bc_forcing_z /= (sum_blank_z + amr_wind::constants::EPS); } } // Target velocity intended for within terrain diff --git a/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.cpp b/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.cpp index 5f9425db1c..ef8e14afaa 100644 --- a/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.cpp +++ b/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.cpp @@ -8,9 +8,8 @@ #include "AMReX_Random.H" #include "amr-wind/utilities/constants.H" -namespace{ -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real -compute_target_theta( +namespace { +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real compute_target_theta( const amrex::Real ux, const amrex::Real uy, const amrex::Real theta, @@ -23,17 +22,14 @@ compute_target_theta( const amrex::Real wspd = std::sqrt(ux * ux + uy * uy); const amrex::Real ustar = wspd * kappa / std::log(1.5 * dx / z0); const amrex::Real thetastar = - theta * ustar * ustar / - (kappa * gravity_mod * monin_obukhov_length); + theta * ustar * ustar / (kappa * gravity_mod * monin_obukhov_length); const amrex::Real surf_temp = - theta - - thetastar / kappa * (std::log(1.5 * dx / z0)); + theta - thetastar / kappa * (std::log(1.5 * dx / z0)); const amrex::Real tTarget = - surf_temp + - thetastar / kappa * (std::log(0.5 * dx / z0)); + surf_temp + thetastar / kappa * (std::log(0.5 * dx / z0)); return tTarget; } -} +} // namespace namespace amr_wind::pde::temperature { DragTempForcing::DragTempForcing(const CFDSim& sim) @@ -114,7 +110,7 @@ void DragTempForcing::operator()( const amrex::Real cd_max = 10.0; const amrex::Real T0 = m_soil_temperature; amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - const int cell_drag = (drag(i,j,k)> 0) ? 1 : 0; + const int cell_drag = (drag(i, j, k) > 0) ? 1 : 0; const amrex::Real z0 = std::max(terrainz0(i, j, k), z0_min); const amrex::Real ux1 = vel(i, j, k, 0); const amrex::Real uy1 = vel(i, j, k, 1); @@ -136,29 +132,33 @@ void DragTempForcing::operator()( thetastar / kappa * (std::log(0.5 * dx[2] / z0) - psi_h_cell); amrex::Real bc_forcing_t = -(tTarget - theta) / dt; // //! West - // amrex::Real tmp_temp_target =compute_target_theta(vel(i-1,j,k,2),vel(i-1,j,k,1),theta, + // amrex::Real tmp_temp_target + // =compute_target_theta(vel(i-1,j,k,2),vel(i-1,j,k,1),theta, // monin_obukhov_length, gravity_mod, dx[0],z0, kappa); // bc_forcing_t += -(tmp_temp_target - theta) / dt * blank(i-1,j,k); // //! East - // tmp_temp_target = compute_target_theta(vel(i+1,j,k,2),vel(i+1,j,k,1),theta, + // tmp_temp_target = + // compute_target_theta(vel(i+1,j,k,2),vel(i+1,j,k,1),theta, // monin_obukhov_length, gravity_mod, dx[0],z0, kappa); // bc_forcing_t += -(tmp_temp_target - theta) / dt * blank(i+1,j,k); // //! South - // tmp_temp_target = compute_target_theta(vel(i,j-1,k,2),vel(i,j-1,k,0),theta, + // tmp_temp_target = + // compute_target_theta(vel(i,j-1,k,2),vel(i,j-1,k,0),theta, // monin_obukhov_length, gravity_mod, dx[1],z0, kappa); // bc_forcing_t += -(tmp_temp_target - theta) / dt * blank(i,j-1,k); // //! North - // tmp_temp_target = compute_target_theta(vel(i,j+1,k,2),vel(i,j+1,k,0),theta, + // tmp_temp_target = + // compute_target_theta(vel(i,j+1,k,2),vel(i,j+1,k,0),theta, // monin_obukhov_length, gravity_mod, dx[1],z0, kappa); // bc_forcing_t += -(tmp_temp_target - theta) / dt * blank(i,j+1,k); - // const amrex::Real sum_blank_t = blank(i,j,k-1)+blank(i - 1, j, k) + blank(i + 1, j, k) + blank(i, j - 1, k) + blank(i, j + 1, k); + // const amrex::Real sum_blank_t = blank(i,j,k-1)+blank(i - 1, j, k) + + // blank(i + 1, j, k) + blank(i, j - 1, k) + blank(i, j + 1, k); // bc_forcing_t /= (sum_blank_t + amr_wind::constants::EPS); 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]); src_term(i, j, k, 0) -= - (Cd * (theta - T0) * blank(i, j, k, 0) + - bc_forcing_t * cell_drag); + (Cd * (theta - T0) * blank(i, j, k, 0) + bc_forcing_t * cell_drag); }); } diff --git a/amr-wind/physics/TerrainDrag.H b/amr-wind/physics/TerrainDrag.H index 86eeb8d411..3495eb88eb 100644 --- a/amr-wind/physics/TerrainDrag.H +++ b/amr-wind/physics/TerrainDrag.H @@ -5,7 +5,6 @@ #include "amr-wind/core/Field.H" #include "amr-wind/CFDSim.H" - namespace amr_wind::terraindrag { /** Terraindrag Flow physics @@ -57,7 +56,6 @@ private: //! Terrain drag force term IntField& m_terrain_drag; - //! Terrain file std::string m_terrain_file{"terrain.amrwind"}; diff --git a/amr-wind/physics/TerrainDrag.cpp b/amr-wind/physics/TerrainDrag.cpp index 376695f79d..9c883028d6 100644 --- a/amr-wind/physics/TerrainDrag.cpp +++ b/amr-wind/physics/TerrainDrag.cpp @@ -142,7 +142,8 @@ void TerrainDrag::initialize_fields(int level, const amrex::Geometry& geom) }); amrex::Gpu::streamSynchronize(); // amrex::ParallelFor( - // blanking, [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { + // blanking, [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept + // { // const amrex::GpuArray cell_blanking{ // levelBlanking[nbx](i, j, k - 1, 0), // levelBlanking[nbx](i - 1, j, k, 0), @@ -152,17 +153,17 @@ void TerrainDrag::initialize_fields(int level, const amrex::Geometry& geom) // levelDrag[nbx](i, j, k, 0) = 0; // for (int index = 0; index <= 4; ++index) { // if ((levelBlanking[nbx](i, j, k, 0) == 0) && - // (cell_blanking[index] == 1) && (levelDrag[nbx](i, j, k, 0) == 0)) { - // levelDrag[nbx](i, j, k, 0) = index + 1; + // (cell_blanking[index] == 1) && (levelDrag[nbx](i, j, k, + // 0) == 0)) { levelDrag[nbx](i, j, k, 0) = index + 1; // } // } // }); - //amrex::Gpu::streamSynchronize(); + // amrex::Gpu::streamSynchronize(); // amrex::ParallelFor( // blanking, // [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { - // if ((levelBlanking[nbx](i, j, k, 0) == 0) && + // if ((levelBlanking[nbx](i, j, k, 0) == 0) && // (levelBlanking[nbx](i, j, k - 1, 0) == 1)) { // levelDrag[nbx](i, j, k, 0) = 1; // } else { @@ -171,7 +172,8 @@ void TerrainDrag::initialize_fields(int level, const amrex::Geometry& geom) // }); // amrex::Gpu::streamSynchronize(); // amrex::ParallelFor( - // blanking, [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { + // blanking, [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept + // { // const amrex::GpuArray cell_blanking{ // levelBlanking[nbx](i - 1, j, k, 0), // levelBlanking[nbx](i + 1, j, k, 0), @@ -188,17 +190,30 @@ void TerrainDrag::initialize_fields(int level, const amrex::Geometry& geom) // amrex::Gpu::streamSynchronize(); amrex::ParallelFor( blanking, [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { - const int blankxp= std::abs(levelBlanking[nbx](i + 1, j, k, 0) - levelBlanking[nbx](i, j, k, 0)); - const int blankxm= std::abs(levelBlanking[nbx](i - 1, j, k, 0) - levelBlanking[nbx](i, j, k, 0)); - const int blankyp= std::abs(levelBlanking[nbx](i, j + 1, k, 0) - levelBlanking[nbx](i, j, k, 0)); - const int blankym= std::abs(levelBlanking[nbx](i, j - 1, k, 0) - levelBlanking[nbx](i, j, k, 0)); - const int blankzp= std::abs(levelBlanking[nbx](i, j, k + 1, 0) - levelBlanking[nbx](i, j, k, 0)); - const int blankzm= std::abs(levelBlanking[nbx](i, j, k - 1, 0) - levelBlanking[nbx](i, j, k, 0)); - levelDrag[nbx](i,j,k,0)= ((k==0 && levelBlanking[nbx](i, j, k, 0) == 1))? - 0: - (blankxp + blankxm + blankyp + blankym + blankzp + blankzm); + const int blankxp = std::abs( + levelBlanking[nbx](i + 1, j, k, 0) - + levelBlanking[nbx](i, j, k, 0)); + const int blankxm = std::abs( + levelBlanking[nbx](i - 1, j, k, 0) - + levelBlanking[nbx](i, j, k, 0)); + const int blankyp = std::abs( + levelBlanking[nbx](i, j + 1, k, 0) - + levelBlanking[nbx](i, j, k, 0)); + const int blankym = std::abs( + levelBlanking[nbx](i, j - 1, k, 0) - + levelBlanking[nbx](i, j, k, 0)); + const int blankzp = std::abs( + levelBlanking[nbx](i, j, k + 1, 0) - + levelBlanking[nbx](i, j, k, 0)); + const int blankzm = std::abs( + levelBlanking[nbx](i, j, k - 1, 0) - + levelBlanking[nbx](i, j, k, 0)); + levelDrag[nbx](i, j, k, 0) = + ((k == 0 && levelBlanking[nbx](i, j, k, 0) == 1)) + ? 0 + : (blankxp + blankxm + blankyp + blankym + blankzp + + blankzm); }); - } void TerrainDrag::post_init_actions() From 106583bbef7418d88649925e551c419ddbfe42d7 Mon Sep 17 00:00:00 2001 From: Harish Date: Mon, 23 Jun 2025 05:58:03 -0400 Subject: [PATCH 03/10] Third --- .../source_terms/DragTempForcing.cpp | 53 +++++++++++-------- .../tke/source_terms/KransAxell.cpp | 53 ++++++++++++++++++- amr-wind/physics/TerrainDrag.cpp | 48 +---------------- 3 files changed, 83 insertions(+), 71 deletions(-) diff --git a/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.cpp b/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.cpp index ef8e14afaa..fca8434b1e 100644 --- a/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.cpp +++ b/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.cpp @@ -131,29 +131,36 @@ void DragTempForcing::operator()( surf_temp + thetastar / kappa * (std::log(0.5 * dx[2] / z0) - psi_h_cell); amrex::Real bc_forcing_t = -(tTarget - theta) / dt; - // //! West - // amrex::Real tmp_temp_target - // =compute_target_theta(vel(i-1,j,k,2),vel(i-1,j,k,1),theta, - // monin_obukhov_length, gravity_mod, dx[0],z0, kappa); - // bc_forcing_t += -(tmp_temp_target - theta) / dt * blank(i-1,j,k); - // //! East - // tmp_temp_target = - // compute_target_theta(vel(i+1,j,k,2),vel(i+1,j,k,1),theta, - // monin_obukhov_length, gravity_mod, dx[0],z0, kappa); - // bc_forcing_t += -(tmp_temp_target - theta) / dt * blank(i+1,j,k); - // //! South - // tmp_temp_target = - // compute_target_theta(vel(i,j-1,k,2),vel(i,j-1,k,0),theta, - // monin_obukhov_length, gravity_mod, dx[1],z0, kappa); - // bc_forcing_t += -(tmp_temp_target - theta) / dt * blank(i,j-1,k); - // //! North - // tmp_temp_target = - // compute_target_theta(vel(i,j+1,k,2),vel(i,j+1,k,0),theta, - // monin_obukhov_length, gravity_mod, dx[1],z0, kappa); - // bc_forcing_t += -(tmp_temp_target - theta) / dt * blank(i,j+1,k); - // const amrex::Real sum_blank_t = blank(i,j,k-1)+blank(i - 1, j, k) + - // blank(i + 1, j, k) + blank(i, j - 1, k) + blank(i, j + 1, k); - // bc_forcing_t /= (sum_blank_t + amr_wind::constants::EPS); + if (cell_drag > 1) { + //! West + amrex::Real tmp_temp_target = compute_target_theta( + vel(i - 1, j, k, 2), vel(i - 1, j, k, 1), theta, + monin_obukhov_length, gravity_mod, dx[0], z0, kappa); + bc_forcing_t += + -(tmp_temp_target - theta) / dt * blank(i - 1, j, k); + //! East + tmp_temp_target = compute_target_theta( + vel(i + 1, j, k, 2), vel(i + 1, j, k, 1), theta, + monin_obukhov_length, gravity_mod, dx[0], z0, kappa); + bc_forcing_t += + -(tmp_temp_target - theta) / dt * blank(i + 1, j, k); + //! South + tmp_temp_target = compute_target_theta( + vel(i, j - 1, k, 2), vel(i, j - 1, k, 0), theta, + monin_obukhov_length, gravity_mod, dx[1], z0, kappa); + bc_forcing_t += + -(tmp_temp_target - theta) / dt * blank(i, j - 1, k); + //! North + tmp_temp_target = compute_target_theta( + vel(i, j + 1, k, 2), vel(i, j + 1, k, 0), theta, + monin_obukhov_length, gravity_mod, dx[1], z0, kappa); + bc_forcing_t += + -(tmp_temp_target - theta) / dt * blank(i, j + 1, k); + const amrex::Real sum_blank_t = + blank(i, j, k - 1) + blank(i - 1, j, k) + blank(i + 1, j, k) + + blank(i, j - 1, k) + blank(i, j + 1, k); + bc_forcing_t /= (sum_blank_t + amr_wind::constants::EPS); + } 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]); diff --git a/amr-wind/equation_systems/tke/source_terms/KransAxell.cpp b/amr-wind/equation_systems/tke/source_terms/KransAxell.cpp index 32a299d126..7fd11aede5 100644 --- a/amr-wind/equation_systems/tke/source_terms/KransAxell.cpp +++ b/amr-wind/equation_systems/tke/source_terms/KransAxell.cpp @@ -6,6 +6,21 @@ #include "amr-wind/wind_energy/MOData.H" #include "amr-wind/utilities/linear_interpolation.H" #include "amr-wind/utilities/constants.H" + +namespace { + +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real compute_target_ustar( + const amrex::Real ux, + const amrex::Real uy, + const amrex::Real dx, + const amrex::Real z0, + const amrex::Real kappa) +{ + const amrex::Real wspd = std::sqrt(ux * ux + uy * uy); + const amrex::Real ustar = wspd * kappa / std::log(1.5 * dx / z0); + return ustar; +} +} // namespace namespace amr_wind::pde::tke { KransAxell::KransAxell(const CFDSim& sim) @@ -174,7 +189,7 @@ void KransAxell::operator()( 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 = + amrex::Real ustar = m * kappa / (std::log(3 * z / cell_z0) - psi_m); const amrex::Real T0 = ref_theta_arr(i, j, k); const amrex::Real hf = std::abs(gravity[2]) / T0 * heat_flux; @@ -184,6 +199,42 @@ void KransAxell::operator()( terrainforcing = (ustar * ustar / (Cmu * Cmu) + rans_b - tke_arr(i, j, k)) / dt; + if (cell_drag > 1) { + //! West + ustar = compute_target_ustar( + vel(i - 1, j, k, 2), vel(i - 1, j, k, 1), dx[0], z0, + kappa); + terrainforcing += + (ustar * ustar / (Cmu * Cmu) - tke_arr(i, j, k)) / dt * + blank_arr(i - 1, j, k); + //! East + ustar = compute_target_ustar( + vel(i + 1, j, k, 2), vel(i + 1, j, k, 1), dx[0], z0, + kappa); + terrainforcing += + (ustar * ustar / (Cmu * Cmu) - tke_arr(i, j, k)) / dt * + blank_arr(i + 1, j, k); + //! South + ustar = compute_target_ustar( + vel(i, j - 1, k, 2), vel(i, j - 1, k, 0), dx[1], z0, + kappa); + terrainforcing += + (ustar * ustar / (Cmu * Cmu) - tke_arr(i, j, k)) / dt * + blank_arr(i, j - 1, k); + //! North + ustar = compute_target_ustar( + vel(i, j + 1, k, 2), vel(i, j + 1, k, 0), dx[1], z0, + kappa); + terrainforcing += + (ustar * ustar / (Cmu * Cmu) - tke_arr(i, j, k)) / dt * + blank_arr(i, j + 1, k); + const amrex::Real sum_blank_tke = + blank_arr(i, j, k - 1) + blank_arr(i - 1, j, k) + + blank_arr(i + 1, j, k) + blank_arr(i, j - 1, k) + + blank_arr(i, j + 1, k); + terrainforcing /= + (sum_blank_tke + amr_wind::constants::EPS); + } amrex::Real bcforcing = 0; if (k == 0) { bcforcing = (1 - blank_arr(i, j, k)) * terrainforcing; diff --git a/amr-wind/physics/TerrainDrag.cpp b/amr-wind/physics/TerrainDrag.cpp index 9c883028d6..7fa3127219 100644 --- a/amr-wind/physics/TerrainDrag.cpp +++ b/amr-wind/physics/TerrainDrag.cpp @@ -141,53 +141,6 @@ void TerrainDrag::initialize_fields(int level, const amrex::Geometry& geom) levelz0[nbx](i, j, k, 0) = roughz0; }); amrex::Gpu::streamSynchronize(); - // amrex::ParallelFor( - // blanking, [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept - // { - // const amrex::GpuArray cell_blanking{ - // levelBlanking[nbx](i, j, k - 1, 0), - // levelBlanking[nbx](i - 1, j, k, 0), - // levelBlanking[nbx](i + 1, j, k, 0), - // levelBlanking[nbx](i, j - 1, k, 0), - // levelBlanking[nbx](i, j + 1, k, 0)}; - // levelDrag[nbx](i, j, k, 0) = 0; - // for (int index = 0; index <= 4; ++index) { - // if ((levelBlanking[nbx](i, j, k, 0) == 0) && - // (cell_blanking[index] == 1) && (levelDrag[nbx](i, j, k, - // 0) == 0)) { levelDrag[nbx](i, j, k, 0) = index + 1; - // } - // } - // }); - // amrex::Gpu::streamSynchronize(); - - // amrex::ParallelFor( - // blanking, - // [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { - // if ((levelBlanking[nbx](i, j, k, 0) == 0) && - // (levelBlanking[nbx](i, j, k - 1, 0) == 1)) { - // levelDrag[nbx](i, j, k, 0) = 1; - // } else { - // levelDrag[nbx](i, j, k, 0) = 0; - // } - // }); - // amrex::Gpu::streamSynchronize(); - // amrex::ParallelFor( - // blanking, [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept - // { - // const amrex::GpuArray cell_blanking{ - // levelBlanking[nbx](i - 1, j, k, 0), - // levelBlanking[nbx](i + 1, j, k, 0), - // levelBlanking[nbx](i, j - 1, k, 0), - // levelBlanking[nbx](i, j + 1, k, 0)}; - // levelDrag[nbx](i, j, k, 0) = 0; - // for (int index = 0; index <= 3; ++index) { - // if ((levelBlanking[nbx](i, j, k, 0) == 0) && - // (cell_blanking[index] == 1) ) { - // levelDrag[nbx](i, j, k, 0) = 1; - // } - // } - // }); - // amrex::Gpu::streamSynchronize(); amrex::ParallelFor( blanking, [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { const int blankxp = std::abs( @@ -214,6 +167,7 @@ void TerrainDrag::initialize_fields(int level, const amrex::Geometry& geom) : (blankxp + blankxm + blankyp + blankym + blankzp + blankzm); }); + amrex::Gpu::streamSynchronize(); } void TerrainDrag::post_init_actions() From 6970b4135df5e9dd7e878957a25cc1578b530a36 Mon Sep 17 00:00:00 2001 From: Harish Date: Mon, 23 Jun 2025 07:42:49 -0400 Subject: [PATCH 04/10] Spelling --- amr-wind/equation_systems/icns/source_terms/DragForcing.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/amr-wind/equation_systems/icns/source_terms/DragForcing.cpp b/amr-wind/equation_systems/icns/source_terms/DragForcing.cpp index f0384058ff..d42fed002c 100644 --- a/amr-wind/equation_systems/icns/source_terms/DragForcing.cpp +++ b/amr-wind/equation_systems/icns/source_terms/DragForcing.cpp @@ -315,7 +315,7 @@ void DragForcing::operator()( // BC forcing pushes nonrelative velocity toward target velocity bc_forcing_x = -(uxTarget - ux1) / dt; bc_forcing_y = -(uyTarget - uy1) / dt; - //! Adding horizonal drag + //! Adding horizontal drag if (drag(i, j, k) > 1) { //! West amrex::GpuArray tmp_wind_target = From 4df09ee75012e2cb643da8dcec58b44d510411a5 Mon Sep 17 00:00:00 2001 From: Harish Date: Thu, 26 Jun 2025 12:11:56 -0400 Subject: [PATCH 05/10] Another try --- .../icns/source_terms/DragForcing.cpp | 24 +++++++++---------- .../source_terms/DragTempForcing.cpp | 8 +++---- .../tke/source_terms/KransAxell.cpp | 12 +++++----- amr-wind/physics/TerrainDrag.cpp | 7 +++--- 4 files changed, 26 insertions(+), 25 deletions(-) diff --git a/amr-wind/equation_systems/icns/source_terms/DragForcing.cpp b/amr-wind/equation_systems/icns/source_terms/DragForcing.cpp index d42fed002c..aeb3fa3ea6 100644 --- a/amr-wind/equation_systems/icns/source_terms/DragForcing.cpp +++ b/amr-wind/equation_systems/icns/source_terms/DragForcing.cpp @@ -347,18 +347,18 @@ void DragForcing::operator()( -(tmp_wind_target[0] - uz1) / dt * blank(i, j + 1, k); bc_forcing_x += -(tmp_wind_target[1] - ux1) / dt * blank(i, j + 1, k); - const amrex::Real sum_blank_x = blank(i, j - 1, k) + - blank(i, j + 1, k) + - blank(i, j, k - 1); - bc_forcing_x /= (sum_blank_x + amr_wind::constants::EPS); - const amrex::Real sum_blank_y = blank(i - 1, j, k) + - blank(i + 1, j, k) + - blank(i, j, k - 1); - bc_forcing_y /= (sum_blank_y + amr_wind::constants::EPS); - const amrex::Real sum_blank_z = - blank(i - 1, j, k) + blank(i + 1, j, k) + - blank(i, j - 1, k) + blank(i, j + 1, k); - bc_forcing_z /= (sum_blank_z + amr_wind::constants::EPS); + // const amrex::Real sum_blank_x = blank(i, j - 1, k) + + // blank(i, j + 1, k) + + // blank(i, j, k - 1); + // bc_forcing_x /= (sum_blank_x + amr_wind::constants::EPS); + // const amrex::Real sum_blank_y = blank(i - 1, j, k) + + // blank(i + 1, j, k) + + // blank(i, j, k - 1); + // bc_forcing_y /= (sum_blank_y + amr_wind::constants::EPS); + // const amrex::Real sum_blank_z = + // blank(i - 1, j, k) + blank(i + 1, j, k) + + // blank(i, j - 1, k) + blank(i, j + 1, k); + // bc_forcing_z /= (sum_blank_z + amr_wind::constants::EPS); } } // Target velocity intended for within terrain diff --git a/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.cpp b/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.cpp index fca8434b1e..0045de55ae 100644 --- a/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.cpp +++ b/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.cpp @@ -156,10 +156,10 @@ void DragTempForcing::operator()( monin_obukhov_length, gravity_mod, dx[1], z0, kappa); bc_forcing_t += -(tmp_temp_target - theta) / dt * blank(i, j + 1, k); - const amrex::Real sum_blank_t = - blank(i, j, k - 1) + blank(i - 1, j, k) + blank(i + 1, j, k) + - blank(i, j - 1, k) + blank(i, j + 1, k); - bc_forcing_t /= (sum_blank_t + amr_wind::constants::EPS); + // const amrex::Real sum_blank_t = + // blank(i, j, k - 1) + blank(i - 1, j, k) + blank(i + 1, j, k) + + // blank(i, j - 1, k) + blank(i, j + 1, k); + // bc_forcing_t /= (sum_blank_t + amr_wind::constants::EPS); } const amrex::Real m = std::sqrt(ux1 * ux1 + uy1 * uy1 + uz1 * uz1); const amrex::Real Cd = diff --git a/amr-wind/equation_systems/tke/source_terms/KransAxell.cpp b/amr-wind/equation_systems/tke/source_terms/KransAxell.cpp index 7fd11aede5..972139c149 100644 --- a/amr-wind/equation_systems/tke/source_terms/KransAxell.cpp +++ b/amr-wind/equation_systems/tke/source_terms/KransAxell.cpp @@ -228,12 +228,12 @@ void KransAxell::operator()( terrainforcing += (ustar * ustar / (Cmu * Cmu) - tke_arr(i, j, k)) / dt * blank_arr(i, j + 1, k); - const amrex::Real sum_blank_tke = - blank_arr(i, j, k - 1) + blank_arr(i - 1, j, k) + - blank_arr(i + 1, j, k) + blank_arr(i, j - 1, k) + - blank_arr(i, j + 1, k); - terrainforcing /= - (sum_blank_tke + amr_wind::constants::EPS); + // const amrex::Real sum_blank_tke = + // blank_arr(i, j, k - 1) + blank_arr(i - 1, j, k) + + // blank_arr(i + 1, j, k) + blank_arr(i, j - 1, k) + + // blank_arr(i, j + 1, k); + // terrainforcing /= + // (sum_blank_tke + amr_wind::constants::EPS); } amrex::Real bcforcing = 0; if (k == 0) { diff --git a/amr-wind/physics/TerrainDrag.cpp b/amr-wind/physics/TerrainDrag.cpp index 7fa3127219..6ae65b73fa 100644 --- a/amr-wind/physics/TerrainDrag.cpp +++ b/amr-wind/physics/TerrainDrag.cpp @@ -162,11 +162,12 @@ void TerrainDrag::initialize_fields(int level, const amrex::Geometry& geom) levelBlanking[nbx](i, j, k - 1, 0) - levelBlanking[nbx](i, j, k, 0)); levelDrag[nbx](i, j, k, 0) = - ((k == 0 && levelBlanking[nbx](i, j, k, 0) == 1)) + ((k == 0) ? 0 : (blankxp + blankxm + blankyp + blankym + blankzp + - blankzm); - }); + blankzm)); + } + ); amrex::Gpu::streamSynchronize(); } From c80655213221cd89713434c382f4655f243d75b3 Mon Sep 17 00:00:00 2001 From: Harish Date: Thu, 26 Jun 2025 12:45:14 -0400 Subject: [PATCH 06/10] Drag --- amr-wind/physics/TerrainDrag.cpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/amr-wind/physics/TerrainDrag.cpp b/amr-wind/physics/TerrainDrag.cpp index 6ae65b73fa..30f2ab8cec 100644 --- a/amr-wind/physics/TerrainDrag.cpp +++ b/amr-wind/physics/TerrainDrag.cpp @@ -162,10 +162,9 @@ void TerrainDrag::initialize_fields(int level, const amrex::Geometry& geom) levelBlanking[nbx](i, j, k - 1, 0) - levelBlanking[nbx](i, j, k, 0)); levelDrag[nbx](i, j, k, 0) = - ((k == 0) + ((k == 0 && levelBlanking[nbx](i, j, k, 0) == 1) ) ? 0 - : (blankxp + blankxm + blankyp + blankym + blankzp + - blankzm)); + : (blankxp + blankxm + blankyp + blankym + blankzp + blankzm); } ); amrex::Gpu::streamSynchronize(); From bf6a5133837346077de9184a0fb7bfd12072d945 Mon Sep 17 00:00:00 2001 From: Harish Date: Tue, 1 Jul 2025 11:48:45 -0400 Subject: [PATCH 07/10] clang and multiple directions --- .../icns/source_terms/DragForcing.cpp | 12 ------------ .../temperature/source_terms/DragTempForcing.cpp | 6 +----- .../equation_systems/tke/source_terms/KransAxell.cpp | 8 +------- 3 files changed, 2 insertions(+), 24 deletions(-) diff --git a/amr-wind/equation_systems/icns/source_terms/DragForcing.cpp b/amr-wind/equation_systems/icns/source_terms/DragForcing.cpp index aeb3fa3ea6..bf9830cc89 100644 --- a/amr-wind/equation_systems/icns/source_terms/DragForcing.cpp +++ b/amr-wind/equation_systems/icns/source_terms/DragForcing.cpp @@ -347,18 +347,6 @@ void DragForcing::operator()( -(tmp_wind_target[0] - uz1) / dt * blank(i, j + 1, k); bc_forcing_x += -(tmp_wind_target[1] - ux1) / dt * blank(i, j + 1, k); - // const amrex::Real sum_blank_x = blank(i, j - 1, k) + - // blank(i, j + 1, k) + - // blank(i, j, k - 1); - // bc_forcing_x /= (sum_blank_x + amr_wind::constants::EPS); - // const amrex::Real sum_blank_y = blank(i - 1, j, k) + - // blank(i + 1, j, k) + - // blank(i, j, k - 1); - // bc_forcing_y /= (sum_blank_y + amr_wind::constants::EPS); - // const amrex::Real sum_blank_z = - // blank(i - 1, j, k) + blank(i + 1, j, k) + - // blank(i, j - 1, k) + blank(i, j + 1, k); - // bc_forcing_z /= (sum_blank_z + amr_wind::constants::EPS); } } // Target velocity intended for within terrain diff --git a/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.cpp b/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.cpp index 0045de55ae..96d228303f 100644 --- a/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.cpp +++ b/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.cpp @@ -131,7 +131,7 @@ void DragTempForcing::operator()( surf_temp + thetastar / kappa * (std::log(0.5 * dx[2] / z0) - psi_h_cell); amrex::Real bc_forcing_t = -(tTarget - theta) / dt; - if (cell_drag > 1) { + if (drag(i, j, k) > 1) { //! West amrex::Real tmp_temp_target = compute_target_theta( vel(i - 1, j, k, 2), vel(i - 1, j, k, 1), theta, @@ -156,10 +156,6 @@ void DragTempForcing::operator()( monin_obukhov_length, gravity_mod, dx[1], z0, kappa); bc_forcing_t += -(tmp_temp_target - theta) / dt * blank(i, j + 1, k); - // const amrex::Real sum_blank_t = - // blank(i, j, k - 1) + blank(i - 1, j, k) + blank(i + 1, j, k) + - // blank(i, j - 1, k) + blank(i, j + 1, k); - // bc_forcing_t /= (sum_blank_t + amr_wind::constants::EPS); } const amrex::Real m = std::sqrt(ux1 * ux1 + uy1 * uy1 + uz1 * uz1); const amrex::Real Cd = diff --git a/amr-wind/equation_systems/tke/source_terms/KransAxell.cpp b/amr-wind/equation_systems/tke/source_terms/KransAxell.cpp index 972139c149..afb855aef4 100644 --- a/amr-wind/equation_systems/tke/source_terms/KransAxell.cpp +++ b/amr-wind/equation_systems/tke/source_terms/KransAxell.cpp @@ -199,7 +199,7 @@ void KransAxell::operator()( terrainforcing = (ustar * ustar / (Cmu * Cmu) + rans_b - tke_arr(i, j, k)) / dt; - if (cell_drag > 1) { + if (drag(i, j, k) > 1) { //! West ustar = compute_target_ustar( vel(i - 1, j, k, 2), vel(i - 1, j, k, 1), dx[0], z0, @@ -228,12 +228,6 @@ void KransAxell::operator()( terrainforcing += (ustar * ustar / (Cmu * Cmu) - tke_arr(i, j, k)) / dt * blank_arr(i, j + 1, k); - // const amrex::Real sum_blank_tke = - // blank_arr(i, j, k - 1) + blank_arr(i - 1, j, k) + - // blank_arr(i + 1, j, k) + blank_arr(i, j - 1, k) + - // blank_arr(i, j + 1, k); - // terrainforcing /= - // (sum_blank_tke + amr_wind::constants::EPS); } amrex::Real bcforcing = 0; if (k == 0) { From ba9a4c5b42bbe115f34879d7635735f124a667bd Mon Sep 17 00:00:00 2001 From: Harish Date: Tue, 1 Jul 2025 11:52:34 -0400 Subject: [PATCH 08/10] clang --- amr-wind/equation_systems/tke/source_terms/KransAxell.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/amr-wind/equation_systems/tke/source_terms/KransAxell.cpp b/amr-wind/equation_systems/tke/source_terms/KransAxell.cpp index afb855aef4..2a1cb2e999 100644 --- a/amr-wind/equation_systems/tke/source_terms/KransAxell.cpp +++ b/amr-wind/equation_systems/tke/source_terms/KransAxell.cpp @@ -199,7 +199,7 @@ void KransAxell::operator()( terrainforcing = (ustar * ustar / (Cmu * Cmu) + rans_b - tke_arr(i, j, k)) / dt; - if (drag(i, j, k) > 1) { + if (drag_arr(i, j, k) > 1) { //! West ustar = compute_target_ustar( vel(i - 1, j, k, 2), vel(i - 1, j, k, 1), dx[0], z0, From f003446b387fa5ecce2bfdcc7dd133ef4752c9d9 Mon Sep 17 00:00:00 2001 From: Harish Date: Tue, 1 Jul 2025 11:57:00 -0400 Subject: [PATCH 09/10] Clang --- amr-wind/physics/TerrainDrag.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/amr-wind/physics/TerrainDrag.cpp b/amr-wind/physics/TerrainDrag.cpp index 30f2ab8cec..7fa3127219 100644 --- a/amr-wind/physics/TerrainDrag.cpp +++ b/amr-wind/physics/TerrainDrag.cpp @@ -162,11 +162,11 @@ void TerrainDrag::initialize_fields(int level, const amrex::Geometry& geom) levelBlanking[nbx](i, j, k - 1, 0) - levelBlanking[nbx](i, j, k, 0)); levelDrag[nbx](i, j, k, 0) = - ((k == 0 && levelBlanking[nbx](i, j, k, 0) == 1) ) + ((k == 0 && levelBlanking[nbx](i, j, k, 0) == 1)) ? 0 - : (blankxp + blankxm + blankyp + blankym + blankzp + blankzm); - } - ); + : (blankxp + blankxm + blankyp + blankym + blankzp + + blankzm); + }); amrex::Gpu::streamSynchronize(); } From 5abea3336841c46c9a204efe5b6f3dd77422e37d Mon Sep 17 00:00:00 2001 From: Harish Date: Tue, 1 Jul 2025 14:59:05 -0400 Subject: [PATCH 10/10] Clang --- .../temperature/source_terms/DragTempForcing.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.cpp b/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.cpp index 96d228303f..e881023a81 100644 --- a/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.cpp +++ b/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.cpp @@ -131,8 +131,8 @@ void DragTempForcing::operator()( surf_temp + thetastar / kappa * (std::log(0.5 * dx[2] / z0) - psi_h_cell); amrex::Real bc_forcing_t = -(tTarget - theta) / dt; + //! West if (drag(i, j, k) > 1) { - //! West amrex::Real tmp_temp_target = compute_target_theta( vel(i - 1, j, k, 2), vel(i - 1, j, k, 1), theta, monin_obukhov_length, gravity_mod, dx[0], z0, kappa);