Cubed-sphere Jacobian and wind interpolation #99
Closed
Conversation
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
…sphere_interpolation
|
Please integrate the develop branch as PR #98 is merged in now. |
Done! |
* develop: Quad2D: Added fused multiply-add to discriminant calculation (ecmwf#102) FieldSet::has() to replace FieldSet::has_field() Config::json() function Create Array using ArraySpec only (containing datatype) Add Datatype to ArraySpec Cosmetic changes to UnstructuredBilinearLonLat Change Quad2D's PointXY arguments to Point2 Make GridToolsArray backend work with mixed indexing types
Codecov Report
@@ Coverage Diff @@
## develop #99 +/- ##
===========================================
+ Coverage 77.93% 78.02% +0.08%
===========================================
Files 781 781
Lines 53503 53726 +223
===========================================
+ Hits 41696 41918 +222
- Misses 11807 11808 +1
Continue to review full report at Codecov.
|
|
@odlomax Again a very good contribution and high quality code! |
|
@wdeconinck Thank you for your kind words, you've been a great help! Squash and merge away! 😀 |
wdeconinck
approved these changes
Apr 8, 2022
| @@ -82,7 +82,7 @@ class Jacobian : public std::array<std::array<double, 2>, 2> { | |||
| return Jacobian{(*this)[1][1], -(*this)[0][1], -(*this)[1][0], (*this)[0][0]} * (1. / determinant()); | |||
| } | |||
|
|
|||
| Jacobian transpose() const { return Jacobian{(*this)[1][1], (*this)[0][1], (*this)[1][0], (*this)[0][0]}; } | |||
| Jacobian transpose() const { return Jacobian{(*this)[0][0], (*this)[1][0], (*this)[0][1], (*this)[1][1]}; } | |||
There was a problem hiding this comment.
Yes, that was a bit of a head-scratcher! 😂
| // wind transform. Then the transform is applied to the entire field, | ||
| // *including* the halo. | ||
|
|
||
| sourceFunctionspace.parallel_for(util::Config("include_halo", true), [&](idx_t idx, idx_t t, idx_t i, idx_t j) { |
There was a problem hiding this comment.
According @MarekWlasak I'm obsessed with std::tie. I think he's right! 😂
wdeconinck
pushed a commit
that referenced
this pull request
Apr 8, 2022
|
Squash-merged in develop with commit 1b08c23 |
wdeconinck
added a commit
that referenced
this pull request
Apr 21, 2022
* release/0.29.0: (49 commits) Version 0.29.0 GitHub Actions CI: enable OpenMP for C/C++ ATLAS-355 Workaround macOS OpenMP problems with 'omp task' ATLAS-354 bamboo using proj/8.2.1 ATLAS-354 Introduce FindPROJ.cmake for proj installations based on autotools ATLAS-354 Compatibility with proj >= 8 Add missing include <string> Restore compatibility with eckit < 1.18.5 Support cubed-sphere wind interpolation (#99) Fix Jacobian::transpose() introduced in PR #93 GitHub actions CI: reduce available MPI_SLOTS for macOS builds to 4 Quad2D: Added fused multiply-add to discriminant calculation (#102) FieldSet::has() to replace FieldSet::has_field() Config::json() function Create Array using ArraySpec only (containing datatype) Add Datatype to ArraySpec Cosmetic changes to UnstructuredBilinearLonLat Change Quad2D's PointXY arguments to Point2 Make GridToolsArray backend work with mixed indexing types Apply clang-format ...
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
This PR implements the
Projection::jacobianmethod for the equiangular cubed-sphere. There are three versions of this function:projection::Jacobian Proection::jacobian(const PointXY& xy)This returns the Jacobian ∂(x, y) / ∂(lon, lat).projection::JacobianProection::jacobian(const PointXY& xy, idx_t t)This returns the Jacobian ∂(x, y) / ∂(lon, lat) for tile t. This is needed to perform transforms on halo points.projection::Jacobian Proection::alphabetaJacobian(const Point2& alphabeta, idx_t t)This returns the Jacobian ∂(alpha, beta) / ∂(lon, lat) for tile t. As above, but rotated into the (alpha, beta) coordinate system.Two tests have been added:
tests/projection/test_cubedsphere_projection.ccThis makes sure that xy + J(x, y) * h(lon, lat) converges quadratically as the step size h(lon, lat) tends to zero.
(added to
tests/interpolation/test_interpolation_cubedsphere.cc). Here, the following operations are performed:alphabetaJacobianmethod.The figures below show the output u and v field from the wind interpolation test.
The two figures below show the absolute interpolation error, i.e. the vector distance between (u_interp, v_interp) and (u_analytic, v_analytic). The colour scale is logarithmic.
error_field_0shows the error when we naively interpolate u and v as if they were scalar fields. Note the excessive errorat the poles (visible near the top of the sphere) where there is a discontinuity in the direction of the (lon, lat) unit vectors.
error_field_1shows the error when we perform the above wind transform. Note the absence of the polar error in the previous field.