Forest Canopy Model

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
   )

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 the forcing term to include the presence of immersed boundary
+ *
+ *  \ingroup icns_src
+ *
+ *
+ */
+class ForestForcing : public MomentumSource::Register<ForestForcing>
+{
+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<amrex::Real>& src_term) const override;
+
+private:
+    const CFDSim& m_sim;
+    const Field& m_velocity;
+};
+
+} // namespace amr_wind::pde::icns
+
+#endif

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<amrex::Real>& 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

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)

amr-wind/physics/ForestDrag.H

@@ -0,0 +1,65 @@
+#ifndef ForestDrag_H
+#define ForestDrag_H
+
+#include "amr-wind/core/Physics.H"
+#include "amr-wind/core/Field.H"
+#include "amr-wind/CFDSim.H"
+
+namespace amr_wind::forestdrag {
+
+namespace {} // namespace
+
+/** Forestdrag Flow physics
+ *  \ingroup physics
+ */
+
+class ForestDrag : public Physics::Register<ForestDrag>
+{
+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 {}
+
+    int return_blank_value(int i, int j, int k);
+
+private:
+    CFDSim& m_sim;
+    const FieldRepo& m_repo;
+    const amrex::AmrCore& m_mesh;
+    Field& m_velocity;
+    //! Forest Drag Force Term
+    Field& m_forest_drag;
+    Field& m_forest_blank;
+    struct Forest
+    {
+        amrex::Real type_forest;
+        amrex::Real x_forest;
+        amrex::Real y_forest;
+        amrex::Real height_forest;
+        amrex::Real diameter_forest;
+        amrex::Real cd_forest;
+        amrex::Real lai_forest;
+        amrex::Real laimax_forest;
+    };
+    amrex::Vector<Forest> m_forests;
+    amrex::Gpu::DeviceVector<Forest> m_d_forests;
+};
+} // namespace amr_wind::forestdrag
+
+#endif

amr-wind/physics/ForestDrag.cpp

