diff --git a/amr-wind/utilities/tagging/CMakeLists.txt b/amr-wind/utilities/tagging/CMakeLists.txt index 163f21b8f0..290dab3c84 100644 --- a/amr-wind/utilities/tagging/CMakeLists.txt +++ b/amr-wind/utilities/tagging/CMakeLists.txt @@ -9,6 +9,8 @@ target_sources(${amr_wind_lib_name} PRIVATE VorticityMagRefinement.cpp OversetRefinement.cpp GeometryRefinement.cpp + TerrainRefinement.cpp + TerrainElevation.cpp # Geometry based refinement options BoxRefiner.cpp diff --git a/amr-wind/utilities/tagging/Polygon.H b/amr-wind/utilities/tagging/Polygon.H new file mode 100644 index 0000000000..d4a7a38e02 --- /dev/null +++ b/amr-wind/utilities/tagging/Polygon.H @@ -0,0 +1,423 @@ +#ifndef AMR_WIND_POLYGON_H_ +#define AMR_WIND_POLYGON_H_ + +#include +#include +#include +#include +#include +#include +#include + +namespace amr_wind::polygon_utils { + +/** + * \class Polygon + * \ingroup amr_utils + * \brief 2D polygon with support for holes and point-in-polygon tests. + * + * Stores all points in a single array, with offsets for each ring (outer and + * holes). This layout is GPU-friendly and matches usage in TerrainRefinement. + * + * Example input for a polygon with one hole: + * \code + * tagging.myrefine.type = TerrainRefinement + * tagging.myrefine.vertical_distance = 100 + * tagging.myrefine.level = 1 + * tagging.myrefine.poly_exterior = 0 0 0 10 10 10 10 0 0 0 + * tagging.myrefine.poly_num_holes = 1 + * tagging.myrefine.poly_hole_0 = 2 2 2 8 8 8 8 2 2 2 + * \endcode + * + * - poly_exterior: List of x y pairs for the outer ring (must be closed). + * - poly_num_holes: Number of holes. + * - poly_hole_N: List of x y pairs for each hole (must be closed). + */ +class Polygon +{ +public: + using Point = amrex::Array; + + Polygon() = default; + + /** \brief Returns true if the polygon has no points */ + bool is_empty() const { return m_points.empty(); } + + /** + * \brief Get the bounding box of the polygon + * + * \param[out] bbox_min Lower left corner of bounding box (x, y) + * \param[out] bbox_max Upper right corner of bounding box (x, y) + */ + void get_bounding_box(Point& bbox_min, Point& bbox_max) const + { + bbox_min = m_bbox_min; + bbox_max = m_bbox_max; + } + + /** + * \brief Returns true if the bounding box has been computed and is valid + * + * The bounding box is considered valid if min < max in both x and y. + * This is typically set by calling compute_bounding_box(). + * + * \return True if bounding box is valid, false otherwise + */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool bbox_valid() const + { + return (m_bbox_min[0] < m_bbox_max[0]) && + (m_bbox_min[1] < m_bbox_max[1]); + } + + /** \brief Add a vertex to the outer ring (must be called before any holes) + */ + void add_outer_vertex(const Point& pt) + { + if (m_ring_offsets.empty()) { + m_ring_offsets.push_back(0); + } + m_points.push_back(pt); + } + + /** \brief Start a new hole (inner ring) */ + void start_hole() + { + m_ring_offsets.push_back(static_cast(m_points.size())); + } + + /** \brief Add a vertex to the current ring (outer or last hole) */ + void add_vertex(const Point& pt) { m_points.push_back(pt); } + + /** \brief Return the number of rings (outer + holes) */ + int num_rings() const { return static_cast(m_ring_offsets.size()); } + + /** \brief Return the number of points in ring i */ + int ring_size(int i) const + { + if (i < 0 || i >= static_cast(m_ring_offsets.size())) { + return 0; + } + int start = m_ring_offsets[i]; + int end = (i + 1 < static_cast(m_ring_offsets.size())) + ? m_ring_offsets[i + 1] + : static_cast(m_points.size()); + return end - start; + } + + /** \brief Return a pointer to the start of ring i */ + const Point* ring_ptr(int i) const + { + if (i < 0 || i >= static_cast(m_ring_offsets.size())) { + return nullptr; + } + return m_points.data() + m_ring_offsets[i]; + } + + /** + * \brief Check if a point is inside the polygon (excluding holes) + * \param pt The point to test + * \return True if inside, false otherwise + */ + bool contains(const Point& pt) const + { + if (is_empty()) { + return false; + } + if (!bounding_box_contains(pt)) { + return false; + } + return is_point_in_polygon( + m_points.data(), m_ring_offsets.data(), num_rings(), num_points(), + pt); + } + + /** + * \brief Standalone check if a point is inside a bounding box + * \param pt The point to test + * \param bbox_min Lower left corner of bounding box + * \param bbox_max Upper right corner of bounding box + * \return True if inside bounding box, false otherwise + */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE static bool + poly_bounding_box_contains( + const Point& pt, const Point& bbox_min, const Point& bbox_max) + { + return ( + pt[0] >= bbox_min[0] && pt[0] <= bbox_max[0] && + pt[1] >= bbox_min[1] && pt[1] <= bbox_max[1]); + } + + /** + * \brief Check if a point is inside the bounding box + * \param pt The point to test + * \return True if inside bounding box, false otherwise + */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool + bounding_box_contains(const Point& pt) const + { + return poly_bounding_box_contains(pt, m_bbox_min, m_bbox_max); + } + + /** \brief Compute the bounding box of the polygon */ + void compute_bounding_box() + { + m_bbox_min = { + std::numeric_limits::max(), + std::numeric_limits::max()}; + m_bbox_max = { + std::numeric_limits::lowest(), + std::numeric_limits::lowest()}; + for (const auto& pt : m_points) { + m_bbox_min[0] = amrex::min(m_bbox_min[0], pt[0]); + m_bbox_min[1] = amrex::min(m_bbox_min[1], pt[1]); + m_bbox_max[0] = amrex::max(m_bbox_max[0], pt[0]); + m_bbox_max[1] = amrex::max(m_bbox_max[1], pt[1]); + } + } + + /** + * \brief Read polygon data from ParmParse + * \param prefix ParmParse prefix for polygon input + */ + void read_from_parmparse(const std::string& prefix) + { + amrex::ParmParse pp(prefix); + + // Outer ring + amrex::Vector outer_coords; + pp.queryarr("poly_exterior", outer_coords); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + outer_coords.size() % 2 == 0, + "Polygon exterior ring must have an even number of coordinates."); + + m_points.clear(); + m_ring_offsets.clear(); + m_ring_offsets.push_back(0); + for (amrex::Long i = 0; + i < static_cast(outer_coords.size()); i += 2) { + m_points.push_back({outer_coords[i], outer_coords[i + 1]}); + } + + // Holes + int num_holes = 0; + pp.query("poly_num_holes", num_holes); + for (int h = 0; h < num_holes; ++h) { + std::string key = "poly_hole_" + std::to_string(h); + amrex::Vector hole_coords; + pp.queryarr(key.c_str(), hole_coords); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + hole_coords.size() % 2 == 0, + "Polygon hole must have an even number of coordinates."); + + m_ring_offsets.push_back(static_cast(m_points.size())); + for (amrex::Long i = 0; + i < static_cast(hole_coords.size()); i += 2) { + m_points.push_back({hole_coords[i], hole_coords[i + 1]}); + } + } + + compute_bounding_box(); + } + + /** + * \brief Print the polygon to the given stream + * \param os Output stream (defaults to amrex::OutStream()) + */ + void print(std::ostream& os = amrex::OutStream()) const + { + if (is_empty()) { + os << "Polygon: (empty)\n"; + os << "Bounding box: not computed\n"; + return; + } + os << "Polygon (outer ring):\n"; + for (int i = 0; i < ring_size(0); ++i) { + const auto& p = ring_ptr(0)[i]; + os << " (" << p[0] << ", " << p[1] << ")\n"; + } + for (int h = 1; h < num_rings(); ++h) { + os << " Hole " << h << ":\n"; + for (int i = 0; i < ring_size(h); ++i) { + const auto& p = ring_ptr(h)[i]; + os << " (" << p[0] << ", " << p[1] << ")\n"; + } + } + + bool bbox_valid = + (m_bbox_min[0] < m_bbox_max[0]) && (m_bbox_min[1] < m_bbox_max[1]); + os << "Bounding box: "; + if (bbox_valid) { + os << "min=(" << m_bbox_min[0] << ", " << m_bbox_min[1] << "), " + << "max=(" << m_bbox_max[0] << ", " << m_bbox_max[1] << ")\n"; + } else { + os << "not computed\n"; + } + } + + /** \brief Access to raw points for GPU */ + const amrex::Vector& points() const { return m_points; } + + /** \brief Get total number of points */ + int num_points() const { return static_cast(m_points.size()); } + + /** \brief Access to ring offsets for GPU */ + const amrex::Vector& ring_offsets() const { return m_ring_offsets; } + + /** + * \brief Standalone ring point-in-polygon test for GPU use + * \param ring Pointer to ring points + * \param n Number of points in ring + * \param pt Point to test + * \return True if inside, false otherwise + */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE static bool + is_point_in_ring(const Point* ring, int n, const Point& pt) + { + int wn = 0; + for (int j = 0; j < n; ++j) { + const auto& p1 = ring[j]; + const auto& p2 = ring[(j + 1) % n]; + amrex::Real cross = (p2[0] - p1[0]) * (pt[1] - p1[1]) - + (pt[0] - p1[0]) * (p2[1] - p1[1]); + if (amrex::Math::abs(cross) < 1e-12 && + pt[0] >= amrex::min(p1[0], p2[0]) && + pt[0] <= amrex::max(p1[0], p2[0]) && + pt[1] >= amrex::min(p1[1], p2[1]) && + pt[1] <= amrex::max(p1[1], p2[1])) { + return false; // on edge = not inside + } + if (p1[1] <= pt[1]) { + if (p2[1] > pt[1] && cross > 0) { + ++wn; + } + } else { + if (p2[1] <= pt[1] && cross < 0) { + --wn; + } + } + } + return wn != 0; + } + + /** + * \brief Standalone check if a point is inside a polygon with holes, given + * raw data pointers + * \param p_poly_points Pointer to all polygon points + * \param p_ring_offsets Pointer to ring offsets + * \param n_rings Number of rings (outer + holes) + * \param n_points Total number of points + * \param pt The point to test + * \return True if strictly inside, false otherwise (false on edge or in + * hole) + */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE static bool is_point_in_polygon( + const Point* p_poly_points, + const int* p_ring_offsets, + int n_rings, + int n_points, + const Point& pt) + { + // Outer ring + int start = static_cast(p_ring_offsets[0]); + int end = (n_rings > 1) ? p_ring_offsets[1] : n_points; + int n = end - start; + // Check if pt is on edge of outer ring + for (int j = 0; j < n; ++j) { + const auto& p1 = p_poly_points[start + j]; + const auto& p2 = p_poly_points[start + ((j + 1) % n)]; + if (is_point_on_segment(pt, p1, p2)) { + return false; + } + } + if (!is_point_in_ring(p_poly_points + start, n, pt)) { + return false; + } + // Check holes + for (int ring_i = 1; ring_i < n_rings; ++ring_i) { + int h_start = p_ring_offsets[ring_i]; + int h_end = + (ring_i + 1 < n_rings) ? p_ring_offsets[ring_i + 1] : n_points; + int h_n = h_end - h_start; + // Check if pt is on edge of hole + for (int j = 0; j < h_n; ++j) { + const auto& p1 = p_poly_points[h_start + j]; + const auto& p2 = p_poly_points[h_start + ((j + 1) % h_n)]; + if (is_point_on_segment(pt, p1, p2)) { + return false; + } + } + if (is_point_in_ring(p_poly_points + h_start, h_n, pt)) { + return false; + } + } + return true; + } + +private: + amrex::Vector m_points; + amrex::Vector m_ring_offsets; + Point m_bbox_min; + Point m_bbox_max; + + /** \brief Helper: check if point pt is on segment p1-p2 */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE static bool + is_point_on_segment(const Point& pt, const Point& p1, const Point& p2) + { + amrex::Real cross = is_left(p1, p2, pt); + if (std::abs(cross) > 1e-12) { + return false; + } + return ( + pt[0] >= amrex::min(p1[0], p2[0]) && + pt[0] <= amrex::max(p1[0], p2[0]) && + pt[1] >= amrex::min(p1[1], p2[1]) && + pt[1] <= amrex::max(p1[1], p2[1])); + } + + /** \brief Compute winding number for a point with respect to ring i */ + int winding_number(int ring_idx, const Point& pt) const + { + if (ring_idx < 0 || + ring_idx >= static_cast(m_ring_offsets.size())) { + return 0; + } + int wn = 0; + int start = m_ring_offsets[ring_idx]; + int end = (ring_idx + 1 < static_cast(m_ring_offsets.size())) + ? m_ring_offsets[ring_idx + 1] + : static_cast(m_points.size()); + int n = end - start; + if (n == 0) { + return 0; + } + for (int j = 0; j < n; ++j) { + const auto& p1 = m_points[start + j]; + const auto& p2 = m_points[start + ((j + 1) % n)]; + if (is_point_on_segment(pt, p1, p2)) { + return 0; + } + if (p1[1] <= pt[1]) { + if (p2[1] > pt[1] && is_left(p1, p2, pt) > 0) { + ++wn; + } + } else { + if (p2[1] <= pt[1] && is_left(p1, p2, pt) < 0) { + --wn; + } + } + } + return wn; + } + + /** \brief Helper for winding number: is pt left of line p0->p1? */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE static amrex::Real + is_left(const Point& p0, const Point& p1, const Point& p2) + { + return (p1[0] - p0[0]) * (p2[1] - p0[1]) - + (p2[0] - p0[0]) * (p1[1] - p0[1]); + } +}; + +} // namespace amr_wind::polygon_utils + +#endif // AMR_WIND_POLYGON_H_ \ No newline at end of file diff --git a/amr-wind/utilities/tagging/TerrainElevation.H b/amr-wind/utilities/tagging/TerrainElevation.H new file mode 100644 index 0000000000..41337803b3 --- /dev/null +++ b/amr-wind/utilities/tagging/TerrainElevation.H @@ -0,0 +1,74 @@ +/* + * NOTE: This class is a work in progress. + * In the future, TerrainElevation may be extended to handle additional + * terrain-related data such as roughness, and could be used by other components + * of the codebase, including TerrainDrag, MetMastForcing, and similar modules. + */ + +#ifndef TERRAINELEVATION_H_ +#define TERRAINELEVATION_H_ + +#include "AMReX.H" +#include +#include +#include + +namespace amr_wind { + +class TerrainElevation +{ +public: + TerrainElevation() = default; + + void load_from_file(const std::string& filename); + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool is_loaded() const + { + return m_loaded; + } + + static std::shared_ptr& get_instance(); + static void ensure_loaded(const std::string& filename); + + // Accessors (const and non-const) + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE const amrex::Vector& + x() const + { + return m_x; + } + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE const amrex::Vector& + y() const + { + return m_y; + } + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE const amrex::Vector& + z() const + { + return m_z; + } + + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Vector& x() + { + return m_x; + } + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Vector& y() + { + return m_y; + } + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Vector& z() + { + return m_z; + } + +private: + amrex::Vector m_x; + amrex::Vector m_y; + amrex::Vector m_z; + bool m_loaded{false}; +}; + +// Singleton accessor +std::shared_ptr& get_terrain_elevation(); + +} // namespace amr_wind + +#endif /* TERRAINELEVATION_H_ */ \ No newline at end of file diff --git a/amr-wind/utilities/tagging/TerrainElevation.cpp b/amr-wind/utilities/tagging/TerrainElevation.cpp new file mode 100644 index 0000000000..425f6e187b --- /dev/null +++ b/amr-wind/utilities/tagging/TerrainElevation.cpp @@ -0,0 +1,31 @@ +#include +#include "amr-wind/utilities/tagging/TerrainElevation.H" +#include "amr-wind/utilities/io_utils.H" + +namespace amr_wind { + +void TerrainElevation::load_from_file(const std::string& filename) +{ + ioutils::read_flat_grid_file(filename, m_x, m_y, m_z); + m_loaded = !m_x.empty() && !m_y.empty() && !m_z.empty(); + amrex::Print() << "TerrainElevation loaded from file: " << filename << "\n"; +} + +std::shared_ptr& TerrainElevation::get_instance() +{ + static std::shared_ptr instance; + return instance; +} + +void TerrainElevation::ensure_loaded(const std::string& filename) +{ + auto& ptr = get_instance(); + if (!ptr) { + ptr = std::make_shared(); + } + if (!ptr->is_loaded()) { + ptr->load_from_file(filename); + } +} + +} // namespace amr_wind diff --git a/amr-wind/utilities/tagging/TerrainRefinement.H b/amr-wind/utilities/tagging/TerrainRefinement.H new file mode 100644 index 0000000000..144aed23d3 --- /dev/null +++ b/amr-wind/utilities/tagging/TerrainRefinement.H @@ -0,0 +1,171 @@ +#ifndef TERRAINREFINEMENT_H +#define TERRAINREFINEMENT_H + +#include "amr-wind/CFDSim.H" +#include "amr-wind/utilities/tagging/RefinementCriteria.H" +#include "amr-wind/utilities/tagging/Polygon.H" +#include +#include +#include + +namespace amr_wind::format_utils { + +template +inline std::string list_to_string( + const std::vector& vec, + const std::string& delimiter = ", ", + const bool use_parentheses = false) +{ + std::ostringstream oss; + if (use_parentheses) { + oss << "("; + } + for (size_t i = 0; i < static_cast(vec.size()); ++i) { + if (i > 0) { + oss << delimiter; + } + oss << vec[i]; + } + if (use_parentheses) { + oss << ")"; + } + return oss.str(); +} + +} // namespace amr_wind::format_utils + +namespace amr_wind { + +class Field; + +/** + * \class TerrainRefinement + * \ingroup amr_utils + * \brief AMR refinement using a given distance to terrain and optional polygon + * region. + * + * Example usage in input file: + * \code + * tagging.labels = terrain_lvl1 + * tagging.terrain_lvl1.type = TerrainRefinement + * tagging.terrain_lvl1.vertical_distance = 100 + * tagging.terrain_lvl1.max_level = 1 + * tagging.terrain_lvl1.grid_buffer_ratio_lo = 0.0 0.5 + * tagging.terrain_lvl1.grid_buffer_ratio_hi = 0.0 0.5 + * tagging.terrain_lvl1.verbose = 0 1 + * tagging.terrain_lvl1.poly_exterior = 0 0 0 10 10 10 10 0 0 0 + * tagging.terrain_lvl1.poly_num_holes = 1 + * tagging.terrain_lvl1.poly_hole_0 = 2 2 2 8 8 8 8 2 2 2 + * tagging.terrain_lvl1.box_lo = 0 0 0 + * tagging.terrain_lvl1.box_hi = 10 10 10 + * \endcode + * + * - Only cells within the polygon (if specified), tagging box, and vertical + * distance above terrain are refined. + * - If no polygon is specified, the refinement applies to the whole tagging + * box. + * - The grid_buffer_ratio_lo and grid_buffer_ratio_hi parameters control the + * buffer region below and above the terrain, respectively, as a fraction of the + * local tile height. Making these ratios larger ensures more cells are tagged + * for refinement. + * + */ + +class TerrainRefinement : public RefinementCriteria::Register +{ +public: + static std::string identifier() { return "TerrainRefinement"; } + + explicit TerrainRefinement(const CFDSim& sim); + + ~TerrainRefinement() override = default; + + /** + * \brief Read input file and initialize boxarray used to refine each level + * \param key ParmParse prefix for this refinement + */ + void initialize(const std::string& key) override; + + /** + * \brief Tag cells for refinement + * \param level AMR level + * \param tags TagBoxArray to set + * \param time Current simulation time + * \param ngrow Number of ghost cells + */ + void operator()( + const int level, + amrex::TagBoxArray& tags, + const amrex::Real time, + const int ngrow) override; + +private: + const CFDSim& m_sim; + + //! Terrain file + std::string m_terrain_file{"terrain.amrwind"}; + + //! Ratios of buffer size below the terrain (as a fraction of tile height) + amrex::Vector m_grid_buffer_ratio_lo; + + //! Ratios of buffer size above the terrain (as a fraction of tile height) + amrex::Vector m_grid_buffer_ratio_hi; + + //! Key name for this refinement criteria + std::string m_key; + + //! Only refine this level if set >= 0 + int m_set_level{-1}; + + //! Minimum level to refine (if m_set_level < 0) + int m_min_level{0}; + + //! Maximum level to refine (if m_set_level < 0) + int m_max_level{32}; + + //! Distance above the terrain to refine + amrex::Real m_vertical_distance{0.0}; + + //! The usual tagging bbox */ + amrex::RealBox m_tagging_box; + + //! Optional: The coordinates of the polygon's exterior ring + amr_wind::polygon_utils::Polygon m_polygon; + + //! Verbose output levels + amrex::Vector m_verbose_levels; + +#ifdef AMREX_USE_GPU + amrex::Gpu::DeviceVector + m_poly_points_dv; + amrex::Gpu::DeviceVector m_ring_offsets_dv; +#endif + + //! Get active levels as a vector of integers + AMREX_FORCE_INLINE amrex::Vector get_active_levels() const; + + //! Check if the level should be tagged + AMREX_FORCE_INLINE bool should_tag_level(int level) const + { + if (m_set_level >= 0) { + return (level == m_set_level); + } + return (level >= m_min_level && level <= m_max_level); + } + + //! Check if the level is verbose + AMREX_FORCE_INLINE bool is_level_verbose(int level) const + { + if (m_verbose_levels.empty()) { + return false; + } + int v = (level < static_cast(m_verbose_levels.size())) + ? m_verbose_levels[level] + : m_verbose_levels.back(); + return v > 0; + } +}; + +} // namespace amr_wind + +#endif /* TERRAINREFINEMENT_H */ diff --git a/amr-wind/utilities/tagging/TerrainRefinement.cpp b/amr-wind/utilities/tagging/TerrainRefinement.cpp new file mode 100644 index 0000000000..82c8221588 --- /dev/null +++ b/amr-wind/utilities/tagging/TerrainRefinement.cpp @@ -0,0 +1,298 @@ +#include "amr-wind/CFDSim.H" +#include "AMReX.H" +#include "AMReX_Gpu.H" +#include "AMReX_ParmParse.H" +#include "amr-wind/physics/TerrainDrag.H" +#include "amr-wind/utilities/tagging/TerrainRefinement.H" +#include "amr-wind/utilities/io_utils.H" +#include "amr-wind/utilities/linear_interpolation.H" +#include "amr-wind/utilities/tagging/TerrainElevation.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) +{ + m_key = key.substr(8); // Store the key name + amrex::ParmParse pp(key); + pp.queryarr("verbose", m_verbose_levels); + pp.queryarr("grid_buffer_ratio_lo", m_grid_buffer_ratio_lo); + pp.queryarr("grid_buffer_ratio_hi", m_grid_buffer_ratio_hi); + + //! Reading the Terrain Coordinates from file + amrex::ParmParse pp_terrain("TerrainDrag"); + pp_terrain.query("terrain_file", m_terrain_file); + amr_wind::TerrainElevation::ensure_loaded(m_terrain_file); + + // Outer radial extent of the cylinder, always read in from input file + pp.get("vertical_distance", m_vertical_distance); + + // Get the levels + pp.query("level", m_set_level); + if (m_set_level >= 0) { + m_min_level = m_set_level; + m_max_level = m_set_level; + } else { + pp.query("min_level", m_min_level); + pp.query("max_level", m_max_level); + } + + // Read polygon from input file + m_polygon.read_from_parmparse(key); + + amrex::Vector box_lo(AMREX_SPACEDIM, 0); + amrex::Vector box_hi(AMREX_SPACEDIM, 0); + if (pp.queryarr("box_lo", box_lo, 0, static_cast(box_lo.size())) == + 1) { + m_tagging_box.setLo(box_lo); + } + if (pp.queryarr("box_hi", box_hi, 0, static_cast(box_hi.size())) == + 1) { + m_tagging_box.setHi(box_hi); + } + +#ifdef AMREX_USE_GPU + m_poly_points_dv.resize(m_polygon.points().size()); + amrex::Gpu::copyAsync( + amrex::Gpu::hostToDevice, m_polygon.points().begin(), + m_polygon.points().end(), m_poly_points_dv.begin()); + + m_ring_offsets_dv.resize(m_polygon.ring_offsets().size()); + amrex::Gpu::copyAsync( + amrex::Gpu::hostToDevice, m_polygon.ring_offsets().begin(), + m_polygon.ring_offsets().end(), m_ring_offsets_dv.begin()); +#endif + + auto active_levels = + amr_wind::format_utils::list_to_string(get_active_levels(), ", ", true); + + amrex::Print() << "Created terrain refinement " << m_key << " for levels " + << active_levels << " and vertical distance " + << m_vertical_distance << " (grid_buffer_ratio_lo: " + << amr_wind::format_utils::list_to_string( + m_grid_buffer_ratio_lo, ", ", true) + << ", grid_buffer_ratio_hi: " + << amr_wind::format_utils::list_to_string( + m_grid_buffer_ratio_hi, ", ", true) + << ")" << std::endl; +} + +void TerrainRefinement::operator()( + int level, amrex::TagBoxArray& tags, amrex::Real /*time*/, int /*ngrow*/) +{ + const bool verbose = is_level_verbose(level); + + if (!should_tag_level(level)) { + return; + } + + // Get the buffer ratios for this level + auto buffer_lo = m_grid_buffer_ratio_lo.empty() + ? 0.0 + : (level < m_grid_buffer_ratio_lo.size() + ? m_grid_buffer_ratio_lo[level] + : m_grid_buffer_ratio_lo.back()); + auto buffer_hi = m_grid_buffer_ratio_hi.empty() + ? 0.0 + : (level < m_grid_buffer_ratio_hi.size() + ? m_grid_buffer_ratio_hi[level] + : m_grid_buffer_ratio_hi.back()); + + // Geometry + const auto& geom = m_sim.repo().mesh().Geom(level); + const auto& prob_lo = geom.ProbLoArray(); + const auto& dx = geom.CellSizeArray(); + + // Get tagging box and vertical distance for device + const auto tagging_box = m_tagging_box; + + // Round up to the nearest positive multiple of dx[2] (if not already a + // multiple) + auto vertical_distance = m_vertical_distance; + if (vertical_distance <= 0.0) { + vertical_distance = dx[2]; + } else { + vertical_distance = std::ceil(vertical_distance / dx[2]) * dx[2]; + } + + // Tag arrays + const auto& tag_arrs = tags.arrays(); + + const auto& mfab = (*m_sim.repo().fields()[0])(level); + + // Prepare polygon data for GPU + const bool polygon_is_empty = m_polygon.is_empty(); + auto n_rings = m_polygon.num_rings(); + auto n_points = m_polygon.points().size(); + +#ifdef AMREX_USE_GPU + auto const* p_poly_points = m_poly_points_dv.data(); + auto const* p_ring_offsets = m_ring_offsets_dv.data(); +#else + auto const* p_poly_points = m_polygon.points().data(); + auto const* p_ring_offsets = m_polygon.ring_offsets().data(); +#endif + + // Get polygon bounding box (lo, hi) + amr_wind::polygon_utils::Polygon::Point bbox_lo = {0.0, 0.0}; + amr_wind::polygon_utils::Polygon::Point bbox_hi = {0.0, 0.0}; + if (!polygon_is_empty) { + if (!m_polygon.bbox_valid()) { + m_polygon.compute_bounding_box(); + } + m_polygon.get_bounding_box(bbox_lo, bbox_hi); + } + + auto terrain_ptr = amr_wind::TerrainElevation::get_instance(); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + terrain_ptr && terrain_ptr->is_loaded(), + "TerrainElevation must be loaded before TerrainRefinement is used."); + + const auto& xterrain = terrain_ptr->x(); + const auto& yterrain = terrain_ptr->y(); + const auto& zterrain = terrain_ptr->z(); + + const auto xterrain_size = xterrain.size(); + const auto yterrain_size = yterrain.size(); + const auto zterrain_size = zterrain.size(); + amrex::Gpu::DeviceVector d_xterrain(xterrain_size); + amrex::Gpu::DeviceVector d_yterrain(yterrain_size); + amrex::Gpu::DeviceVector d_zterrain(zterrain_size); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, xterrain.begin(), xterrain.end(), + d_xterrain.begin()); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, yterrain.begin(), yterrain.end(), + d_yterrain.begin()); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, zterrain.begin(), zterrain.end(), + d_zterrain.begin()); + const auto* xterrain_ptr = d_xterrain.data(); + const auto* yterrain_ptr = d_yterrain.data(); + const auto* zterrain_ptr = d_zterrain.data(); + + for (amrex::MFIter mfi(mfab, amrex::TilingIfNotGPU()); mfi.isValid(); + ++mfi) { + const auto& bx = mfi.tilebox(); + int nbx = mfi.LocalIndex(); + + // Compute buffer height based on the vertical distance and buffer ratio + // This buffer is used to ensure that we tag cells that are slightly + // above the terrain, allowing for a smoother transition in terrain + // refinement. + // The buffer is a fraction of the box height in the z-direction + // and is defined as a ratio of the vertical distance. + // This ensures that the buffer is proportional to the vertical distance + // and adapts to the size of the box in the z-direction. + const auto box_height = bx.length(2) * dx[2]; + auto buffer_z_lo = buffer_lo * box_height; + buffer_z_lo = std::max(buffer_z_lo, 0.5 * dx[2]); + auto buffer_z_hi = buffer_hi * box_height; + + // Compute tilebox bounds in physical space + auto tile_lo_x = prob_lo[0] + (bx.smallEnd(0) + 0.0) * dx[0]; + auto tile_hi_x = prob_lo[0] + (bx.bigEnd(0) + 1.0) * dx[0]; + auto tile_lo_y = prob_lo[1] + (bx.smallEnd(1) + 0.0) * dx[1]; + auto tile_hi_y = prob_lo[1] + (bx.bigEnd(1) + 1.0) * dx[1]; + + // Only check polygon if not empty + if (!polygon_is_empty) { + // Check for overlap with polygon bounding box + if (tile_hi_x < bbox_lo[0] || tile_lo_x > bbox_hi[0] || + tile_hi_y < bbox_lo[1] || tile_lo_y > bbox_hi[1]) { + if (verbose) { + amrex::Print() + << m_key << " - Level " << level + << " (vertical_distance=" << vertical_distance + << "): Skipping tile: [" << tile_lo_x << "," + << tile_hi_x << "] x [" << tile_lo_y << "," << tile_hi_y + << "] " + << "does not overlap polygon bbox [" << bbox_lo[0] + << "," << bbox_hi[0] << "] x [" << bbox_lo[1] << "," + << bbox_hi[1] << "]\n"; + } + continue; // Skip this tile, no overlap + } + } + + amrex::ParallelFor( + bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + // 1. Early exit if already tagged + if (tag_arrs[nbx](i, j, k) == amrex::TagBox::SET) { + return; + } + + 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 amrex::Real terrainHt = interp::bilinear( + xterrain_ptr, xterrain_ptr + xterrain_size, yterrain_ptr, + yterrain_ptr + yterrain_size, zterrain_ptr, coord[0], + coord[1]); + + const auto cellHt = coord[2] - terrainHt; + + const amr_wind::polygon_utils::Polygon::Point center_point{ + coord[0], coord[1]}; + + // 2. Check vertical distance + if ((cellHt < -buffer_z_lo) || + (cellHt > std::max(vertical_distance, buffer_z_hi))) { + return; + } + + // 3. Check tagging box + if (!tagging_box.contains(coord)) { + return; + } + + // 4. Check polygon bounding box + if (!polygon_is_empty) { + if (!amr_wind::polygon_utils::Polygon:: + poly_bounding_box_contains( + center_point, bbox_lo, bbox_hi)) { + return; + } + } + + // 5. Point-in-polygon + bool in_poly = false; + if (!polygon_is_empty) { + in_poly = + amr_wind::polygon_utils::Polygon::is_point_in_polygon( + p_poly_points, p_ring_offsets, n_rings, + static_cast(n_points), center_point); + + } else { + in_poly = true; + } + + if (in_poly) { + tag_arrs[nbx](i, j, k) = amrex::TagBox::SET; + } + }); + } + + amrex::Gpu::streamSynchronize(); +} + +amrex::Vector TerrainRefinement::get_active_levels() const +{ + const int max_level = m_sim.repo().mesh().maxLevel(); + amrex::Vector levels; + int min_lev = std::max(0, m_set_level >= 0 ? m_set_level : m_min_level); + int max_lev = + std::min(max_level, m_set_level >= 0 ? m_set_level : m_max_level); + + for (int lev = min_lev; lev <= max_lev; ++lev) { + levels.push_back(lev); + } + return levels; +} + +} // namespace amr_wind \ No newline at end of file diff --git a/docs/sphinx/user/inputs_tagging.rst b/docs/sphinx/user/inputs_tagging.rst index 6dcaa914a9..f4bf06005a 100644 --- a/docs/sphinx/user/inputs_tagging.rst +++ b/docs/sphinx/user/inputs_tagging.rst @@ -41,6 +41,7 @@ Each section must contain the keyword ``type`` that is one of the refinement typ ``FieldRefinement`` Refinement based on error metric for field or its gradient ``OversetRefinement`` Refinement around fringe/field interface ``GeometryRefinement`` Refinement using geometric shapes +``TerrainRefinement`` Refinement using terrain fields and polygon regions ``QCriterionRefinement`` Refinement using Q-Criterion ``VorticityMagRefinement`` Refinement using vorticity ========================== =================================================================== @@ -130,7 +131,7 @@ block, and 2. ``cylinder`` -- refines the region inside a cylindrical block. Names of the input subsections that define specific geometries for refinement. -.. input_param:: tagging.GeoemtryRefinement.level +.. input_param:: tagging.GeometryRefinement.level **type:** Integer, optional, default: -1 @@ -212,6 +213,121 @@ The axis and the extents along the axis are defined by two position vectors optional ``inner_radius`` can be specified to restrict tagging to an annulus between the inner and outer radii. +Refinement using terrain and polygons +````````````````````````````````````` + +This section controls refinement using terrain fields in the domain along +with user-specified polygon regions. + +.. input_param:: tagging.TerrainRefinement.level + + **type:** Integer, optional, default: -1 + + If ``level`` is provided and is greater than or equal to 0, then the + refinement based on this section is only performed at that specific level. + +.. input_param:: tagging.TerrainRefinement.min_level + + **type:** Integer, optional, default: 0 + + If ``level`` is not specified, this option specifies the minimum level + where this refinement is active. + +.. input_param:: tagging.TerrainRefinement.max_level + + **type:** Integer, optional, default: ``mesh.maxLevel()`` + + If ``level`` is not specified, this option specifies the maximum level + where this refinement is active. + +.. input_param:: tagging.TerrainRefinement.vertical_distance + + **type:** Real, required + + Distance (in z) above the terrain to refine. Tagging is added between + the terrain and the specified height above it, and it is also + bound laterally by the polygon parameters listed below. + +.. input_param:: tagging.TerrainRefinement.grid_buffer_ratio_lo + + **type:** List of reals, optional, default = 0.0 + + Ratio of the subgrid (tile) height to use as a buffer **below** the terrain surface. + This buffer ensures that cells below the nominal terrain tagging region are also tagged, + which can help with smoother transitions and prevent under-refinement near the terrain interface. + The buffer size is computed as ``grid_buffer_ratio_lo * box_height`` for each tile, + where ``box_height`` is the physical height of the tile in the z-direction. + **Making this ratio larger will ensure that more cells below the terrain are tagged for refinement.** + **If multiple values are provided, they will be applied to different levels of the mesh.** + +.. input_param:: tagging.TerrainRefinement.grid_buffer_ratio_hi + + **type:** List of reals, optional, default = 0.0 + + Similar to ``grid_buffer_ratio_lo``, but this parameter specifies the ratio of the subgrid (tile) + height to use as a buffer **above** the terrain surface. + +.. input_param:: tagging.TerrainRefinement.poly_exterior + + **type:** List of reals, optional + + Coordinates (x and y) of the polygon exterior. Within this polygonal + region, the refinement is applied. The coordinates are input as pairs of x and y locations, + so there must be an even number of entries for this argument. + +.. input_param:: tagging.TerrainRefinement.poly_num_holes + + **type:** Integer, optional, default = 0 + + The number of holes within the polygonal region. This parameter allows the user to carve out + sections within the polygon exterior where refinement is not desired. Each hole is defined + in the same manner as the polygon exterior with coordinates constructing a polygon boundary. + +.. input_param:: tagging.TerrainRefinement.poly_hole_n + + **type:** List of reals, optional + + Coordinates (x and y) of the polygon hole boundary, where ``n`` is the 0-based index of the + hole being specified (e.g., ``poly_hole_0`` would be used to define the first polygon hole). + The coordinates are input as pairs of x and y locations, so there must be an even number + of entries for this argument. + +.. input_param:: tagging.TerrainRefinement.box_lo + + **type:** Vector, optional + + List of the low corner values for a bounding box where the tagging + will be active. By default the bounding box will span the entire domain. + +.. input_param:: tagging.TerrainRefinement.box_hi + + **type:** Vector, optional + + List of the high corner values for a bounding box where the tagging + will be active. By default the bounding box will span the entire domain. + +.. input_param:: tagging.TerrainRefinement.verbose + + **type:** List of ints, optional + + If provided, this parameter controls the verbosity of the tagging output on a per-level basis. + Each entry in the list corresponds to a mesh level, where a nonzero value enables verbose output for that level. + If fewer values are provided than the number of levels, the last value is used for all higher levels. + For example, if ``verbose = 1 0 0``, then only level 0 will print tagging information, + while levels 1 and above will not print any information. + +Example:: + + tagging.terr1.type = TerrainRefinement + tagging.terr1.vertical_distance = 200 + tagging.terr1.max_level = 1 # Refine only levels 1 and below + tagging.terr1.poly_exterior = 10 10 10 20 20 20 20 10 + tagging.terr1.poly_num_holes = 1 + tagging.terr1.poly_hole_0 = 5 5 5 10 10 10 10 5 + tagging.terr1.grid_buffer_ratio_lo = 0.0 0.75 + tagging.terr1.grid_buffer_ratio_hi = 0.0 0.75 + tagging.terr1.verbose = 0 1 + Refinement using Q-Criterion ````````````````````````````````````` diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 1aa09bd243..7507cadf16 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -274,6 +274,7 @@ add_test_re(rankine) add_test_re(rankine-sym) add_test_re(terrain_box) add_test_re(terrain_box_amr) +add_test_re(terrain_refinement_polygon_amr) add_test_re(forest_drag) add_test_re(box_refinement) add_test_re(cylinder_refinement) diff --git a/test/test_files/terrain_refinement_polygon_amr/terrain.amrwind b/test/test_files/terrain_refinement_polygon_amr/terrain.amrwind new file mode 100644 index 0000000000..9a12085262 --- /dev/null +++ b/test/test_files/terrain_refinement_polygon_amr/terrain.amrwind @@ -0,0 +1,290 @@ +16 +16 +300.0 +326.667 +353.333 +380.0 +406.667 +433.333 +460.0 +486.667 +513.333 +540.0 +566.667 +593.333 +620.0 +646.667 +673.333 +700.0 +300.0 +326.667 +353.333 +380.0 +406.667 +433.333 +460.0 +486.667 +513.333 +540.0 +566.667 +593.333 +620.0 +646.667 +673.333 +700.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +200.0 +200.0 +200.0 +200.0 +200.0 +200.0 +200.0 +200.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +200.0 +200.0 +200.0 +200.0 +200.0 +200.0 +200.0 +200.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +200.0 +200.0 +200.0 +200.0 +200.0 +200.0 +200.0 +200.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +200.0 +200.0 +200.0 +200.0 +200.0 +200.0 +200.0 +200.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +200.0 +200.0 +200.0 +200.0 +200.0 +200.0 +200.0 +200.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +200.0 +200.0 +200.0 +200.0 +200.0 +200.0 +200.0 +200.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +200.0 +200.0 +200.0 +200.0 +200.0 +200.0 +200.0 +200.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +200.0 +200.0 +200.0 +200.0 +200.0 +200.0 +200.0 +200.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 \ No newline at end of file diff --git a/test/test_files/terrain_refinement_polygon_amr/terrain_refinement_polygon_amr.inp b/test/test_files/terrain_refinement_polygon_amr/terrain_refinement_polygon_amr.inp new file mode 100644 index 0000000000..19f026008d --- /dev/null +++ b/test/test_files/terrain_refinement_polygon_amr/terrain_refinement_polygon_amr.inp @@ -0,0 +1,86 @@ +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# SIMULATION STOP # +#.......................................# +time.stop_time = -100.0 # Max (simulated) time to evolve +time.max_step = 0 # Max number of time steps + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# TIME STEP COMPUTATION # +#.......................................# +time.fixed_dt = 1 # Use this constant dt if > 0 +time.cfl = 0.9 # CFL factor + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# INPUT AND OUTPUT # +#.......................................# +time.plot_interval = 1000000 # Steps between plot files +time.checkpoint_interval = -1 # Steps between checkpoint files + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# PHYSICS # +#.......................................# +ConstValue.density.value = 1.0 +ConstValue.velocity.value = 1.0 0.0 0.0 + +incflo.use_godunov = 1 +incflo.diffusion_type = 2 +incflo.do_initial_proj = 0 +incflo.initial_iterations = 0 +transport.viscosity = 1.0e-15 +transport.laminar_prandtl = 0.7 +transport.turbulent_prandtl = 0.3333 +turbulence.model = Laminar +io.int_outputs = terrain_blank terrain_drag +incflo.physics = FreeStream TerrainDrag +ICNS.source_terms = DragForcing +DragForcing.is_laminar = 1 +TerrainDrag.terrain_file = "terrain.amrwind" +amr.n_cell = 56 56 52 # Grid cells at coarsest AMRlevel +geometry.prob_lo = 0 0 0 +geometry.prob_hi = 1024 1024 1024 +geometry.is_periodic = 0 0 0 # Periodicity x y z (0/1) + +# Boundary conditions +xlo.type = "mass_inflow" +xlo.density = 1.0 +xlo.velocity = 1.0 0.0 0.0 +xhi.type = "pressure_outflow" +ylo.type = "slip_wall" +yhi.type = "slip_wall" +zlo.type = "slip_wall" +zhi.type = "slip_wall" + +incflo.verbose = 0 # incflo_level + +amr.blocking_factor_z = 4 +amr.max_level = 3 +amr.grid_eff = 1 +amr.n_error_buf = 1 + +tagging.labels = terrain_lvl0 terrain_lvl1 terrain_lvl1_whole_terrain terrain_lvl2 + +tagging.terrain_lvl0.type = TerrainRefinement +tagging.terrain_lvl0.verbose = 1 +tagging.terrain_lvl0.level = 0 +tagging.terrain_lvl0.vertical_distance = 200 +tagging.terrain_lvl0.poly_exterior = 898 461 892 422 883 384 870 347 853 311 833 278 809 246 783 217 754 191 722 167 689 147 653 130 616 117 578 108 539 102 500 100 461 102 422 108 384 117 347 130 311 147 278 167 246 191 217 217 191 246 167 278 147 311 130 347 117 384 108 422 102 461 100 500 102 539 108 578 117 616 130 653 147 689 167 722 191 754 217 783 246 809 278 833 311 853 347 870 384 883 422 892 461 898 500 900 539 898 578 892 616 883 653 870 689 853 722 833 754 809 783 783 809 754 833 722 853 689 870 653 883 616 892 578 898 539 900 500 898 461 +tagging.terrain_lvl0.poly_num_holes = 1 +tagging.terrain_lvl0.poly_hole_0 = 699 520 696 539 691 558 685 577 676 594 666 611 655 627 641 641 627 655 611 666 594 676 577 685 558 691 539 696 520 699 500 700 480 699 461 696 442 691 423 685 406 676 389 666 373 655 359 641 345 627 334 611 324 594 315 577 309 558 304 539 301 520 300 500 301 480 304 461 309 442 315 423 324 406 334 389 345 373 359 359 373 345 389 334 406 324 423 315 442 309 461 304 480 301 500 300 520 301 539 304 558 309 577 315 594 324 611 334 627 345 641 359 655 373 666 389 676 406 685 423 691 442 696 461 699 480 700 500 699 520 + +tagging.terrain_lvl1.type = TerrainRefinement +tagging.terrain_lvl1.verbose = 0 +tagging.terrain_lvl1.level = 1 +tagging.terrain_lvl1.vertical_distance = 100 +tagging.terrain_lvl1.poly_exterior = 700 500 699 480 696 461 691 442 685 423 676 406 666 389 655 373 641 359 627 345 611 334 594 324 577 315 558 309 539 304 520 301 500 300 480 301 461 304 442 309 423 315 406 324 389 334 373 345 359 359 345 373 334 389 324 406 315 423 309 442 304 461 301 480 300 500 301 520 304 539 309 558 315 577 324 594 334 611 345 627 359 641 373 655 389 666 406 676 423 685 442 691 461 696 480 699 500 700 520 699 539 696 558 691 577 685 594 676 611 666 627 655 641 641 655 627 666 611 676 594 685 577 691 558 696 539 699 520 700 500 + +# This will apply to the whole terrain +tagging.terrain_lvl1_whole_terrain.type = TerrainRefinement +tagging.terrain_lvl1_whole_terrain.verbose = 0 +tagging.terrain_lvl1_whole_terrain.level = 1 +tagging.terrain_lvl1_whole_terrain.vertical_distance = 50 + +tagging.terrain_lvl2.type = TerrainRefinement +tagging.terrain_lvl2.verbose = 1 +tagging.terrain_lvl2.level = 2 +tagging.terrain_lvl2.vertical_distance = 10 +tagging.terrain_lvl2.poly_exterior = 699 520 696 539 691 558 685 577 676 594 666 611 655 627 641 641 627 655 611 666 594 676 577 685 558 691 539 696 520 699 500 700 480 699 461 696 442 691 423 685 406 676 389 666 373 655 359 641 345 627 334 611 324 594 315 577 309 558 304 539 301 520 300 500 301 480 304 461 309 442 315 423 324 406 334 389 345 373 359 359 373 345 389 334 406 324 423 315 442 309 461 304 480 301 500 300 520 301 539 304 558 309 577 315 594 324 611 334 627 345 641 359 655 373 666 389 676 406 685 423 691 442 696 461 699 480 700 500 699 520 diff --git a/unit_tests/utilities/CMakeLists.txt b/unit_tests/utilities/CMakeLists.txt index d90c23b575..43da352376 100644 --- a/unit_tests/utilities/CMakeLists.txt +++ b/unit_tests/utilities/CMakeLists.txt @@ -14,6 +14,7 @@ target_sources(${amr_wind_unit_test_exe_name} PRIVATE test_field_norms.cpp test_tensor_ops.cpp test_post_processing_time.cpp + test_polygon.cpp test_time_averaging.cpp ) diff --git a/unit_tests/utilities/test_polygon.cpp b/unit_tests/utilities/test_polygon.cpp new file mode 100644 index 0000000000..b5f58cdd04 --- /dev/null +++ b/unit_tests/utilities/test_polygon.cpp @@ -0,0 +1,232 @@ +#include "amr-wind/utilities/tagging/Polygon.H" +#include +#include + +using namespace amr_wind::polygon_utils; + +TEST(Polygon, SimpleSquare) +{ + Polygon poly; + poly.add_outer_vertex({0.0, 0.0}); + poly.add_outer_vertex({0.0, 10.0}); + poly.add_outer_vertex({10.0, 10.0}); + poly.add_outer_vertex({10.0, 0.0}); + poly.add_outer_vertex({0.0, 0.0}); + poly.compute_bounding_box(); + + const auto* pts = poly.ring_ptr(0); + int n_outer = poly.ring_size(0); + const auto* all_pts = poly.points().data(); + const auto* ring_offsets = poly.ring_offsets().data(); + int n_rings = poly.num_rings(); + int n_points = poly.num_points(); + + std::cout << "SimpleSquare: Testing point (5,5) (inside)\n"; + EXPECT_TRUE(poly.contains({5.0, 5.0})); + EXPECT_TRUE(Polygon::is_point_in_ring(pts, n_outer, {5.0, 5.0})); + EXPECT_TRUE( + Polygon::is_point_in_polygon( + all_pts, ring_offsets, n_rings, n_points, {5.0, 5.0})); + + std::cout << "SimpleSquare: Testing point (0,5) (on edge)\n"; + EXPECT_FALSE(poly.contains({0.0, 5.0})); + EXPECT_FALSE(Polygon::is_point_in_ring(pts, n_outer, {0.0, 5.0})); + EXPECT_FALSE( + Polygon::is_point_in_polygon( + all_pts, ring_offsets, n_rings, n_points, {0.0, 5.0})); + + std::cout << "SimpleSquare: Testing point (-1,5) (outside)\n"; + EXPECT_FALSE(poly.contains({-1.0, 5.0})); + EXPECT_FALSE(Polygon::is_point_in_ring(pts, n_outer, {-1.0, 5.0})); + EXPECT_FALSE( + Polygon::is_point_in_polygon( + all_pts, ring_offsets, n_rings, n_points, {-1.0, 5.0})); +} + +TEST(Polygon, PolygonWithHole) +{ + Polygon poly; + // Outer square + poly.add_outer_vertex({0.0, 0.0}); + poly.add_outer_vertex({0.0, 10.0}); + poly.add_outer_vertex({10.0, 10.0}); + poly.add_outer_vertex({10.0, 0.0}); + poly.add_outer_vertex({0.0, 0.0}); + poly.start_hole(); + // Hole: smaller square + poly.add_vertex({3.0, 3.0}); + poly.add_vertex({3.0, 7.0}); + poly.add_vertex({7.0, 7.0}); + poly.add_vertex({7.0, 3.0}); + poly.add_vertex({3.0, 3.0}); + poly.compute_bounding_box(); + + const auto* all_pts = poly.points().data(); + const auto* ring_offsets = poly.ring_offsets().data(); + int n_rings = poly.num_rings(); + int n_points = poly.num_points(); + + std::cout << "PolygonWithHole: Testing point (2,2) (inside outer, outside " + "hole)\n"; + EXPECT_TRUE(poly.contains({2.0, 2.0})); + EXPECT_TRUE( + Polygon::is_point_in_polygon( + all_pts, ring_offsets, n_rings, n_points, {2.0, 2.0})); + + std::cout << "PolygonWithHole: Testing point (5,5) (inside hole)\n"; + EXPECT_FALSE(poly.contains({5.0, 5.0})); + EXPECT_FALSE( + Polygon::is_point_in_polygon( + all_pts, ring_offsets, n_rings, n_points, {5.0, 5.0})); + + std::cout << "PolygonWithHole: Testing point (10,5) (on exterior edge)\n"; + EXPECT_FALSE(poly.contains({10.0, 5.0})); + EXPECT_FALSE( + Polygon::is_point_in_polygon( + all_pts, ring_offsets, n_rings, n_points, {10.0, 5.0})); + + std::cout << "PolygonWithHole: Testing point (3,5) (on hole edge)\n"; + EXPECT_FALSE(poly.contains({3.0, 5.0})); + EXPECT_FALSE( + Polygon::is_point_in_polygon( + all_pts, ring_offsets, n_rings, n_points, {3.0, 5.0})); + + std::cout << "PolygonWithHole: Testing point (-1,5) (outside)\n"; + EXPECT_FALSE(poly.contains({-1.0, 5.0})); + EXPECT_FALSE( + Polygon::is_point_in_polygon( + all_pts, ring_offsets, n_rings, n_points, {-1.0, 5.0})); +} + +TEST(Polygon, ComplexPolygonSelfIntersecting) +{ + Polygon poly; + // Bowtie (self-intersecting) + poly.add_outer_vertex({0.0, 0.0}); + poly.add_outer_vertex({5.0, 10.0}); + poly.add_outer_vertex({10.0, 0.0}); + poly.add_outer_vertex({0.0, 10.0}); + poly.add_outer_vertex({10.0, 10.0}); + poly.add_outer_vertex({0.0, 0.0}); + poly.compute_bounding_box(); + + const auto* all_pts = poly.points().data(); + const auto* ring_offsets = poly.ring_offsets().data(); + int n_rings = poly.num_rings(); + int n_points = poly.num_points(); + + std::cout << "ComplexPolygonSelfIntersecting: Testing point (5,5) (center, " + "ambiguous)\n"; + EXPECT_FALSE(poly.contains({5.0, 5.0})); + EXPECT_FALSE( + Polygon::is_point_in_polygon( + all_pts, ring_offsets, n_rings, n_points, {5.0, 5.0})); +} + +TEST(Polygon, DegenerateCases) +{ + Polygon poly; + // Single point + poly.add_outer_vertex({1.0, 1.0}); + poly.compute_bounding_box(); + const auto* all_pts = poly.points().data(); + const auto* ring_offsets = poly.ring_offsets().data(); + int n_rings = poly.num_rings(); + int n_points = poly.num_points(); + + std::cout << "DegenerateCases: Testing single point polygon at (1,1)\n"; + EXPECT_FALSE(poly.contains({1.0, 1.0})); + EXPECT_FALSE( + Polygon::is_point_in_polygon( + all_pts, ring_offsets, n_rings, n_points, {1.0, 1.0})); + + // Line + Polygon poly2; + poly2.add_outer_vertex({0.0, 0.0}); + poly2.add_outer_vertex({1.0, 1.0}); + poly2.compute_bounding_box(); + + const auto* all_pts2 = poly2.points().data(); + const auto* ring_offsets2 = poly2.ring_offsets().data(); + int n_rings2 = poly2.num_rings(); + int n_points2 = poly2.num_points(); + + std::cout << "DegenerateCases: Testing line polygon from (0,0) to (1,1)\n"; + EXPECT_FALSE(poly2.contains({0.5, 0.5})); + EXPECT_FALSE( + Polygon::is_point_in_polygon( + all_pts2, ring_offsets2, n_rings2, n_points2, {0.5, 0.5})); +} + +TEST(Polygon, BoundingBox) +{ + Polygon poly; + poly.add_outer_vertex({-2.0, 3.0}); + poly.add_outer_vertex({4.0, 5.0}); + poly.add_outer_vertex({1.0, -1.0}); + poly.add_outer_vertex({-2.0, 3.0}); + poly.compute_bounding_box(); + + Polygon::Point bbox_lo = {0.0, 0.0}; + Polygon::Point bbox_hi = {0.0, 0.0}; + poly.get_bounding_box(bbox_lo, bbox_hi); + + std::cout << "BoundingBox: Testing point (0,0) (inside bbox)\n"; + EXPECT_TRUE(poly.bounding_box_contains({0.0, 0.0})); + EXPECT_TRUE( + Polygon::poly_bounding_box_contains({0.0, 0.0}, bbox_lo, bbox_hi)); + std::cout << "BoundingBox: Testing point (10,10) (outside bbox)\n"; + EXPECT_FALSE(poly.bounding_box_contains({10.0, 10.0})); + EXPECT_FALSE( + Polygon::poly_bounding_box_contains({10.0, 10.0}, bbox_lo, bbox_hi)); +} + +TEST(Polygon, MultipleHoles) +{ + Polygon poly; + // Outer square + poly.add_outer_vertex({0.0, 0.0}); + poly.add_outer_vertex({0.0, 10.0}); + poly.add_outer_vertex({10.0, 10.0}); + poly.add_outer_vertex({10.0, 0.0}); + poly.add_outer_vertex({0.0, 0.0}); + poly.start_hole(); + // Hole 1 + poly.add_vertex({2.0, 2.0}); + poly.add_vertex({2.0, 4.0}); + poly.add_vertex({4.0, 4.0}); + poly.add_vertex({4.0, 2.0}); + poly.add_vertex({2.0, 2.0}); + poly.start_hole(); + // Hole 2 + poly.add_vertex({6.0, 6.0}); + poly.add_vertex({6.0, 8.0}); + poly.add_vertex({8.0, 8.0}); + poly.add_vertex({8.0, 6.0}); + poly.add_vertex({6.0, 6.0}); + poly.compute_bounding_box(); + + const auto* all_pts = poly.points().data(); + const auto* ring_offsets = poly.ring_offsets().data(); + int n_rings = poly.num_rings(); + int n_points = poly.num_points(); + + std::cout + << "MultipleHoles: Testing point (5,5) (inside outer, outside holes)\n"; + EXPECT_TRUE(poly.contains({5.0, 5.0})); + EXPECT_TRUE( + Polygon::is_point_in_polygon( + all_pts, ring_offsets, n_rings, n_points, {5.0, 5.0})); + + std::cout << "MultipleHoles: Testing point (3,3) (inside hole 1)\n"; + EXPECT_FALSE(poly.contains({3.0, 3.0})); + EXPECT_FALSE( + Polygon::is_point_in_polygon( + all_pts, ring_offsets, n_rings, n_points, {3.0, 3.0})); + + std::cout << "MultipleHoles: Testing point (7,7) (inside hole 2)\n"; + EXPECT_FALSE(poly.contains({7.0, 7.0})); + EXPECT_FALSE( + Polygon::is_point_in_polygon( + all_pts, ring_offsets, n_rings, n_points, {7.0, 7.0})); +} \ No newline at end of file