Skip to content
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.

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

Terrain-Aware Immersed Forcing Method #1115

Merged
merged 101 commits into from
Oct 8, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
101 commits
Select commit Hold shift + click to select a range
d3087d5
Initial Commit of Terrain Model
hgopalan Apr 4, 2024
69320fa
Update for clang format consistency
hgopalan Apr 8, 2024
ec4086f
Coding modifications for consistency
hgopalan Apr 8, 2024
fc3724a
Clang Formatting
hgopalan Apr 8, 2024
7797847
GPU Compatibility
hgopalan Apr 8, 2024
c3ad937
Clang format
hgopalan Apr 8, 2024
bd87c48
Code update for GPU compatibility
hgopalan Apr 8, 2024
25ffb13
Update for Clang Format
hgopalan Apr 8, 2024
1b1b612
GPU compatiblity for Drag Force
hgopalan Apr 8, 2024
b735987
GPU Compatbility for Vertical Location
hgopalan Apr 8, 2024
7fca03d
GPU compatibility
hgopalan Apr 8, 2024
903ae8e
GPU Compatible Local Copy
hgopalan Apr 8, 2024
201602c
Clang
hgopalan Apr 8, 2024
6055ba7
GPU Updated
hgopalan Apr 8, 2024
7e74c79
Clang
hgopalan Apr 8, 2024
425d455
GPU Compatibility
hgopalan Apr 9, 2024
cd03664
Merge branch 'Exawind:main' into hgopalan-terrain
hgopalan Apr 9, 2024
5b401df
Merge branch 'Exawind:main' into hgopalan-terrain
hgopalan Apr 12, 2024
9dccb0d
Merge branch 'Exawind:main' into hgopalan-terrain
hgopalan May 28, 2024
88dfaf9
Merge branch 'Exawind:main' into hgopalan-terrain
hgopalan Jun 27, 2024
d871f53
Merge branch 'Exawind:main' into hgopalan-terrain
hgopalan Jul 2, 2024
db46f4a
This update modifies the code-base to make it more convenient to make…
hgopalan Jul 2, 2024
b225703
Python Utilities to create the terrain file
hgopalan Jul 2, 2024
383ce04
Updates for the clang format
hgopalan Jul 2, 2024
2e7d011
Removing the unused variables for Temperature
hgopalan Jul 2, 2024
673d25f
Formatting
hgopalan Jul 2, 2024
7e8a3db
fstate commented out
hgopalan Jul 2, 2024
287d4b1
fstate
hgopalan Jul 2, 2024
37abe4b
Lint-clang-tidy formatting
hgopalan Jul 2, 2024
41beae0
Minnor commend changes
hgopalan Jul 2, 2024
b270562
Clang Formatting
hgopalan Jul 2, 2024
702b762
Updating int to const unsigned
hgopalan Jul 2, 2024
6ba2347
Clang
hgopalan Jul 2, 2024
2e4e780
integer issues
hgopalan Jul 2, 2024
c5549c3
Formatting
hgopalan Jul 2, 2024
25840f0
unsigned
hgopalan Jul 2, 2024
34416cc
clang-tidy error for unsigned int
hgopalan Jul 2, 2024
45b6b3f
Clang Format
hgopalan Jul 2, 2024
2e22d48
Removing comment
hgopalan Jul 2, 2024
a9e2595
Added a test case for Terrain
hgopalan Jul 3, 2024
74b163e
Removed unnecessary initialization of blanking fields and converted t…
hgopalan Jul 3, 2024
890e8fe
Added Unit Test
hgopalan Jul 4, 2024
8ac13e1
Clang Formatting for Unit Test
hgopalan Jul 4, 2024
ad84719
Update test_abl_terrain.cpp
hgopalan Jul 4, 2024
63968e3
Added a terrain case.
hgopalan Jul 6, 2024
2473677
Commit to address some nomenclature changes and moving python files
hgopalan Jul 8, 2024
0e679d4
Merge branch 'main' into hgopalan-terrain
hgopalan Jul 8, 2024
c7f463b
Merge branch 'hgopalan-terrain' of https://github.com/hgopalan/amr-wi…
hgopalan Jul 8, 2024
4ecbf15
Clang Formatting
hgopalan Jul 8, 2024
02505aa
Fixing CMake Issue
hgopalan Jul 8, 2024
f81b41a
Coupled testing
hgopalan Jul 8, 2024
a481179
Merge branch 'main' into hgopalan-terrain
hgopalan Jul 9, 2024
cd73d96
Merge branch 'main' into hgopalan-terrain
hgopalan Jul 9, 2024
1693c6a
Merge branch 'Exawind:main' into hgopalan-terrain
hgopalan Jul 9, 2024
64c5d5c
Modification forc converting blanking to int field
hgopalan Jul 9, 2024
7a54646
Merge branch 'hgopalan-terrain' of https://github.com/hgopalan/amr-wi…
hgopalan Jul 9, 2024
cd01d02
Updated for consistency and formatting
hgopalan Jul 9, 2024
18ee504
Clang formatting
hgopalan Jul 9, 2024
13557f0
Formatting
hgopalan Jul 9, 2024
fbc496e
Merge branch 'Exawind:main' into hgopalan-terrain
hgopalan Jul 14, 2024
7da4b68
Merge branch 'Exawind:main' into hgopalan-terrain
hgopalan Jul 16, 2024
64c285b
Rename terrainBox.inp to terrainbox.inp
hgopalan Jul 16, 2024
ef6194b
Merge branch 'Exawind:main' into hgopalan-terrain
hgopalan Jul 16, 2024
e8e50eb
Fixing terrain_box
hgopalan Jul 16, 2024
e005a1d
Removed terrain tools from utilities. Modified forcing terms to avoid…
hgopalan Jul 17, 2024
0254ceb
Merge branch 'main' into hgopalan-terrain
marchdf Jul 19, 2024
d1a5381
Reformulated the boundary conditions to improve the surface layer pre…
hgopalan Jul 22, 2024
dea867a
Update AMRex conflict
hgopalan Jul 22, 2024
406628d
Wall Function Update
hgopalan Jul 22, 2024
f824119
Drag Force Information
hgopalan Jul 22, 2024
bebab49
Updating amrex submob
hgopalan Jul 22, 2024
332dfef
Clang Format for Wall Function
hgopalan Jul 22, 2024
f23d79a
Clang Formatting
hgopalan Jul 22, 2024
ff98ab5
Merge branch 'Exawind:main' into hgopalan-terrain
hgopalan Jul 22, 2024
36e7fe7
Clang formatting for deltaT
hgopalan Jul 23, 2024
19af0be
Kosovic Model Update
hgopalan Jul 23, 2024
512e87f
Adjusting for clang tidy
hgopalan Jul 23, 2024
03be705
Update CMakeLists.txt
hgopalan Aug 7, 2024
87d2b87
Merge branch 'Exawind:main' into hgopalan-terrain
hgopalan Aug 7, 2024
62a97a9
Update with main
hgopalan Aug 7, 2024
532bcd4
Added commit for failed test
hgopalan Aug 7, 2024
af6b0ae
Merge branch 'Exawind:main' into hgopalan-terrain
hgopalan Aug 9, 2024
c686a84
Merge branch 'Exawind:main' into hgopalan-terrain
hgopalan Aug 14, 2024
5099d33
Update CMakeLists.txt
hgopalan Aug 26, 2024
150adf0
Update CMakeLists.txt
hgopalan Aug 26, 2024
06dd7a3
Merge branch 'Exawind:main' into hgopalan-terrain
hgopalan Aug 26, 2024
bef083c
Update CMakeLists.txt
hgopalan Aug 26, 2024
1ed8544
Added damping for all the horizontal directions and made them separate.
hgopalan Aug 26, 2024
cb3843b
Clang formatting
hgopalan Aug 26, 2024
f9b249f
Adjusing the formatting for unit test
hgopalan Aug 26, 2024
a277b77
Merge branch 'Exawind:main' into hgopalan-terrain
hgopalan Sep 11, 2024
f64dce8
Merge branch 'Exawind:main' into hgopalan-terrain
hgopalan Sep 20, 2024
c1efb05
Merge branch 'Exawind:main' into hgopalan-terrain
hgopalan Oct 8, 2024
61630e2
Some code clean up for readability.
hgopalan Oct 8, 2024
c8fc54a
Updated comments for Doxygen style. Changing int to float for z0_min.
hgopalan Oct 8, 2024
cd583e6
Cosmetic changes and fixes for unit test error.
hgopalan Oct 8, 2024
4c279fc
Clang formatting.
hgopalan Oct 8, 2024
119a156
Silly error
hgopalan Oct 8, 2024
c274b0a
Merge branch 'Exawind:main' into hgopalan-terrain
hgopalan Oct 8, 2024
2161761
Removing group causing warning message.
hgopalan Oct 8, 2024
42d211c
Comments removed
hgopalan Oct 8, 2024
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
3 changes: 2 additions & 1 deletion amr-wind/equation_systems/icns/source_terms/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -13,5 +13,6 @@ target_sources(${amr_wind_lib_name} PRIVATE
ABLMesoForcingMom.cpp
BurggrafFlowForcing.cpp
RayleighDamping.cpp
NonLinearSGSTerm.cpp
NonLinearSGSTerm.cpp
DragForcing.cpp
)
55 changes: 55 additions & 0 deletions amr-wind/equation_systems/icns/source_terms/DragForcing.H
Original file line number Diff line number Diff line change
@@ -0,0 +1,55 @@
#ifndef DRAGFORCING_H
Fixed Show fixed Hide fixed
#define DRAGFORCING_H

#include "amr-wind/equation_systems/icns/MomentumSource.H"
#include "amr-wind/core/SimTime.H"
#include "amr-wind/CFDSim.H"

namespace amr_wind::pde::icns {

/** Adds the forcing term to include the presence of immersed boundary
*
* \ingroup icns_src
*
*
*/
class DragForcing : public MomentumSource::Register<DragForcing>
{
public:
static std::string identifier() { return "DragForcing"; }

explicit DragForcing(const CFDSim& sim);

~DragForcing() override;

void operator()(
const int lev,
const amrex::MFIter& mfi,
const amrex::Box& bx,
const FieldState fstate,
const amrex::Array4<amrex::Real>& src_term) const override;

private:
const SimTime& m_time;
const CFDSim& m_sim;
const amrex::AmrCore& m_mesh;
const Field& m_velocity;
amrex::Gpu::DeviceVector<amrex::Real> m_device_vel_ht;
amrex::Gpu::DeviceVector<amrex::Real> m_device_vel_vals;
amrex::Real m_drag_coefficient{10.0};
amrex::Real m_sponge_strength{1.0};
amrex::Real m_sponge_density{1.0};
amrex::Real m_sponge_distance_west{-1000};
amrex::Real m_sponge_distance_east{1000};
amrex::Real m_sponge_distance_south{-1000};
amrex::Real m_sponge_distance_north{1000};
int m_sponge_west{0};
int m_sponge_east{1};
int m_sponge_south{0};
int m_sponge_north{1};
bool m_is_laminar{false};
};

} // namespace amr_wind::pde::icns

#endif
181 changes: 181 additions & 0 deletions amr-wind/equation_systems/icns/source_terms/DragForcing.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,181 @@
#include "amr-wind/equation_systems/icns/source_terms/DragForcing.H"
#include "amr-wind/utilities/IOManager.H"

#include "AMReX_Gpu.H"
#include "AMReX_Random.H"
#include "amr-wind/wind_energy/ABL.H"

namespace amr_wind::pde::icns {

DragForcing::DragForcing(const CFDSim& sim)
: m_time(sim.time())
, m_sim(sim)
, m_mesh(sim.mesh())
, m_velocity(sim.repo().get_field("velocity"))
{
amrex::ParmParse pp("DragForcing");
pp.query("drag_coefficient", m_drag_coefficient);
pp.query("sponge_strength", m_sponge_strength);
pp.query("sponge_density", m_sponge_density);
pp.query("sponge_distance_west", m_sponge_distance_west);
pp.query("sponge_distance_east", m_sponge_distance_east);
pp.query("sponge_distance_south", m_sponge_distance_south);
pp.query("sponge_distance_north", m_sponge_distance_north);
pp.query("sponge_west", m_sponge_west);
pp.query("sponge_east", m_sponge_east);
pp.query("sponge_south", m_sponge_south);
pp.query("sponge_north", m_sponge_north);
pp.query("is_laminar", m_is_laminar);
const auto& phy_mgr = m_sim.physics_manager();
if (phy_mgr.contains("ABL")) {
const auto& abl = m_sim.physics_manager().get<amr_wind::ABL>();
const auto& fa_velocity = abl.abl_statistics().vel_profile_coarse();
m_device_vel_ht.resize(fa_velocity.line_centroids().size());
m_device_vel_vals.resize(fa_velocity.line_average().size());
amrex::Gpu::copy(
amrex::Gpu::hostToDevice, fa_velocity.line_centroids().begin(),
fa_velocity.line_centroids().end(), m_device_vel_ht.begin());
amrex::Gpu::copy(
amrex::Gpu::hostToDevice, fa_velocity.line_average().begin(),
fa_velocity.line_average().end(), m_device_vel_vals.begin());
} else {
m_sponge_strength = 0.0;
}
}

DragForcing::~DragForcing() = default;

void DragForcing::operator()(
const int lev,
const amrex::MFIter& mfi,
const amrex::Box& bx,
const FieldState fstate,
const amrex::Array4<amrex::Real>& src_term) const
{
const auto& vel =
m_velocity.state(field_impl::dof_state(fstate))(lev).const_array(mfi);
const bool is_terrain =
this->m_sim.repo().int_field_exists("terrain_blank");
if (!is_terrain) {
amrex::Abort("Need terrain blanking variable to use this source term");
}
auto* const m_terrain_blank =
&this->m_sim.repo().get_int_field("terrain_blank");
const auto& blank = (*m_terrain_blank)(lev).const_array(mfi);
auto* const m_terrain_drag =
&this->m_sim.repo().get_int_field("terrain_drag");
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& geom = m_mesh.Geom(lev);
const auto& dx = geom.CellSizeArray();
const auto& prob_lo = geom.ProbLoArray();
const auto& prob_hi = geom.ProbHiArray();
const amrex::Real drag_coefficient = m_drag_coefficient;
const amrex::Real sponge_strength = m_sponge_strength;
const amrex::Real sponge_density = m_sponge_density;
const amrex::Real start_east = prob_hi[0] - m_sponge_distance_east;
const amrex::Real start_west = prob_lo[0] - m_sponge_distance_west;
const amrex::Real start_north = prob_hi[1] - m_sponge_distance_north;
const amrex::Real start_south = prob_lo[1] - m_sponge_distance_south;
const int sponge_east = m_sponge_east;
const int sponge_west = m_sponge_west;
const int sponge_south = m_sponge_south;
const int sponge_north = m_sponge_north;
// Copy Data
const auto* device_vel_ht = m_device_vel_ht.data();
const auto* device_vel_vals = m_device_vel_vals.data();
const unsigned vsize = m_device_vel_ht.size();
const auto& dt = m_time.delta_t();
const bool is_laminar = m_is_laminar;
const amrex::Real scale_factor = (dx[2] < 1.0) ? 1.0 : 1.0 / dx[2];
const amrex::Real Cd =
(is_laminar && dx[2] < 1) ? drag_coefficient : drag_coefficient / dx[2];
const amrex::Real z0_min = 1e-4;
const auto tiny = std::numeric_limits<amrex::Real>::epsilon();
const amrex::Real kappa = 0.41;
const amrex::Real cd_max = 1000.0;
amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
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];
amrex::Real xstart_damping = 0;
amrex::Real ystart_damping = 0;
amrex::Real xend_damping = 0;
amrex::Real yend_damping = 0;
amrex::Real xi_end = (x - start_east) / (prob_hi[0] - start_east);
amrex::Real xi_start = (start_west - x) / (start_west - prob_lo[0]);
xi_start = sponge_west * std::max(xi_start, 0.0);
xi_end = sponge_east * std::max(xi_end, 0.0);
xstart_damping = sponge_west * sponge_strength * xi_start * xi_start;
xend_damping = sponge_east * sponge_strength * xi_end * xi_end;
amrex::Real yi_end = (y - start_north) / (prob_hi[1] - start_north);
amrex::Real yi_start = (start_south - y) / (start_south - prob_lo[1]);
yi_start = sponge_south * std::max(yi_start, 0.0);
yi_end = sponge_north * std::max(yi_end, 0.0);
ystart_damping = sponge_strength * yi_start * yi_start;
yend_damping = sponge_strength * yi_end * yi_end;
amrex::Real spongeVelX = 0.0;
amrex::Real spongeVelY = 0.0;
amrex::Real spongeVelZ = 0.0;
amrex::Real residual = 1000;
amrex::Real height_error = 0.0;
for (unsigned ii = 0; ii < vsize; ++ii) {
Copy link
Contributor

Choose a reason for hiding this comment

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

There might be a better way to find the first closest index to the reference height x3 using std::lower_bound and std::upper_bound or some such instead of this loop?

Copy link
Contributor

Choose a reason for hiding this comment

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

Please make sure that std::lower_bound is GPU safe. As many std functions are not.

Copy link
Contributor

Choose a reason for hiding this comment

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

Completely agreed, forgot to update in the PR, but I clarified with @hgopalan offline to not use the std functions right after I had this comment here.

height_error = std::abs(z - device_vel_ht[ii]);
if (height_error < residual) {
residual = height_error;
const unsigned ix = 3 * ii;
const unsigned iy = 3 * ii + 1;
const unsigned iz = 3 * ii + 2;
spongeVelX = device_vel_vals[ix];
spongeVelY = device_vel_vals[iy];
spongeVelZ = device_vel_vals[iz];
}
}
amrex::Real Dxz = 0.0;
amrex::Real Dyz = 0.0;
amrex::Real bc_forcing_x = 0;
amrex::Real bc_forcing_y = 0;
const amrex::Real ux1 = vel(i, j, k, 0);
const amrex::Real uy1 = vel(i, j, k, 1);
const amrex::Real uz1 = vel(i, j, k, 2);
const amrex::Real m = std::sqrt(ux1 * ux1 + uy1 * uy1 + uz1 * uz1);
if (drag(i, j, k) == 1 && (!is_laminar)) {
const amrex::Real ux2 = vel(i, j, k + 1, 0);
const amrex::Real uy2 = vel(i, j, k + 1, 1);
const amrex::Real m2 = std::sqrt(ux2 * ux2 + uy2 * uy2);
const amrex::Real z0 = std::max(terrainz0(i, j, k), z0_min);
const amrex::Real ustar = m2 * kappa / std::log(1.5 * dx[2] / z0);
const amrex::Real uTarget =
ustar / kappa * std::log(0.5 * dx[2] / z0);
const amrex::Real uxTarget =
uTarget * ux2 / (tiny + std::sqrt(ux2 * ux2 + uy2 * uy2));
const amrex::Real uyTarget =
uTarget * uy2 / (tiny + std::sqrt(ux2 * ux2 + uy2 * uy2));
bc_forcing_x = -(uxTarget - ux1) / dt;
bc_forcing_y = -(uyTarget - uy1) / dt;
Dxz = -ustar * ustar * ux1 /
(tiny + std::sqrt(ux1 * ux1 + uy1 * uy1)) / dx[2];
Dyz = -ustar * ustar * uy1 /
(tiny + std::sqrt(ux1 * ux1 + uy1 * uy1)) / dx[2];
}
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) +
(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) +
(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) +
(xstart_damping + xend_damping + ystart_damping + yend_damping) *
(uz1 - sponge_density * spongeVelZ));
});
}

} // namespace amr_wind::pde::icns
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
target_sources(${amr_wind_lib_name} PRIVATE
ABLMesoForcingTemp.cpp BodyForce.cpp
HurricaneTempForcing.cpp
DragTempForcing.cpp
)
Original file line number Diff line number Diff line change
@@ -0,0 +1,36 @@
#ifndef DRAGTEMPFORCING_H
#define DRAGTEMPFORCING_H

