diff --git a/amr-wind/equation_systems/icns/source_terms/CMakeLists.txt b/amr-wind/equation_systems/icns/source_terms/CMakeLists.txt index 959893163c..1d9c83c943 100644 --- a/amr-wind/equation_systems/icns/source_terms/CMakeLists.txt +++ b/amr-wind/equation_systems/icns/source_terms/CMakeLists.txt @@ -15,4 +15,5 @@ target_sources(${amr_wind_lib_name} PRIVATE RayleighDamping.cpp NonLinearSGSTerm.cpp DragForcing.cpp + ForestForcing.cpp ) diff --git a/amr-wind/equation_systems/icns/source_terms/ForestForcing.H b/amr-wind/equation_systems/icns/source_terms/ForestForcing.H new file mode 100644 index 0000000000..63b3743c23 --- /dev/null +++ b/amr-wind/equation_systems/icns/source_terms/ForestForcing.H @@ -0,0 +1,39 @@ +#ifndef FORESTFORCING_H +#define FORESTFORCING_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 a forcing term to represent forest terrain + * + * \ingroup icns_src + * + * + */ +class ForestForcing : public MomentumSource::Register +{ +public: + static std::string identifier() { return "ForestForcing"; } + + explicit ForestForcing(const CFDSim& sim); + + ~ForestForcing() override; + + void operator()( + const int lev, + const amrex::MFIter& mfi, + const amrex::Box& bx, + const FieldState fstate, + const amrex::Array4& src_term) const override; + +private: + const CFDSim& m_sim; + const Field& m_velocity; +}; + +} // namespace amr_wind::pde::icns + +#endif diff --git a/amr-wind/equation_systems/icns/source_terms/ForestForcing.cpp b/amr-wind/equation_systems/icns/source_terms/ForestForcing.cpp new file mode 100644 index 0000000000..79a723701b --- /dev/null +++ b/amr-wind/equation_systems/icns/source_terms/ForestForcing.cpp @@ -0,0 +1,42 @@ +#include "amr-wind/equation_systems/icns/source_terms/ForestForcing.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 { + +ForestForcing::ForestForcing(const CFDSim& sim) + : m_sim(sim), m_velocity(sim.repo().get_field("velocity")) +{} + +ForestForcing::~ForestForcing() = default; + +void ForestForcing::operator()( + const int lev, + const amrex::MFIter& mfi, + const amrex::Box& bx, + const FieldState fstate, + const amrex::Array4& src_term) const +{ + const auto& vel = + m_velocity.state(field_impl::dof_state(fstate))(lev).const_array(mfi); + const bool has_forest = this->m_sim.repo().field_exists("forest_drag"); + if (!has_forest) { + amrex::Abort("Need a forest to use this source term"); + } + auto* const m_forest_drag = &this->m_sim.repo().get_field("forest_drag"); + const auto& forest_drag = (*m_forest_drag)(lev).const_array(mfi); + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + const amrex::Real ux = vel(i, j, k, 0); + const amrex::Real uy = vel(i, j, k, 1); + const amrex::Real uz = vel(i, j, k, 2); + const amrex::Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz); + src_term(i, j, k, 0) -= forest_drag(i, j, k) * ux * windspeed; + src_term(i, j, k, 1) -= forest_drag(i, j, k) * uy * windspeed; + src_term(i, j, k, 2) -= forest_drag(i, j, k) * uz * windspeed; + }); +} + +} // namespace amr_wind::pde::icns diff --git a/amr-wind/physics/CMakeLists.txt b/amr-wind/physics/CMakeLists.txt index 841a4eefdd..7bc1ecab64 100644 --- a/amr-wind/physics/CMakeLists.txt +++ b/amr-wind/physics/CMakeLists.txt @@ -19,6 +19,7 @@ target_sources(${amr_wind_lib_name} TerrainDrag.cpp ActuatorSourceTagging.cpp Intermittency.cpp + ForestDrag.cpp ) add_subdirectory(multiphase) diff --git a/amr-wind/physics/ForestDrag.H b/amr-wind/physics/ForestDrag.H new file mode 100644 index 0000000000..9433981edf --- /dev/null +++ b/amr-wind/physics/ForestDrag.H @@ -0,0 +1,120 @@ +#ifndef ForestDrag_H +#define ForestDrag_H + +#include "amr-wind/core/Physics.H" +#include "amr-wind/core/Field.H" +#include "amr-wind/CFDSim.H" +#include "amr-wind/utilities/index_operations.H" +#include "amr-wind/utilities/integrals.H" +#include "amr-wind/utilities/constants.H" + +namespace amr_wind::forestdrag { + +struct Forest +{ + int m_id{-1}; + amrex::Real m_type_forest; + amrex::Real m_x_forest; + amrex::Real m_y_forest; + amrex::Real m_height_forest; + amrex::Real m_diameter_forest; + amrex::Real m_cd_forest; + amrex::Real m_lai_forest; + amrex::Real m_laimax_forest; + + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real lm() const + { + amrex::Real treelaimax = 0.0; + if (m_type_forest == 2) { + const amrex::Real treeZm = m_laimax_forest * m_height_forest; + const int n = 100; + const amrex::Real int_0_fh = amr_wind::utils::trapz( + 0, m_height_forest - constants::TIGHT_TOL, n, + [=](const amrex::Real x) noexcept { + const amrex::Real ratio = + (m_height_forest - treeZm) / (m_height_forest - x); + const auto exponent = (x < treeZm) ? 6.0 : 0.5; + return std::pow(ratio, exponent) * + std::exp(exponent * (1 - ratio)); + }); + treelaimax = m_lai_forest / int_0_fh; + } + return treelaimax; + } + + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real + area_fraction(const amrex::Real z, const amrex::Real treelaimax) const + { + amrex::Real af = 0.0; + if (m_type_forest == 1) { + af = m_lai_forest / m_height_forest; + } else if (m_type_forest == 2) { + const auto treeZm = m_laimax_forest * m_height_forest; + const auto ratio = + (m_height_forest - treeZm) / (m_height_forest - z); + const auto exponent = (z < treeZm) ? 6.0 : 0.5; + af = treelaimax * std::pow(ratio, exponent) * + std::exp(exponent * (1 - ratio)); + } + return af; + } + + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::RealBox real_bounding_box( + const amrex::GpuArray& prob_lo) const + { + const amrex::Real search_tol = 1.1; + const amrex::Real search_radius = 0.5 * search_tol * m_diameter_forest; + const auto x0 = m_x_forest - search_radius; + const auto y0 = m_y_forest - search_radius; + const auto z0 = prob_lo[2]; + const auto x1 = m_x_forest + search_radius; + const auto y1 = m_y_forest + search_radius; + const auto z1 = m_height_forest; + return {x0, y0, z0, x1, y1, z1}; + } + + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Box + bounding_box(const amrex::Geometry& geom) const + { + return utils::realbox_to_box( + real_bounding_box(geom.ProbLoArray()), geom); + } +}; + +/** Forestdrag Flow physics + * \ingroup physics + */ + +class ForestDrag : public Physics::Register +{ +public: + static std::string identifier() { return "ForestDrag"; } + + explicit ForestDrag(CFDSim& sim); + + ~ForestDrag() override = default; + + void + initialize_fields(int /*level*/, const amrex::Geometry& /*geom*/) override; + + void pre_init_actions() override {} + + void post_init_actions() override {} + + void post_regrid_actions() override; + + void pre_advance_work() override {} + + void post_advance_work() override {} + + amrex::Vector read_forest(const int level) const; + +private: + CFDSim& m_sim; + Field& m_forest_drag; + Field& m_forest_id; + std::string m_forest_file{"forest.amrwind"}; +}; +} // namespace amr_wind::forestdrag + +#endif diff --git a/amr-wind/physics/ForestDrag.cpp b/amr-wind/physics/ForestDrag.cpp new file mode 100644 index 0000000000..be6e3b2490 --- /dev/null +++ b/amr-wind/physics/ForestDrag.cpp @@ -0,0 +1,126 @@ +#include "amr-wind/physics/ForestDrag.H" +#include "amr-wind/CFDSim.H" +#include "AMReX_iMultiFab.H" +#include "AMReX_MultiFabUtil.H" +#include "AMReX_ParmParse.H" +#include "AMReX_ParReduce.H" +#include "amr-wind/utilities/trig_ops.H" +#include "amr-wind/utilities/IOManager.H" + +namespace amr_wind::forestdrag { + +ForestDrag::ForestDrag(CFDSim& sim) + : m_sim(sim) + , m_forest_drag(sim.repo().declare_field("forest_drag", 1, 1, 1)) + , m_forest_id(sim.repo().declare_field("forest_id", 1, 1, 1)) +{ + + amrex::ParmParse pp(identifier()); + pp.query("forest_file", m_forest_file); + + m_sim.io_manager().register_output_var("forest_drag"); + m_sim.io_manager().register_output_var("forest_id"); + + m_forest_drag.setVal(0.0); + m_forest_id.setVal(-1.0); + m_forest_id.set_default_fillpatch_bc(m_sim.time()); + m_forest_drag.set_default_fillpatch_bc(m_sim.time()); +} + +void ForestDrag::initialize_fields(int level, const amrex::Geometry& geom) +{ + BL_PROFILE("amr-wind::" + this->identifier() + "::initialize_fields"); + + const auto forests = read_forest(level); + amrex::Gpu::DeviceVector d_forests(forests.size()); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, forests.begin(), forests.end(), + d_forests.begin()); + + const auto& dx = geom.CellSizeArray(); + const auto& prob_lo = geom.ProbLoArray(); + auto& drag = m_forest_drag(level); + auto& fst_id = m_forest_id(level); + drag.setVal(0.0); + fst_id.setVal(-1.0); +#ifdef AMREX_USE_OMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for (amrex::MFIter mfi(m_forest_drag(level)); mfi.isValid(); ++mfi) { + const auto& vbx = mfi.growntilebox(); + for (int nf = 0; nf < static_cast(forests.size()); nf++) { + const auto bxi = vbx & forests[nf].bounding_box(geom); + if (!bxi.isEmpty()) { + const auto& levelDrag = drag.array(mfi); + const auto& levelId = fst_id.array(mfi); + const auto* forests_ptr = d_forests.data(); + amrex::ParallelFor( + bxi, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + const auto x = prob_lo[0] + (i + 0.5) * dx[0]; + const auto y = prob_lo[1] + (j + 0.5) * dx[1]; + const auto z = prob_lo[2] + (k + 0.5) * dx[2]; + const auto& fst = forests_ptr[nf]; + const auto radius = std::sqrt( + (x - fst.m_x_forest) * (x - fst.m_x_forest) + + (y - fst.m_y_forest) * (y - fst.m_y_forest)); + if (z <= fst.m_height_forest && + radius <= (0.5 * fst.m_diameter_forest)) { + const auto treelaimax = fst.lm(); + levelId(i, j, k) = fst.m_id; + levelDrag(i, j, k) += + fst.m_cd_forest * + fst.area_fraction(z, treelaimax); + } + }); + } + } + } +} + +void ForestDrag::post_regrid_actions() +{ + const int nlevels = m_sim.repo().num_active_levels(); + for (int lev = 0; lev < nlevels; ++lev) { + initialize_fields(lev, m_sim.repo().mesh().Geom(lev)); + } +} + +amrex::Vector ForestDrag::read_forest(const int level) const +{ + BL_PROFILE("amr-wind::" + this->identifier() + "::read_forest"); + + std::ifstream file(m_forest_file, std::ios::in); + if (!file.good()) { + amrex::Abort("Cannot find file " + m_forest_file); + } + + //! TreeType xc yc height diameter cd lai laimax + amrex::Vector forests; + const auto& geom = m_sim.repo().mesh().Geom(level); + const auto& ba = m_sim.repo().mesh().boxArray(level); + int cnt = 0; + amrex::Real value1, value2, value3, value4, value5, value6, value7, value8; + while (file >> value1 >> value2 >> value3 >> value4 >> value5 >> value6 >> + value7 >> value8) { + Forest f; + f.m_id = cnt; + f.m_type_forest = value1; + f.m_x_forest = value2; + f.m_y_forest = value3; + f.m_height_forest = value4; + f.m_diameter_forest = value5; + f.m_cd_forest = value6; + f.m_lai_forest = value7; + f.m_laimax_forest = value8; + + // Only keep a forest if the rank owns it + const auto bx = f.bounding_box(geom); + if (ba.intersects(bx)) { + forests.push_back(f); + } + cnt++; + } + file.close(); + return forests; +} +} // namespace amr_wind::forestdrag diff --git a/amr-wind/physics/TerrainDrag.cpp b/amr-wind/physics/TerrainDrag.cpp index 92523bfb49..77f0d74923 100644 --- a/amr-wind/physics/TerrainDrag.cpp +++ b/amr-wind/physics/TerrainDrag.cpp @@ -136,6 +136,7 @@ void TerrainDrag::initialize_fields(int level, const amrex::Geometry& geom) } }); } + void TerrainDrag::post_regrid_actions() { const int nlevels = m_sim.repo().num_active_levels(); diff --git a/amr-wind/utilities/sampling/FieldNorms.H b/amr-wind/utilities/sampling/FieldNorms.H index 53290ba123..4f571e3e9b 100644 --- a/amr-wind/utilities/sampling/FieldNorms.H +++ b/amr-wind/utilities/sampling/FieldNorms.H @@ -35,7 +35,7 @@ public: //! calculate the L2 norm of a given field and component static amrex::Real - l2_norm(amr_wind::Field& field, const int comp, const bool use_mask); + l2_norm(const amr_wind::Field& field, const int comp, const bool use_mask); const amrex::Vector& var_names() const { return m_var_names; } diff --git a/amr-wind/utilities/sampling/FieldNorms.cpp b/amr-wind/utilities/sampling/FieldNorms.cpp index 7766144ab0..acdd8995bc 100644 --- a/amr-wind/utilities/sampling/FieldNorms.cpp +++ b/amr-wind/utilities/sampling/FieldNorms.cpp @@ -34,8 +34,8 @@ void FieldNorms::initialize() prepare_ascii_file(); } -amrex::Real -FieldNorms::l2_norm(amr_wind::Field& field, const int comp, const bool use_mask) +amrex::Real FieldNorms::l2_norm( + const amr_wind::Field& field, const int comp, const bool use_mask) { amrex::Real nrm = 0.0; diff --git a/docs/sphinx/theory/theory.rst b/docs/sphinx/theory/theory.rst index d1677a3c79..447b7f96cb 100644 --- a/docs/sphinx/theory/theory.rst +++ b/docs/sphinx/theory/theory.rst @@ -557,6 +557,37 @@ the orange arrow below). :align: center :width: 30% + +Forest Model +-------------- +The forest model provides an option to include the drag from forested regions to be included in the momentum equation. The +drag force is calculated as follows: + +.. math:: + + F_i= - C_d L(x,y,z) U_i | U_i | + + +Here :math:`C_d` is the coefficient of drag for the forested region and :math:`L(x,y,z)` is the leaf area density (LAD) for the +forested region. A three-dimensional model for the LAD is usually unavailable and is also cumbersome to use if there are thousands +of trees. Two different models are available as an alternative: + +.. math:: + L=\frac{LAI}{h} + +.. math:: + L(z)=L_m \left(\frac{h - z_m}{h - z}\right)^n exp\left[n \left(1 -\frac{h - z_m}{h - z}\right )\right] + +Here :math:`LAI` is the leaf area index and is available from measurements, :math:`h` is the height of the tree, :math:`z_m` is the location +of the maximum LAD, :math:`L_m` is the maximum value of LAD at :math:`z_m` and :math:`n` is a model constant with values 6 (below :math:`z_m`) and 0.5 +(above :math:`z_m`), respectively. :math:`L_m` is computed by integrating the following equation: + +.. math:: + LAI = \int_{0}^{h} L(z) dz + +The simplified model with uniform LAD is recommended for forested regions with no knowledge of the individual trees. LAI values can be used from +climate model look-up tables for different regions around the world if no local remote sensing data is available. + Navigating source code ------------------------ diff --git a/docs/sphinx/user/inputs_incflo.rst b/docs/sphinx/user/inputs_incflo.rst index 7c3267e315..d6fa77f5c1 100644 --- a/docs/sphinx/user/inputs_incflo.rst +++ b/docs/sphinx/user/inputs_incflo.rst @@ -15,6 +15,7 @@ as initial conditions and discretization options. Current implemented physics include FreeStream, SyntheticTurbulence, ABL, Actuator, RayleighTaylor, BoussinesqBubble, TaylorGreenVortex, and ScalarAdvection (which is an example of using a passive scalar advection). For multiphase simulations, the MultiPhase physics must be specified, and for forcing wave profiles into the domain, the OceanWaves physics must be specified as well. For immersed boundary forcing method TerrainDrag must be specified and the folder should include a terrain file (default name: `terrain.amrwind`, user control is `TerrainDrag.terrain_file`) file. + For including forested regions ForestDrag must be specified and the folder should include a forest file (default name: `forest.amrwind`, user control is `ForestDrag.forest_file`) file. .. input_param:: incflo.density diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 8799e0eeee..a2825f98f3 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -267,6 +267,7 @@ add_test_re(rankine) add_test_re(rankine-sym) add_test_re(terrain_box) add_test_re(terrain_box_amr) +add_test_re(forest_drag) add_test_re(box_refinement) add_test_re(cylinder_refinement) add_test_re(freestream_godunov_inout) diff --git a/test/test_files/forest_drag/forest.amrwind b/test/test_files/forest_drag/forest.amrwind new file mode 100644 index 0000000000..785959c20c --- /dev/null +++ b/test/test_files/forest_drag/forest.amrwind @@ -0,0 +1,4 @@ +1 1024 512 45 200 0.2 6 0.8 +1 1024 1024 35 200 0.2 6 0.8 +1 1024 1224 75 200 0.2 6 0.8 +2 1024 1524 120 200 0.2 10 0.8 diff --git a/test/test_files/forest_drag/forest_drag.inp b/test/test_files/forest_drag/forest_drag.inp new file mode 100644 index 0000000000..30a50f8dc8 --- /dev/null +++ b/test/test_files/forest_drag/forest_drag.inp @@ -0,0 +1,93 @@ +# Ref: Kirkil et.al."Implementation and Evaluation of Dynamic Subfilter-Scale Stress Models for Large-Eddy Simulation Using WRF", MWR, 2012. +# Geometry +geometry.prob_lo = 0.0 0.0 0.0 +geometry.prob_hi = 2048 2048 1024 +geometry.is_periodic = 1 1 0 +# Grid +amr.n_cell = 64 64 32 +amr.max_level = 0 +# Simulation control parameters +time.stop_time = 3600 +time.max_step = 10 +time.fixed_dt = 1.0 +time.cfl = 0.9895 +time.init_shrink = 0.1 +time.regrid_interval = -1 +time.plot_interval = 86400 +time.checkpoint_interval = 14400 +# incflo +incflo.physics = ABL ForestDrag +incflo.density = 1.225 +incflo.gravity = 0. 0. -9.81 # Gravitational force (3D) +incflo.velocity = 10.0 0.0 0.0 +incflo.verbose = 0 +incflo.initial_iterations = 8 +incflo.do_initial_proj = true +incflo.constant_density = true +incflo.use_godunov = true +incflo.godunov_type = "weno_z" +incflo.diffusion_type = 2 +# transport equation parameters +transport.model = ConstTransport +transport.viscosity = 0.0 +transport.laminar_prandtl = 0.7 +transport.turbulent_prandtl = 0.333 +# turbulence equation parameters +turbulence.model = Kosovic +# Changes the length scale from filter width to kappa * z /phiM and blends with LES length scale using +# an exponential function +Kosovic.surfaceRANS = true +# A power coefficient of two for the exponential term gave best result +Kosovic.surfaceRANSExp = 2 +# Writing the non-linear terms for check +Kosovic.writeTerms = false +# Height at which the length scale switches from RANS to LES. Recommended location is 1.5 dz +Kosovic.switchLoc = 48 +# Reference Monin-Obukhov length for the RANS length scale at surface. +# Can be replaced with the computed length scale in the code. However, adds code complications. +Kosovic.refMOL = -1e30 +# Experimental feature to turn off the LES model and use the numerical dissipation as a sub-grid scale model +# Initial results showed a value of 100 m improved the results by 3-5% for neutral stratification when dx<16 m +Kosovic.LESOff = 100 +# Atmospheric boundary layer +ABL.Uperiods = 12.0 +ABL.Vperiods = 12.0 +ABL.cutoff_height = 50.0 +ABL.deltaU = 1.0 +ABL.deltaV = 1.0 +ABL.kappa = .41 +ABL.normal_direction = 2 +ABL.perturb_ref_height = 50.0 +ABL.perturb_velocity = true +ABL.perturb_temperature = false +ABL.reference_temperature = 300. +ABL.stats_output_format = native +ABL.surface_roughness_z0 = 0.1 +ABL.temperature_heights = 0 800 900 1800 2700 +ABL.temperature_values = 300 300 308 311 314 +ABL.wall_shear_stress_type = local +ABL.surface_temp_flux = 0.0 +# Source +ICNS.source_terms = BoussinesqBuoyancy CoriolisForcing GeostrophicForcing RayleighDamping NonLinearSGSTerm ForestForcing +BoussinesqBuoyancy.reference_temperature = 300.0 +BoussinesqBuoyancy.thermal_expansion_coeff = 0.003333 +#CoriolisForcing.east_vector = 1.0 0.0 0.0 +#CoriolisForcing.latitude = 45.0 +#CoriolisForcing.north_vector = 0.0 1.0 0.0 +#CoriolisForcing.rotational_time_period = 86400.0 +CoriolisForcing.east_vector = 1.0 0.0 0.0 +CoriolisForcing.north_vector = 0.0 1.0 0.0 +CoriolisForcing.latitude = 90.0 +CoriolisForcing.rotational_time_period = 125663.706143592 +GeostrophicForcing.geostrophic_wind = 10 0.0 0.0 +RayleighDamping.reference_velocity = 10 0.0 0.0 +RayleighDamping.length_sloped_damping = 150 +RayleighDamping.length_complete_damping = 50 +RayleighDamping.time_scale = 9.0 +# BC +zhi.type = "slip_wall" +zhi.temperature_type = "fixed_gradient" +zhi.temperature = 0.003 +zlo.type = "wall_model" +# Hypre + diff --git a/unit_tests/aw_test_utils/MeshTest.H b/unit_tests/aw_test_utils/MeshTest.H index f43a0f9e9d..9a5184bbdf 100644 --- a/unit_tests/aw_test_utils/MeshTest.H +++ b/unit_tests/aw_test_utils/MeshTest.H @@ -6,6 +6,7 @@ #include "aw_test_utils/AmrexTest.H" #include "aw_test_utils/AmrTestMesh.H" #include "amr-wind/core/SimTime.H" +#include "amr-wind/utilities/constants.H" #include "amr-wind/utilities/tagging/CartBoxRefinement.H" namespace amr_wind_tests { diff --git a/unit_tests/wind_energy/CMakeLists.txt b/unit_tests/wind_energy/CMakeLists.txt index be874a0e37..7d03879304 100644 --- a/unit_tests/wind_energy/CMakeLists.txt +++ b/unit_tests/wind_energy/CMakeLists.txt @@ -9,6 +9,7 @@ target_sources(${amr_wind_unit_test_exe_name} PRIVATE test_abl_bc.cpp test_abl_src_timetable.cpp test_abl_terrain.cpp + test_abl_forest.cpp ) if (AMR_WIND_ENABLE_NETCDF) diff --git a/unit_tests/wind_energy/test_abl_forest.cpp b/unit_tests/wind_energy/test_abl_forest.cpp new file mode 100644 index 0000000000..6e40b637e4 --- /dev/null +++ b/unit_tests/wind_energy/test_abl_forest.cpp @@ -0,0 +1,83 @@ +#include "aw_test_utils/MeshTest.H" +#include "aw_test_utils/iter_tools.H" +#include "aw_test_utils/test_utils.H" +#include "amr-wind/physics/ForestDrag.H" +#include "amr-wind/core/field_ops.H" +#include "amr-wind/utilities/sampling/FieldNorms.H" + +namespace { +void write_forest(const std::string& fname) +{ + std::ofstream os(fname); + //! Write forest + os << "1 512 256 45 200 0.2 6 0.8 \n"; + os << "1 512 512 35 200 0.2 6 0.8 \n"; + os << "1 512 612 75 200 0.2 6 0.8 \n"; + os << "2 512 762 120 200 0.2 10 0.8 \n"; +} + +} // namespace + +namespace amr_wind_tests { + +// Testing the terrain drag reading to ensure that terrain is properly setup +class ForestTest : public MeshTest +{ +protected: + void populate_parameters() override + { + MeshTest::populate_parameters(); + // Make computational domain like ABL mesh + { + amrex::ParmParse pp("amr"); + amrex::Vector ncell{{32, 32, 16}}; + pp.addarr("n_cell", ncell); + pp.add("blocking_factor", 2); + } + + { + amrex::ParmParse pp("geometry"); + amrex::Vector probhi{{1024, 1024, 512}}; + pp.addarr("prob_hi", probhi); + } + } + std::string m_forest_fname = "forest.amrwind"; +}; + +TEST_F(ForestTest, forest) +{ + // Write target wind file + write_forest(m_forest_fname); + populate_parameters(); + initialize_mesh(); + auto& pde_mgr = sim().pde_manager(); + pde_mgr.register_icns(); + sim().init_physics(); + amrex::ParmParse pp("incflo"); + amrex::Vector physics{"forestDrag"}; + pp.addarr("physics", physics); + amr_wind::forestdrag::ForestDrag forest_drag(sim()); + const int nlevels = sim().repo().num_active_levels(); + for (int lev = 0; lev < nlevels; ++lev) { + const auto& geom = sim().repo().mesh().Geom(lev); + forest_drag.initialize_fields(lev, geom); + } + + constexpr amrex::Real n_forests = 3.0; + const auto& f_id = sim().repo().get_field("forest_id"); + const amrex::Real max_id = amr_wind::field_ops::global_max_magnitude(f_id); + EXPECT_EQ(max_id, n_forests); + + constexpr amrex::Real expected_max_drag = 0.050285714285714288; + const auto& f_drag = sim().repo().get_field("forest_drag"); + const amrex::Real max_drag = + amr_wind::field_ops::global_max_magnitude(f_drag); + EXPECT_NEAR(max_drag, expected_max_drag, amr_wind::constants::TIGHT_TOL); + + constexpr amrex::Real expected_norm_drag = 0.0030635155406915832; + const auto norm_drag = + amr_wind::field_norms::FieldNorms::l2_norm(f_drag, 0, false); + EXPECT_NEAR(norm_drag, expected_norm_drag, amr_wind::constants::TIGHT_TOL); +} + +} // namespace amr_wind_tests