add ConservativeRegridding extension #2400
Conversation
juliasloan25
force-pushed
the
js/kp/conservative
branch
2 times, most recently
from
0809c2a to
2579639
Compare
|
Thanks to @imreddyTeja for helping with debuging the zero area polygons. The zero area polygons are line segments on a lonlat grid. If you convert to Cartesian coordinates, and plot them on a XYZ grid, then the vertices of the elements have positive area. In particular, they correspond to polygons at the north and south poles. See the plot below where I plotted the vertices of one of the zero area polygons as yellow nodes. |
juliasloan25
changed the title
juliasloan25
force-pushed
the
js/kp/conservative
branch
from
f6e640a to
6c4e945
Compare
test/ClimaCoreConservativeRegriddingExt/test_conservative_regridding.jl
Outdated
Show resolved
Hide resolved
juliasloan25
force-pushed
the
js/kp/conservative
branch
from
a34c123 to
cb97188
Compare
juliasloan25
force-pushed
the
js/kp/conservative
branch
2 times, most recently
from
eeb08ee to
52748ca
Compare
juliasloan25
force-pushed
the
js/kp/conservative
branch
3 times, most recently
from
ce2b29b to
e998d45
Compare
juliasloan25
force-pushed
the
js/kp/conservative
branch
from
235a37f to
f5571ef
Compare
| # Check for zero area polygons (all latitude or longitude values are the same) | ||
| for polygon in vertices | ||
| if allequal(first.(polygon)) || allequal(last.(polygon)) | ||
| @error "Zero area polygon found in vertices" polygon | ||
| end | ||
| end |
There was a problem hiding this comment.
I think the error message can be better here. They are only zero area polygons, because the polygons are on a longitude-latitude grid. If you view the polygons on a sphere, then they are not zero area. It might also be helpful to point out that the current implementation of ConservativeRegridding does not work with these zero area polygons in the docstring or error message.
| CartesianIndex(i, j, 1, 1, e) # f and v are 1 for SpectralElementSpace2D | ||
| for e in 1:Nh | ||
| for (i, j) in [(1, 1), (1, Nq), (Nq, Nq), (Nq, 1), (1, 1)] | ||
| ] # repeat the first coordinate pair at the end |
There was a problem hiding this comment.
Does this only work with a SpectralElementSpace2D? If it is, can you add a type annotation for that?
| function Remapping.get_value_per_element!( | ||
| value_per_element, | ||
| field, | ||
| ones_field, | ||
| ) |
There was a problem hiding this comment.
Is this function meant to be used by ClimaCore users or is it internal? If it is the latter, can you rename this to _get_value_per_element!? If it the former, can you add a check for ones_field that the field contains all ones or refactor this to remove that argument?
| Create a regridder between two ClimaCore Spaces. | ||
| This works by finding the vertices of the elements of the source and |
There was a problem hiding this comment.
| Create a regridder between two ClimaCore Spaces. | |
| This works by finding the vertices of the elements of the source and | |
| Create a regridder between two ClimaCore spaces. | |
| This works by finding the vertices of the elements of the source and |
| function Remapping.integrate_each_element(field) | ||
| space = axes(field) |
There was a problem hiding this comment.
Can you add a type annotation for field?
| - `value_per_element_src`: Pre-allocated buffer for source values | ||
| - `value_per_element_dst`: Pre-allocated buffer for destination values | ||
| - `ones_src`: Pre-allocated field of ones on the source space |
There was a problem hiding this comment.
Is it possible to encapsulate this in the regridder object? It seems error prone to expect an user to make this NamedTuple. Alternatively, is it possible to make a helper function for constructing this?
| function ConservativeRegridding.Regridder( | ||
| dst_space::Spaces.SpectralElementSpace2D, | ||
| src_space::Spaces.SpectralElementSpace2D; | ||
| kwargs..., |
There was a problem hiding this comment.
What are the keyword arguments? I couldn't find them.
| function ConservativeRegridding.regrid!( | ||
| dst_field::Fields.Field, | ||
| regridder::ConservativeRegridding.Regridder, | ||
| src_field::Fields.Field, | ||
| ) |
There was a problem hiding this comment.
I am not sure if it is worth doing, but if you think the backend (what interpolator is being used) will change, you can put this regridder in another struct and have the user be exposed to that as the public interface. Then, developers can change what is happening behind the scene without breaking anything.
2 participants
Prototypes regridding using ConservativeRegridding with FV approximation
Closes #2398
Content
get_element_vertices,integrate_each_element,get_value_per_element!,set_value_per_element!regrid!for ClimaCore fields - think about intermediate fieldsNew method ambiguities
This PR increases the number of method ambiguities because ConservativeRegridding.jl pulls in GeometryOps, leading to more ambiguities of
broadcasted. For example, in GeometryOps, they definebroadcasted(f, a::AbstractArray{T}, b::GeometryOps.UnitSpherical.UnitSphericalPoint)and in ClimaCore, we definebroadcasted(op::ClimaCore.Operators.FiniteDifferenceOperator, args...). Thank you to @ph-kev for finding this.A note on zero-area elements when using an odd number of elements
When using a space with odd
h_elem(but not even) we encounter zero-area polygons. This is because when we have odd h_elem, there are elements centered around the north and south poles. Their vertices each lie along a line of constant latitude. When we then compute planar area of this element, we get 0 area, but if we were to compute the correct spherical area it would be nonzero.ConservativeRegridding.jl currently computes simple planar areas, so it currently fails when using a space with even
h_elem. This should be resolved once we switch to computing spherical areas (coming soon).Thank you to @imreddyTeja and @ph-kev for tracking this down!
Prototype regridding latitude
Source field
Source field with one value per element
Destination field (remapped)