@@ -0,0 +1,140 @@
+#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 {
+
+namespace {} // namespace
+
+ForestDrag::ForestDrag(CFDSim& sim)
+    : m_sim(sim)
+    , m_repo(sim.repo())
+    , m_mesh(sim.mesh())
+    , m_velocity(sim.repo().get_field("velocity"))
+    , m_forest_drag(sim.repo().declare_field("forest_drag", 1, 1, 1))
+    , m_forest_blank(sim.repo().declare_field("forest_blank", 1, 1, 1))
+{
+    std::string forestfile("forest.amrwind");
+    amrex::ParmParse pp(identifier());
+    pp.query("forest_file", forestfile);
+    std::ifstream file(forestfile, std::ios::in);
+    if (!file.good()) {
+        amrex::Abort("Cannot find file " +forestfile);
+    }
+    //! TreeType xc yc height diameter cd lai laimax
+    amrex::Real value1, value2, value3, value4, value5, value6, value7, value8;
+    while (file >> value1 >> value2 >> value3 >> value4 >> value5 >> value6 >>
+           value7 >> value8) {
+        Forest f;
+        f.type_forest = value1;
+        f.x_forest = value2;
+        f.y_forest = value3;
+        f.height_forest = value4;
+        f.diameter_forest = value5;
+        f.cd_forest = value6;
+        f.lai_forest = value7;
+        f.laimax_forest = value8;
+        m_forests.push_back(f);
+    }
+    file.close();
+    m_d_forests.resize(m_forests.size());
+    m_sim.io_manager().register_output_var("forest_drag");
+    m_sim.io_manager().register_output_var("forest_blank");
+}
+
+void ForestDrag::post_init_actions()
+{
+    BL_PROFILE("amr-wind::" + this->identifier() + "::post_init_actions");
+    const auto& geom_vec = m_sim.repo().mesh().Geom();
+    const int nlevels = m_sim.repo().num_active_levels();
+    for (int level = 0; level < nlevels; ++level) {
+        const auto& geom = geom_vec[level];
+        const auto& dx = geom.CellSizeArray();
+        const auto& prob_lo = geom.ProbLoArray();
+        auto& velocity = m_velocity(level);
+        auto& drag = m_forest_drag(level);
+        auto& blank = m_forest_blank(level);
+        amrex::Gpu::copy(
+            amrex::Gpu::hostToDevice, m_forests.begin(), m_forests.end(),
+            m_d_forests.begin());
+        const auto* forests_ptr = m_d_forests.data();
+        const auto forestSize = m_d_forests.size();
+        for (unsigned ii = 0; ii < forestSize; ++ii) {
+            amrex::Real treelaimax = 0;
+            amrex::Real treeZm = 0.0;
+            const amrex::Real type_forest = forests_ptr[ii].type_forest;
+            const amrex::Real x_forest = forests_ptr[ii].x_forest;
+            const amrex::Real y_forest = forests_ptr[ii].y_forest;
+            const amrex::Real height_forest = forests_ptr[ii].height_forest;
+            const amrex::Real diameter_forest = forests_ptr[ii].diameter_forest;
+            const amrex::Real cd_forest = forests_ptr[ii].cd_forest;
+            const amrex::Real lai_forest = forests_ptr[ii].lai_forest;
+            const amrex::Real laimax_forest = forests_ptr[ii].laimax_forest;
+            if (type_forest == 2) {
+                amrex::Real ztree = 0;
+                amrex::Real expFun = 0;
+                amrex::Real ratio = 0;
+                treeZm = laimax_forest * height_forest;
+                const amrex::Real dz = height_forest / 100;
+                while (ztree <= height_forest) {
+                    ratio = (height_forest - treeZm) / (height_forest - ztree);
+                    if (ztree < treeZm) {
+                        expFun = expFun + std::pow(ratio, 6.0) *
+                                              std::exp(6 * (1 - ratio));
+                    } else {
+                        expFun = expFun + std::pow(ratio, 0.5) *
+                                              std::exp(0.5 * (1 - ratio));
+                    }
+                    ztree = ztree + dz;
+                }
+                treelaimax = lai_forest / (expFun * dz);
+            }
+            for (amrex::MFIter mfi(velocity); mfi.isValid(); ++mfi) {
+                const auto& vbx = mfi.validbox();
+                auto levelDrag = drag.array(mfi);
+                auto levelBlank = blank.array(mfi);
+                amrex::ParallelFor(
+                    vbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+                        // compute the source term
+                        amrex::Real af = 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 radius = std::sqrt(
+                            (x - x_forest) * (x - x_forest) +
+                            (y - y_forest) * (y - y_forest));
+                        if (z <= height_forest &&
+                            radius <= (0.5 * diameter_forest)) {
+                            if (type_forest == 1) {
+                                af = lai_forest / height_forest;
+                            } else if (type_forest == 2) {
+                                const amrex::Real ratio =
+                                    (height_forest - treeZm) /
+                                    (height_forest - z);
+                                if (z < treeZm) {
+                                    af = treelaimax * std::pow(ratio, 6.0) *
+                                         std::exp(6 * (1 - ratio));
+                                } else if (z <= height_forest) {
+                                    af = treelaimax * std::pow(ratio, 0.5) *
+                                         std::exp(0.5 * (1 - ratio));
+                                }
+                            }
+                            levelBlank(i, j, k) = ii + 1;
+                            levelDrag(i, j, k) = cd_forest * af;
+                        }
+                    });
+            }
+        }
+    }
+}
+
+void ForestDrag::pre_init_actions()
+{
+    BL_PROFILE("amr-wind::" + this->identifier() + "::pre_init_actions");
+}
+} // namespace amr_wind::forestdrag

docs/sphinx/theory/theory.rst

@@ -543,6 +543,36 @@ the orange arrow below).
    :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
 ------------------------
 

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
 

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)

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