Skip to content

Terrain Drag Marking Update #1668

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 13 commits into
base: main
Choose a base branch
from
Open

Terrain Drag Marking Update #1668

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
27b0cb7
First
hgopalan Jun 23, 2025
f42a948
Merge branch 'Exawind:main' into multi_dir_drag
hgopalan Jun 23, 2025
85eff5a
Second
hgopalan Jun 23, 2025
106583b
Third
hgopalan Jun 23, 2025
6970b41
Spelling
hgopalan Jun 23, 2025
4df09ee
Another try
hgopalan Jun 26, 2025
c806552
Drag
hgopalan Jun 26, 2025
cf0060c
Merge branch 'Exawind:main' into multi_dir_drag
hgopalan Jun 26, 2025
eb69556
Merge branch 'main' into multi_dir_drag
hgopalan Jun 27, 2025
bf6a513
clang and multiple directions
hgopalan Jul 1, 2025
ba9a4c5
clang
hgopalan Jul 1, 2025
f003446
Clang
hgopalan Jul 1, 2025
5abea33
Clang
hgopalan Jul 1, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
70 changes: 63 additions & 7 deletions amr-wind/equation_systems/icns/source_terms/DragForcing.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,27 @@
#include "amr-wind/utilities/constants.H"

namespace {

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::GpuArray<amrex::Real, 2>
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,
Expand Down Expand Up @@ -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];
Expand Down Expand Up @@ -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;
Expand All @@ -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) {
Expand All @@ -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<amrex::Real, 2> 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 +=
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

for loop opportunity (see comments below)

-(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.;
Expand All @@ -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));
});
Expand Down
55 changes: 52 additions & 3 deletions amr-wind/equation_systems/temperature/source_terms/DragTempForcing.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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);
Expand All @@ -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(
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Another for loop opportunity

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);
});
}

Expand Down
54 changes: 50 additions & 4 deletions amr-wind/equation_systems/tke/source_terms/KransAxell.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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;
Expand All @@ -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(
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

can we do a for loop? if you cast as iv and then loop on dir and to iv +/- thedimensionvector(dir) it should do this automatically.

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;
Expand Down Expand Up @@ -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<int>(has_terrain) *
(sponge_forcing - bcforcing);
Expand Down
32 changes: 25 additions & 7 deletions amr-wind/physics/TerrainDrag.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think these includes are necessary?

namespace amr_wind::terraindrag {

namespace {} // namespace
Expand Down Expand Up @@ -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(
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

can we do a for loop? see comment above

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();
}
Expand Down
Loading