diff --git a/amr-wind/equation_systems/icns/source_terms/DragForcing.cpp b/amr-wind/equation_systems/icns/source_terms/DragForcing.cpp index 80ecaa7f88..bf9830cc89 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,39 @@ void DragForcing::operator()( // BC forcing pushes nonrelative velocity toward target velocity bc_forcing_x = -(uxTarget - ux1) / dt; bc_forcing_y = -(uyTarget - uy1) / dt; + //! Adding horizontal 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); + } } // Target velocity intended for within terrain amrex::Real target_u = 0.; @@ -302,21 +358,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..e881023a81 100644 --- a/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.cpp +++ b/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.cpp @@ -6,7 +6,30 @@ #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 namespace amr_wind::pde::temperature { DragTempForcing::DragTempForcing(const CFDSim& sim) @@ -87,6 +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 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 +130,38 @@ 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 + if (drag(i, j, k) > 1) { + 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 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)); + (Cd * (theta - T0) * blank(i, j, k, 0) + 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..2a1cb2e999 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) @@ -164,16 +179,17 @@ 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); 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; @@ -183,6 +199,36 @@ void KransAxell::operator()( terrainforcing = (ustar * ustar / (Cmu * Cmu) + rans_b - tke_arr(i, j, k)) / dt; + 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, + 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); + } amrex::Real bcforcing = 0; if (k == 0) { bcforcing = (1 - blank_arr(i, j, k)) * terrainforcing; @@ -211,7 +257,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.cpp b/amr-wind/physics/TerrainDrag.cpp index dd8096aa15..7fa3127219 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 @@ -142,12 +143,29 @@ 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 { - 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(); }