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. |
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
Otherwise, problems after regridding or on restart
Major changes in this commitWhen running with AMR cases and regridding, I have noticed that one cannot rely on the TerrainDrag derived fields. On regridding, they are being re-initialized after Refinement operator() is called. This is not the expected behavior. This is almost a good thing, as it always bothered me to have to rely on these. So instead:
Other changes:
Resultsamr.n_cell = 128 128 112 # Grid cells at coarsest AMRlevel
amr.max_level = 4
amr.blocking_factor_x = 8
amr.blocking_factor_y = 8
amr.max_grid_size_x = 32
amr.max_grid_size_y = 32
amr.blocking_factor_z = 4
amr.max_grid_size_z = 24
amr.grid_eff = 0.95 # Could probably be relaxed now
# With high enough grid_buffer_ratio_lo in high mesh resolution area, no need to buffer
amr.n_error_buf = 0
amr.n_error_buf_z = 0
tagging.labels = terrain_level_0_2000m terrain_level_1_500m terrain_level_2_50m terrain_level_2_300m terrain_level_3_100m
tagging.terrain_level_0_2000m.type = TerrainRefinement
tagging.terrain_level_0_2000m.vertical_distance = 2000
tagging.terrain_level_0_2000m.max_level = 0
tagging.terrain_level_0_2000m.grid_buffer_ratio_lo = 0.0
tagging.terrain_level_0_2000m.grid_buffer_ratio_hi = 0.0
tagging.terrain_level_1_500m.type = TerrainRefinement
tagging.terrain_level_1_500m.vertical_distance = 400
tagging.terrain_level_1_500m.max_level = 1
tagging.terrain_level_1_500m.grid_buffer_ratio_lo = 0.0
tagging.terrain_level_1_500m.grid_buffer_ratio_hi = 0.0
tagging.terrain_level_2_50m.type = TerrainRefinement
tagging.terrain_level_2_50m.vertical_distance = 1
tagging.terrain_level_2_50m.max_level = 2
tagging.terrain_level_2_50m.grid_buffer_ratio_lo = 0 0 0.75
tagging.terrain_level_2_50m.grid_buffer_ratio_hi = 0 0 0.75
tagging.terrain_level_2_50m.poly_exterior = 3060 12235 ...
tagging.terrain_level_2_300m.type = TerrainRefinement
tagging.terrain_level_2_300m.vertical_distance = 200
tagging.terrain_level_2_300m.max_level = 2
tagging.terrain_level_2_300m.grid_buffer_ratio_lo = 0 0 0.75
tagging.terrain_level_2_300m.grid_buffer_ratio_hi = 0 0 0.75
tagging.terrain_level_2_300m.poly_exterior = 1871 4323 ...
tagging.terrain_level_3_100m.type = TerrainRefinement
tagging.terrain_level_3_100m.vertical_distance = 100
# tagging.terrain_level_3_100m.verbose = 0 0 1
tagging.terrain_level_3_100m.max_level = 3
tagging.terrain_level_3_100m.grid_buffer_ratio_lo = 0 0 0.75
tagging.terrain_level_3_100m.grid_buffer_ratio_hi = 0 0 0.75
tagging.terrain_level_3_100m.poly_exterior = 1871 4323 ...By increasing the Finer area Edges area From bellow (terrain refinement)
@hgopalan what do you think of moving the terrain height reading / data storing outside terrain trag, like i did? Also, is the singleton pattern bad? I guess resources are hogged for the whole computation time, which is not the case whrn you only read it in operator(). |
…proved consistency - Allow grid_buffer_ratio_lo and grid_buffer_ratio_hi to be specified as lists for per-level control - Support per-level verbose output via a list of ints, with safe fallback for missing levels - Align level selection logic (set_level, min_level, max_level) with GeometryRefinement for consistency - Add utility for formatting lists as strings - Improve documentation and add notes for future extensibility of TerrainElevation - Fix potential out-of-bounds access in verbose logic
Influence of regridding before first time-step, when restarting from a coarser checkpointThe refinement parameters are exactly the same. Without restart
With restart
I still need to understand why ... Still problems with restartSame inp as before, but starting from a coarser mesh.
Cranking up grid_buffer_ratio on all levels and heigh
|
| static std::shared_ptr<TerrainElevation>& get_instance(); | ||
| static void ensure_loaded(const std::string& filename); | ||
|
|
||
| // Accessors (const and non-const) |
There was a problem hiding this comment.
This all feels like overkill. How about we just use a struct with public members? Also none of these should be DEVICE since we should not be accessing a Vector on device.
There was a problem hiding this comment.
You are surely right. Sorry, but my CPP skills are not the best, and i didn't have time yet to fully familiarize myself with your codebase.
| } | ||
|
|
||
| // Get the buffer ratios for this level | ||
| auto buffer_lo = m_grid_buffer_ratio_lo.empty() |
There was a problem hiding this comment.
A bunch of these autos in this file can be const auto
|
This PR is stale because it has been open 30 days with no activity. Remove stale label or comment or this will be closed in 7 days. |
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: