Add regional interpolation #215
Conversation
|
It seems the hash test is failing depending on the compiler, I'll do it differently. |
|
The new test is working. |
|
I have extended the regional interpolation test : accuracy and adjoint tests, for both 1D and 2D fields. |
There was a problem hiding this comment.
This is a very nice addition @benjaminmenetrier and will be very useful!
This is just a first review without having played around with this feature yet, and sorry that it has taken this long.
Just looking through the code I could find some performance / memory issues that need to be addressed. Please have a look.
| std::vector<bool> toFind = {true, !colocatedX, !colocatedY, !colocatedX && !colocatedY}; | ||
| std::vector<size_t> valueToFind = {indexI*sourceNy+indexJ, (indexI+1)*sourceNy+indexJ, | ||
| indexI*sourceNy+(indexJ+1), (indexI+1)*sourceNy+(indexJ+1)}; | ||
| std::vector<int> foundIndex(4, -1); |
There was a problem hiding this comment.
Since these are always size 4, it will be much better to use std::array<type,4> instead of std::vector<type>
| } else { | ||
| indexJ = sourceNy-1; | ||
| } | ||
| std::cout << "WARNING: point outside of the domain" << std::endl; |
There was a problem hiding this comment.
Please use "Log::info()" or "Log::error()"
| std::vector<int> targetRecvCounts_; | ||
| std::vector<int> targetRecvDispls_; | ||
| std::vector<size_t> sourceSendMapping_; | ||
| std::vector<atlas::interpolation::element::InterpElement> horInterp_; |
There was a problem hiding this comment.
std::vector<atlas::interpolation::element::InterpElement> in current implementation of InterpElement is essentially the same as
std::vector<std::vector<std::pair<size_t,double>>>
The issue is std::vector<std::vector< ... >> because that infers multiple allocations on the heap within tight loop and slow memory access.
I'm thinking to avoid the InterpElement class altogether (and remove the file etc.)
What would be better is multiple arrays:
std::vector< std::array<size_t,4> > stencil_;
std::vector< std::array<double,4> > weights_;
std::vector< size_t > stencil_size_;and then resize all these arrays to targetSize_ before looping over target points, and fill the values at the correct place in these arrays, without using the temporary operations variable. This will also allow this loop to become OpenMP multi-threadable.
These 3 vectors are now each one large memory allocation because the memory of std::array<type,4> is known at compile time.
We need the stencil_size_ array to avoid going into possibly over-allocated entries in case of linear interpolation or point value, which will almost never be the case.
There was a problem hiding this comment.
Great suggestion. For regional DA, we often interpolate from a grid to another one with exactly the same rectangular domain, and just a different cell size, so stencil_size_ might be lower than 4 quite often.
| xy(gp.r, XX) = compute_x(gp.i, gp.j); | ||
| xy(gp.r, YY) = compute_y(gp.j); | ||
| if (regional) { | ||
| std::vector<double> lonlatVec(2); |
There was a problem hiding this comment.
| std::vector<double> lonlatVec(2); | |
| std::array<double,2> lonlatVec; |
This avoids heap allocation within the loop over grid points, and should give very nice speedup.
|
Thank you very much @wdeconinck for your very relevant comments. I realize that I am not well aware of memory usage in C++ and should be more careful about it. |
This PR adds a new interpolation for regional structured source grids, available with the factory tag "regional-linear-2d". It uses a generic sparse matrix element class
InterpElementadded ininterpolation/element, that might also be used for vertical interpolation in the future.An associated test has been added, which defines field values on a coarse source grid, interpolates on a refined target grid, and compute the target field hash on each MPI task (2 tasks requested for this test).
The PR also includes a bugfix in
functionspace/detail/StructuredColumns_setup.cc: it seems that the fieldxyis actually corresponding to thelonlatvalues for all structured grids. However, this was not true indevelopfor regional grids. I have implemented a minimal bugfix where the correct values oflonlatare specified in thexyfield for regional grids too.