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,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<amrex::Real, AMREX_SPACEDIM>& 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<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 {}
+
+    amrex::Vector<Forest> 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

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<Forest> 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<int>(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<Forest> 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<Forest> 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

amr-wind/physics/TerrainDrag.cpp

@@ -132,6 +132,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();

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<std::string>& var_names() const { return m_var_names; }
 

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;