Skip to content

Commit

Permalink
First push
Browse files Browse the repository at this point in the history
  • Loading branch information
@hgopalan
hgopalan committed Jan 5, 2025
1 parent e13a6a4 commit 1480a10
Show file tree
Hide file tree
Showing 7 changed files with 61 additions and 14 deletions.
12 changes: 7 additions & 5 deletions amr-wind/equation_systems/icns/source_terms/DragForcing.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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();
Expand Down Expand Up @@ -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));
});
Expand Down
5 changes: 4 additions & 1 deletion amr-wind/equation_systems/temperature/source_terms/DragTempForcing.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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];
Expand All @@ -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));
});
}

Expand Down
7 changes: 5 additions & 2 deletions amr-wind/equation_systems/tke/source_terms/KransAxell.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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;
});
}
}
Expand Down
3 changes: 3 additions & 0 deletions amr-wind/physics/TerrainDrag.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
26 changes: 26 additions & 0 deletions amr-wind/physics/TerrainDrag.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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();
Expand Down Expand Up @@ -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 {
Expand All @@ -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<int>((z <= terrainHt) && (z > prob_lo[2]));
levelvf[nbx](i,j,k,0)=static_cast<float>(levelBlanking[nbx](i,j,k,0));
levelheight[nbx](i, j, k, 0) =
std::max(std::abs(z - terrainHt), 0.5 * dx[2]);

Expand All @@ -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()
Expand Down
9 changes: 8 additions & 1 deletion amr-wind/turbulence/LES/Kosovic.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,6 +59,8 @@ void Kosovic<Transport>::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);
Expand Down Expand Up @@ -89,6 +91,9 @@ void Kosovic<Transport>::update_turbulent_viscosity(
const auto& blank_arr =
has_terrain ? (*m_terrain_blank)(lev).const_array(mfi)
: amrex::Array4<int>();
const auto& vf_arr = has_terrain
? (*m_terrain_vf)(lev).const_array(mfi)
: amrex::Array4<double>();
amrex::ParallelFor(
bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
const amrex::Real rho = rho_arr(i, j, k);
Expand All @@ -107,7 +112,9 @@ void Kosovic<Transport>::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 =
Expand Down
13 changes: 8 additions & 5 deletions amr-wind/turbulence/RANS/KLAxell.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -122,9 +122,12 @@ void KLAxell<Transport>::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 =
Expand Down Expand Up @@ -182,15 +185,15 @@ void KLAxell<Transport>::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);
Expand Down

0 comments on commit 1480a10

Please sign in to comment.