Added terrain refinement option in tagging #1633
Conversation
|
With a polygon as a boundary, one can have rather good control!
Missing
|
|
The implementation is ugly, and my C++ skills are very rusty. One can do a lot better (polygon class, DRY principle), but here is a complete implementation with inner rings.
|
Note to maintainersI know this is very rough and not according to your standards. I haven't programmed in C++ in 15 years 🫨. It also solves some annoying problems referenced in #1629 (like not being able to capture a flat surface (no blank)). For the user, it is also very intuitive. |
hgopalan
requested review from
hgopalan,
mchurchf,
marchdf,
moprak-nrel and
mbkuhn
Update with real example
Python code import numpy as np
from shapely.geometry import Polygon, Point
import IPython.display
def get_amr_wind_polygon(geom: Polygon, precision: int = 0) -> tuple[str, list[str]]:
assert isinstance(geom, Polygon)
assert geom.geom_type == "Polygon"
outer_coords = np.asarray(geom.exterior.coords)[:, :2].round(precision)
if precision == 0:
outer_coords = outer_coords.astype("int")
poly_outer = " ".join(map(str, outer_coords.flatten()))
poly_inners = []
if len(geom.interiors):
for ring_geom in geom.interiors:
ring_coords = np.asarray(ring_geom.coords)[:, :2].round(precision)
if precision == 0:
ring_coords = ring_coords.astype("int")
poly_inners.append(
" ".join(map(str, ring_coords.flatten()))
)
return poly_outer, poly_inners
level = 1
geom = Point(500, 500).buffer(400)
if 1:
geom = geom.difference(Point(500, 500).buffer(200))
IPython.display.display(geom)
poly_outer, poly_inners = get_amr_wind_polygon(geom)
print(f"tagging.terrain_lvl{level}.type = TerrainRefinement")
print(f"tagging.terrain_lvl{level}.vertical_distance = 100")
print(f"tagging.terrain_lvl{level}.level = {level}")
print(f"tagging.terrain_lvl{level}.poly_exterior = {poly_outer}")
if poly_inners:
print(f"tagging.terrain_lvl{level}.poly_num_holes = {len(poly_inners)}")
for i, l in enumerate(poly_inners):
print(f"tagging.terrain_lvl{level}.poly_hole_{i} = {l}")
level = 2
geom = Point(500, 500).buffer(200)
IPython.display.display(geom)
poly_outer, poly_inners = get_amr_wind_polygon(geom)
print(f"tagging.terrain_lvl{level}.type = TerrainRefinement")
print(f"tagging.terrain_lvl{level}.vertical_distance = 40")
print(f"tagging.terrain_lvl{level}.level = {level}")
print(f"tagging.terrain_lvl{level}.poly_exterior = {poly_outer}")
if poly_inners:
print(f"tagging.terrain_lvl{level}.poly_num_holes = {len(poly_inners)}")
for i, l in enumerate(poly_inners):
print(f"tagging.terrain_lvl{level}.poly_hole_{i} = {l}") |
|
Basically, one needs to assert after mesh creation the refinement level at the location of some points. I am so sorry that I cannot do that. This is probably trivial... |
There was a problem hiding this comment.
This is the same file as used for the terrain drag test (taken from test/test_files/terrain_box.).
Can you provide more details on why is this assertion necessary? |
There was a problem hiding this comment.
I know this is ugly as ***. This was a quick proof of concept, taken from internet. I had used this mehod before in Openfoam.
As we probably won't have super complex polygons, I think this is fast enough ...
AMReX probably already has a class for this, but I haven't found it.
There was a problem hiding this comment.
Basically, this should just generate a mesh. I have cranked the grid_eff to 1 for testing purpose. But I have to say that for my test case in complex terrain (see #1630), I kinda feel that this is needed. Also taken from test/test_files/terrain_box.
To create a unit test I mean. The generated mesh should have 2 refinement levels. I don't know how to find the nearest cell of a point in a mesh, test for the maximum refinement level, or actually write unit testing in cpp. This is why I said that for a maintainer, this is probably rather trivial. |
|
It fails on GPU. I would like to change this to draft, as I don't want to trigger costly CI. |
![github-advanced-security[bot]](https://avatars.githubusercontent.com/in/57789?s=60&v=4){kind=small}
|
This is my attempt to fix the GPU compilation. Note that I have not tested it and it may not even compile. But it gives you the idea of how to port it to GPU. Please feel free to use it in whatever way you prefer. |
Thank you so much! I have to admit, storing a list of lists for the inner rings was very stupid. |
Polygon testsI had CoPilot generate the test-suite in python, to better see the process. # Polygon test visualization and logic using matplotlib and shapely
import matplotlib.pyplot as plt
from shapely.geometry import Polygon, Point, LinearRing
from shapely.errors import TopologicalError
def bounding_box_contains(poly, pt):
minx, miny, maxx, maxy = poly.bounds
x, y = pt
return (minx <= x <= maxx) and (miny <= y <= maxy)
def plot_polygon(ax, poly, holes=None, color="lightblue", edgecolor="blue", alpha=0.3):
if holes is None:
holes = []
x, y = poly.exterior.xy
ax.fill(x, y, color=color, alpha=alpha, edgecolor=edgecolor, linewidth=2)
for hole in holes:
hx, hy = hole.xy
ax.fill(hx, hy, color="white", alpha=1, edgecolor="red", linewidth=2)
def plot_points(ax, points, results, labels=None):
for i, (pt, res) in enumerate(zip(points, results)):
color = "green" if res else "red"
ax.plot(pt[0], pt[1], "o", color=color, markersize=10)
if labels:
ax.text(pt[0], pt[1], f"{labels[i]}", fontsize=9, ha="left")
def run_test(title, poly, test_points, expected, holes=None, show_bbox=False):
fig, ax = plt.subplots()
ax.set_title(title)
plot_polygon(ax, poly, holes)
results = [poly.contains(Point(p)) for p in test_points]
plot_points(ax, test_points, results, labels=[str(i) for i in range(len(test_points))])
for i, (pt, exp, res) in enumerate(zip(test_points, expected, results)):
print(f"Test point {pt}: expected {exp}, got {res}")
assert res == exp, f"Test failed for point {pt}: expected {exp}, got {res}"
if show_bbox:
minx, miny, maxx, maxy = poly.bounds
ax.plot([minx, maxx, maxx, minx, minx], [miny, miny, maxy, maxy, miny], "k--", lw=2)
ax.set_aspect("equal")
plt.show()
# 1. SimpleSquare
poly1 = Polygon([(0, 0), (0, 10), (10, 10), (10, 0), (0, 0)])
run_test("SimpleSquare", poly1, [(5, 5), (0, 5), (-1, 5)], [True, False, False])
# 2. PolygonWithHole
outer = [(0, 0), (0, 10), (10, 10), (10, 0), (0, 0)]
hole = [(3, 3), (3, 7), (7, 7), (7, 3), (3, 3)]
poly2 = Polygon(outer, [hole])
run_test(
"PolygonWithHole",
poly2,
[(2, 2), (5, 5), (10, 5), (3, 5), (-1, 5)],
[True, False, False, False, False],
holes=[LinearRing(hole)],
)
# 3. ComplexPolygonSelfIntersecting (bowtie)
bowtie = [(0, 0), (5, 10), (10, 0), (0, 10), (10, 10), (0, 0)]
try:
poly3 = Polygon(bowtie)
# Shapely may throw TopologicalError for self-intersecting polygons
run_test("ComplexPolygonSelfIntersecting", poly3, [(5, 5)], [False])
except TopologicalError as e:
print("Self-intersecting polygon could not be constructed:", e)
# 4. DegenerateCases: single point and line
try:
poly4 = Polygon([(1, 1)])
run_test("DegenerateCases: Single Point", poly4, [(1, 1)], [False])
except ValueError:
print("DegenerateCases: Single Point - Polygon not valid, contains always False")
try:
poly5 = Polygon([(0, 0), (1, 1)])
run_test("DegenerateCases: Line", poly5, [(0.5, 0.5)], [False])
except ValueError:
print("DegenerateCases: Line - Polygon not valid, contains always False")
# 5. BoundingBox
poly6 = Polygon([(-2,3), (4,5), (1,-1), (-2,3)])
test_points = [(0,0), (10,10)]
expected = [
bounding_box_contains(poly6, (0,0)),
bounding_box_contains(poly6, (10,10))
]
def run_bbox_test(title, poly, test_points, expected):
fig, ax = plt.subplots()
ax.set_title(title)
plot_polygon(ax, poly)
results = [bounding_box_contains(poly, p) for p in test_points]
plot_points(ax, test_points, results, labels=[str(i) for i in range(len(test_points))])
for i, (pt, exp, res) in enumerate(zip(test_points, expected, results)):
print(f"Test point {pt}: expected {exp}, got {res}")
assert res == exp, f"Test failed for point {pt}: expected {exp}, got {res}"
minx, miny, maxx, maxy = poly.bounds
ax.plot([minx, maxx, maxx, minx, minx], [miny, miny, maxy, maxy, miny], 'k--', lw=2)
ax.set_aspect('equal')
plt.show()
run_bbox_test("BoundingBox", poly6, test_points, expected)
# 6. MultipleHoles
outer = [(0, 0), (0, 10), (10, 10), (10, 0), (0, 0)]
hole1 = [(2, 2), (2, 4), (4, 4), (4, 2), (2, 2)]
hole2 = [(6, 6), (6, 8), (8, 8), (8, 6), (6, 6)]
poly7 = Polygon(outer, [hole1, hole2])
run_test(
"MultipleHoles",
poly7,
[(5, 5), (3, 3), (7, 7)],
[True, False, False],
holes=[LinearRing(hole1), LinearRing(hole2)],
)
|
|
Thank you for putting up the PR to close the feature request. |
Thanks. I have to say that CoPilot helped me a lot to get back into CPP programming. The OpenFoam library makes it almost too easy because every possible container type and parsing is implemented! |
Bugs
|
Results on my project's domain
I just dont know how to remove these. Hopefully, the happen far away enough from my ROI.
|
|
For those errors that you're seeing, what is the vertical distance specified for the levels that are having the problem? And what is the dz for those levels? |
|
|
||
| tagging.terrain_lvl1.type = TerrainRefinement | ||
| tagging.terrain_lvl1.vertical_distance = 200 | ||
| tagging.terrain_lvl1.level = 1 |
There was a problem hiding this comment.
I'm noticing that the first level in the tagging is 1. This argument specifies which level the tagging should be applied to, not which level will be created by the tagging. It's confusing. But the tagging of level = 0 will determine where the cells of level = 1 show up. This might be affecting the other case that you are having challenges with.
There was a problem hiding this comment.
Yes, I know it is confusing. In my mind, it is the final target grid level, not the tagging process. I actually wonder if I should round up the vertical distance by the current grid dz when tagging the cells. but somehow, I don't think this is what causes the problems. These cells are close tobthe terrain, and if you see the cells around, they are clearly refined all the way up.
Could it be the ngrow paremeter that I ignore?
There was a problem hiding this comment.
I will post tomorrow my inp file for the complex case. The refinements are not the same as in the example.
Thank you so much for the assistance!
Would you like me to run a case with the standard grid efficiency of 70%, to show the differences?
There was a problem hiding this comment.
I have changed the implementation to be more inline with the rest of the code.
Level now referes to the grid level the grid refinement should apply to.
There was a problem hiding this comment.
@mbkuhn , with my updates, this disappeared.
Changes:
- Use the same convention as in the rest of the code for level
- Round up the vertical distance to the nearest multiple of the current grid size
I will try to get rid of the annoying mis-refinements on some edges, but i am afraid there is just no way with AMReX meshing. I am posting pictures and setup here:
On top of my head, I recall running with 5 levels, the finest being 10m, and an aspect ratio of 4. I will need to check my terrain refinement level, but I think that the 2 big patches that fail are level 2. level 1 is 600m high, level 1 300m across the whole domain. 3 4 and 5 are just around the wind farm. |
|
I wonder if I should make sure that the rank actually contains the polygon, like it is done if the forest implementation. At creation time, the polygon bbox is computed. It is only in the tagging process that the bbox is lost (someone with skills could probably easily fix that, and add a bbox check) https://github.com/ews-ffarella/amr-wind/blob/main/amr-wind/physics/ForestDrag.cpp#L118 |
Level 0 is allowed, and will refine the base grid!
…amr-wind into 1-terrain-based-refinement
…amr-wind into 1-terrain-based-refinement
Summary
This PR introduces polygon-based terrain refinement for AMR tagging, enabling users to specify arbitrary 2D polygonal regions (with optional holes) for mesh refinement above terrain. The new feature at least compiles on GPU (thanks @WeiqunZhang) and supports multiple refinement levels and flexible region definitions via the input file.
Polygonal Terrain Refinement
The new
Polygonclass allows users to define arbitrary 2D polygons (with holes) for refinement regions.poly_exteriorfor the outer ringpoly_num_holesandpoly_hole_Nfor holesFlexible Tagging
The new
TerrainRefinementclass supports:box_lo,box_hi)Of course, using
Geoswould be a lot better, but this would be an overkill!See example picture:
The black area is the terrain mask
The advantages of this method are:
box_lo,box_hi, etc as this is the case withGeometryRefinement.Future idea(s)
terrain_height. In flat terrain, but with forest, one could just use the z-coordinate as a proxy.Pull request type
Please check the type of change introduced:
Checklist
The following is included:
unit_tests/utilities/test_polygon.cpptest/test_files/terrain_refinement_polygon_amr/terrain_refinement_polygon_amr.inpPolygon.HandTerrainRefinement.HThis PR was tested by running:
See CI summary on my repo.
Additional background
Often, one has a collection of interesting objects (measurements, turbines, planning area) for which high resolution is needed.
Refining these using geometric entities like boxes and cylinders is very messy and inaccurate (requires a lot of pre-processing effort to be sure not to refine crazy large area). This PR aims at providing an easy and flexible way of refining the domain.
Example levels: