Added terrain refinement option in tagging

amr-wind/utilities/tagging/CMakeLists.txt

@@ -9,6 +9,7 @@ target_sources(${amr_wind_lib_name} PRIVATE
   VorticityMagRefinement.cpp
   OversetRefinement.cpp
   GeometryRefinement.cpp
+  TerrainRefinement.cpp
 
   # Geometry based refinement options
   BoxRefiner.cpp

amr-wind/utilities/tagging/TerrainRefinement.H

@@ -0,0 +1,57 @@
+#ifndef FIELDREFINEMENT_H
+#define FIELDREFINEMENT_H
+
+#include "amr-wind/CFDSim.H"
+#include "amr-wind/utilities/tagging/poly_ops.H"
+#include "amr-wind/utilities/tagging/RefinementCriteria.H"
+
+namespace amr_wind {
+// class CFDSim;
+class Field;
+// class IntField;
+
+/** AMR refinement using a given field (e.g., density)
+ *  \ingroup amr_utils
+ *
+ *  ```
+ *  tagging.labels = t1
+ *  tagging/t1.type = TerrainRefinement
+ *  tagging/t1.level = 3
+ *  tagging/t1.vertical_distance = 1´500
+ *  ```
+ */
+class TerrainRefinement : public RefinementCriteria::Register<TerrainRefinement>
+{
+public:
+    static std::string identifier() { return "TerrainRefinement"; }
+
+    explicit TerrainRefinement(const CFDSim& sim);
+
+    ~TerrainRefinement() override = default;
+
+    //! Read input file and initialize boxarray used to refine each level
+    void initialize(const std::string& key) override;
+
+    void operator()(
+        const int level,
+        amrex::TagBoxArray& tags,
+        const amrex::Real time,
+        const int ngrow) override;
+
+private:
+    const CFDSim& m_sim;
+
+    Field* m_terrain_height{nullptr};
+    IntField* m_terrain_blank{nullptr};
+    int m_max_lev{0};
+    amrex::Real m_vertical_distance{0.0};
+    amrex::RealBox m_tagging_box;
+
+    //! Boxarrays for each level in AMR hierarchy
+    amrex::Vector<amr_wind::polygon_utils::Point> m_poly_outer;
+    amrex::Vector<amrex::Vector<amr_wind::polygon_utils::Point>> m_poly_rings;
+};
+
+} // namespace amr_wind
+
+#endif /* TERRAINREFINEMENT_H */

amr-wind/utilities/tagging/TerrainRefinement.cpp

@@ -0,0 +1,159 @@
+#include "amr-wind/utilities/tagging/TerrainRefinement.H"
+#include "amr-wind/CFDSim.H"
+
+#include "AMReX.H"
+#include "AMReX_ParmParse.H"
+
+namespace amr_wind {
+
+TerrainRefinement::TerrainRefinement(const CFDSim& sim)
+    : m_sim(sim), m_tagging_box(m_sim.repo().mesh().Geom(0).ProbDomain())
+{}
+
+void TerrainRefinement::initialize(const std::string& key)
+{
+    amrex::ParmParse pp(key);
+
+    const auto& repo = m_sim.repo();
+
+    const bool is_terrain = repo.field_exists("terrain_height");
+    if (!is_terrain) {
+        amrex::Abort("Need terrain blanking variable to use this refinement");
+    }
+    m_terrain_height = &(m_sim.repo().get_field("terrain_height"));
+    m_terrain_blank = &(m_sim.repo().get_int_field("terrain_blank"));
+
+    // Outer radial extent of the cylinder, always read in from input file
+    pp.get("vertical_distance", m_vertical_distance);
+    pp.get("level", m_max_lev);
+    if (m_max_lev <= 0) {
+        amrex::Abort("TerrainRefinement: level should be strictly above 0");
+    }
+
+    amrex::Vector<amrex::Real> poly_outer;
+    amrex::Vector<amrex::Real> poly_inners;
+    pp.queryarr("poly_outer", poly_outer);
+
+    if (poly_outer.size() > 0) {
+        pp.getarr("poly_outer", poly_outer);
+        const int n = poly_outer.size();
+        const int n_points = static_cast<int>(poly_outer[0]);
+        // amrex::Print() << "Number of points:" << n_points << std::endl;
+        const int n_expected = n_points * 2 + 1;
+        if (n_expected != n) {
+            amrex::Abort(
+                "Expected a list of " + std::to_string(n_expected) +
+                " numbers, found " + std::to_string(n - 1) + "!");
+        }
+        m_poly_outer.resize(n_points);
+        for (int i = 0; i < n_points; ++i) {
+            const auto pt_x = poly_outer[1 + 2 * i];
+            const auto pt_y = poly_outer[1 + 2 * i + 1];
+            m_poly_outer[i] = amr_wind::polygon_utils::Point({pt_x, pt_y});
+            // amrex::Print() << pt_x << " " << pt_y << std::endl;
+        }
+        pp.queryarr("poly_inners", poly_inners);
+        if (poly_inners.size() > 0) {
+            const int n_rings = static_cast<int>(poly_inners[0]);
+            m_poly_rings.resize(n_rings);
+            // amrex::Print() << "Number of inner rings:" << n_rings <<
+            // std::endl;
+            int offset = 1;
+            for (int ring_i = 0; ring_i < n_rings; ++ring_i) {
+                const int n_pts = static_cast<int>(poly_inners[offset]);
+                m_poly_rings[ring_i].resize(n_pts);
+                // amrex::Print() << "Number of inner points:" << n_pts <<
+                // std::endl;
+                offset += 1;
+                for (int pt_i = 0; pt_i < n_pts; ++pt_i) {
+                    const auto pt_x = poly_inners[offset + 2 * pt_i];
+                    const auto pt_y = poly_inners[offset + 2 * pt_i + 1];
+                    m_poly_rings[ring_i][pt_i] =
+                        amr_wind::polygon_utils::Point({pt_x, pt_y});
+                    // amrex::Print() << m_poly_rings[ring_i][pt_i].x <<
+                    // std::endl;
+                }
+                offset += 2 * n_pts;
+            }
+            if ((n_rings > 0) && (offset != poly_inners.size())) {
+                amrex::Abort(
+                    "Expected a list of " + std::to_string(offset) +
+                    " numbers, found " + std::to_string(poly_inners.size()) +
+                    "!");
+            }
+        }
+    }
+
+    amrex::Vector<amrex::Real> box_lo(AMREX_SPACEDIM, 0);
+    amrex::Vector<amrex::Real> box_hi(AMREX_SPACEDIM, 0);
+    if (pp.queryarr("box_lo", box_lo, 0, static_cast<int>(box_lo.size())) ==
+        1) {
+        m_tagging_box.setLo(box_lo);
+    }
+    if (pp.queryarr("box_hi", box_hi, 0, static_cast<int>(box_hi.size())) ==
+        1) {
+        m_tagging_box.setHi(box_hi);
+    }
+
+    amrex::Print() << "Created terrain refinement with level " << m_max_lev
+                   << " and vertical distance " << m_vertical_distance
+                   << std::endl;
+}
+
+void TerrainRefinement::operator()(
+    int level, amrex::TagBoxArray& tags, amrex::Real time, int /*ngrow*/)
+{
+    const bool do_tag = level < m_max_lev;
+    if (!do_tag) {
+        return;
+    }
+    const auto& repo = m_sim.repo();
+    const auto& geom = repo.mesh().Geom(level);
+    const auto& prob_lo = geom.ProbLoArray();
+    const auto& dx = geom.CellSizeArray();
+    const auto tagging_box = m_tagging_box;
+
+    const auto& tag_arrs = tags.arrays();
+
+    // const auto& farrs = mfab.const_arrays();
+    const auto& mfab = (*m_terrain_height)(level);
+    const auto& mterrain_h_arrs = mfab.const_arrays();
+    const auto& mterrain_b_arrs = (*m_terrain_blank)(level).const_arrays();
+
+    amrex::ParallelFor(
+        mfab, [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept {
+            const amrex::Real z = prob_lo[2] + (k + 0.5) * dx[2];
+            const amrex::Real terrainHt = mterrain_h_arrs[nbx](i, j, k);
+            const amrex::Real cellHt = z - terrainHt;
+
+            const amrex::RealVect coord = {AMREX_D_DECL(
+                prob_lo[0] + (i + 0.5) * dx[0], prob_lo[1] + (j + 0.5) * dx[1],
+                prob_lo[2] + (k + 0.5) * dx[2])};
+
+            const auto testPt =
+                amr_wind::polygon_utils::Point({coord[0], coord[1]});
+            bool in_poly = m_poly_outer.size()
+                               ? amr_wind::polygon_utils::is_point_in_polygon(
+                                     m_poly_outer, testPt)
+                               : true;
+            if (in_poly) {
+                for (int ring_i = 0; ring_i < m_poly_rings.size(); ++ring_i) {
+                    const bool in_ring =
+                        amr_wind::polygon_utils::is_point_in_polygon(
+                            m_poly_rings[ring_i], testPt);
+                    if (in_ring) {
+                        in_poly = false;
+                        break;
+                    }
+                }
+            }
+
+            if (((cellHt >= -0.5) * dx[2] && (cellHt <= m_vertical_distance)) &&
+                (mterrain_b_arrs[nbx](i, j, k) < 1) && in_poly &&
+                (tagging_box.contains(coord))) {
+                tag_arrs[nbx](i, j, k) = amrex::TagBox::SET;
+            }
+        });
+}
+
+} // namespace amr_wind

amr-wind/utilities/tagging/poly_ops.H

@@ -0,0 +1,84 @@
+#ifndef POLY_OPS_H
+#define POLY_OPS_H
+
+#include "AMReX.H"
+#include <cmath>
+#include "AMReX_REAL.H"
+#include "AMReX_Gpu.H"
+
+// Taken from https://www.geeksforgeeks.org/point-in-polygon-in-cpp/
+namespace amr_wind::polygon_utils {
+
+struct Point
+{
+    amrex::Real x, y;
+};
+
+// Function to compute the cross product of vectors (p1p2)
+// and (p1p3)
+
+AMREX_FORCE_INLINE amrex::Real
+cross_product(const Point& p1, const Point& p2, const Point& p3)
+{
+    return (p2.x - p1.x) * (p3.y - p1.y) - (p2.y - p1.y) * (p3.x - p1.x);
+}
+
+// Function to check if point p lies on segment p1p2
+AMREX_FORCE_INLINE bool
+is_point_on_segment(const Point& p, const Point& p1, const Point& p2)
+{
+    // Check if point p lies on the line segment p1p2 and
+    // within the bounding box of p1p2
+    return cross_product(p1, p2, p) == 0 && p.x >= std::min(p1.x, p2.x) &&
+           p.x <= std::max(p1.x, p2.x) && p.y >= std::min(p1.y, p2.y) &&
+           p.y <= std::max(p1.y, p2.y);
+}
+
+// Function to compute the winding number of a point with
+// respect to a polygon
+int winding_number(const amrex::Vector<Point>& polygon, const Point& point)
+{
+    int n = polygon.size();
+    int windingNumber = 0;
+
+    // Iterate through each edge of the polygon
+    for (int i = 0; i < n; i++) {
+        Point p1 = polygon[i];
+        // Next vertex in the polygon
+        Point p2 = polygon[(i + 1) % n];
+
+        // Check if the point lies on the current edge
+        if (is_point_on_segment(point, p1, p2)) {
+            // Point is on the polygon boundary
+            return 0;
+        }
+
+        // Calculate the cross product to determine winding
+        // direction
+        if (p1.y <= point.y) {
+            if (p2.y > point.y && cross_product(p1, p2, point) > 0) {
+                windingNumber++;
+            }
+        } else {
+            if (p2.y <= point.y && cross_product(p1, p2, point) < 0) {
+                windingNumber--;
+            }
+        }
+    }
+    // Return the winding number
+    return windingNumber;
+}
+
+// Function to check if a point is inside a polygon using
+// the winding number algorithm
+bool is_point_in_polygon(
+    const amrex::Vector<Point>& polygon, const Point& point)
+{
+    // Compute the winding number for the point with respect
+    // to the polygon
+    return winding_number(polygon, point) != 0;
+}
+
+} // namespace amr_wind::polygon_utils
+
+#endif /* POLY_OPS_H */