#include "amr-wind/equation_systems/temperature/TemperatureSource.H"
#include "amr-wind/core/SimTime.H"
#include "amr-wind/CFDSim.H"

namespace amr_wind::pde::temperature {

class DragTempForcing : public TemperatureSource::Register<DragTempForcing>
{
public:
static std::string identifier() { return "DragTempForcing"; }

explicit DragTempForcing(const CFDSim& sim);

~DragTempForcing() override;

void operator()(
const int lev,
const amrex::MFIter& mfi,
const amrex::Box& bx,
const FieldState /*fstate*/,
const amrex::Array4<amrex::Real>& src_term) const override;

private:
const CFDSim& m_sim;
const amrex::AmrCore& m_mesh;
const Field& m_velocity;
const Field& m_temperature;
amrex::Real m_drag_coefficient{1.0};
amrex::Real m_reference_temperature{300.0};
};

} // namespace amr_wind::pde::temperature
#endif
Original file line number Diff line number Diff line change
@@ -0,0 +1,63 @@

#include "amr-wind/equation_systems/temperature/source_terms/DragTempForcing.H"
#include "amr-wind/utilities/IOManager.H"

#include "AMReX_ParmParse.H"
#include "AMReX_Gpu.H"
#include "AMReX_Random.H"

namespace amr_wind::pde::temperature {

DragTempForcing::DragTempForcing(const CFDSim& sim)
: m_sim(sim)
, m_mesh(sim.mesh())
, m_velocity(sim.repo().get_field("velocity"))
, m_temperature(sim.repo().get_field("temperature"))
{
amrex::ParmParse pp("DragTempForcing");
marchdf marked this conversation as resolved.
Show resolved Hide resolved
pp.query("drag_coefficient", m_drag_coefficient);
pp.query("reference_temperature", m_reference_temperature);
}

DragTempForcing::~DragTempForcing() = default;

void DragTempForcing::operator()(
const int lev,
const amrex::MFIter& mfi,
const amrex::Box& bx,
const FieldState fstate,
const amrex::Array4<amrex::Real>& src_term) const
{
const auto& vel =
m_velocity.state(field_impl::dof_state(fstate))(lev).const_array(mfi);
const auto& temperature =
m_temperature.state(field_impl::dof_state(fstate))(lev).const_array(
mfi);
const bool is_terrain =
this->m_sim.repo().int_field_exists("terrain_blank");
if (!is_terrain) {
amrex::Abort("Need terrain blanking variable to use this source term");
}
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& geom = m_mesh.Geom(lev);
const auto& dx = geom.CellSizeArray();
const amrex::Real drag_coefficient = m_drag_coefficient / dx[2];
const amrex::Real reference_temperature = m_reference_temperature;
const auto tiny = std::numeric_limits<amrex::Real>::epsilon();
const amrex::Real cd_max = 10.0;
amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
const amrex::Real ux1 = vel(i, j, k, 0);
const amrex::Real uy1 = vel(i, j, k, 1);
const amrex::Real uz1 = vel(i, j, k, 2);
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 * (temperature(i, j, k, 0) - reference_temperature) *
blank(i, j, k, 0));
});
}

} // namespace amr_wind::pde::temperature
1 change: 1 addition & 0 deletions amr-wind/physics/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@ target_sources(${amr_wind_lib_name}
ScalarAdvection.cpp
VortexDipole.cpp
BurggrafFlow.cpp
TerrainDrag.cpp
ActuatorSourceTagging.cpp
Intermittency.cpp
)
Expand Down
Loading