Exporter: allow extrapolation into excisions #6056
Conversation
nilsvu
force-pushed
the
extrapolate_excisions
branch
4 times, most recently
from
e8dd821
to
25febf1
Compare
| const double radial_distance_this_block = | ||
| x_logical->get(direction.dimension()) * direction.sign(); | ||
| if (radial_distance_this_block < 1.) { | ||
| continue; |
There was a problem hiding this comment.
If I follow the logic correctly, if we hit this case the point is not inside the excision sphere. Should that be an error or an early return?
| * The point will be reported to belong to one of the blocks that abut the | ||
| * excision sphere and will have a logical coordinate outside the range [-1, 1] | ||
| * in the radial direction. This is useful for extrapolating data into excision | ||
| * regions. |
There was a problem hiding this comment.
What happens if input_point is not inside excision_sphere? Either document the result or state it is considered invalid to call the function in that case.
| continue; | ||
| } | ||
| // The checks above should leave only 1 valid block, so return that | ||
| return make_id_pair(domain::BlockId(block_id), |
There was a problem hiding this comment.
Could there be two valid blocks for points on the extrapolation of the boundary between two blocks into the sphere? The blocks are being pulled from an unordered container, so always choosing the first one isn't deterministic.
src/IO/Exporter/Exporter.cpp
Outdated
| for (size_t s = 0; s < block_logical_coords->size(); ++s) { | ||
| auto& block_logical_coord = (*block_logical_coords)[s]; | ||
| // Process only points inside the excision sphere. These points were found | ||
| // by `block_logical_coordinates` above, have a radial logical coordinate |
There was a problem hiding this comment.
block_logical_coordinates doesn't appear anywhere above.
| * \cite Etienne2008re adjusted for distorted excisions: we choose uniformly | ||
| * spaced radial anchor points spaced as $\Delta r = 0.3 r_\mathrm{AH}$ in the | ||
| * grid frame (where the excision is spherical), then map the anchor points to | ||
| * the distorted frame (where we have the target point) and do a 7th order |
There was a problem hiding this comment.
The code has constexpr size_t extrapolation_order = 8;
|
Thanks for your review @wthrowe! I pushed a fixup. |
|
I'm investigating the NaNs in some of the unit tests, but the fixup is ready for your review @wthrowe |
| // range [-1, 1] | ||
| // range [-1, 1]. Also discard if the point is on the upper angular face of | ||
| // the block, to disambiguate which of the two blocks sharing the face the | ||
| // point is in. |
There was a problem hiding this comment.
This doesn't work. The blocks don't have aligned frames, so some interfaces are two lower sides or two upper sides.
src/IO/Exporter/Exporter.hpp
Outdated
| * the distorted frame (where we have the target point) and do a 7th order | ||
| * polynomial extrapolation into the excision region. | ||
| * the distorted frame (where we have the target point) and do a polynomial | ||
| * extrapolation into the excision region. |
There was a problem hiding this comment.
It's a 7th order polynomial with 8 collocation points. However in the polynomial_interpolation function the name order is degree + 1, so I used that naming here as well:
To avoid confusion I renamed to
num_extrapolation_anchors.
nilsvu
force-pushed
the
extrapolate_excisions
branch
from
a27ab51
to
3369b36
Compare
nilsvu
force-pushed
the
extrapolate_excisions
branch
from
3369b36
to
779a034
Compare
|
@wthrowe I pushed a new fixup. I dropped the changes to |
|
OK, looks good. Squash and fix the clang-tidy complaint. |
This can be useful to fill the excision region with (non-constraint-satisfying) data so it can be imported into moving puncture codes.
nilsvu
force-pushed
the
extrapolate_excisions
branch
from
779a034
to
af2f718
Compare
|
I don't see a connection 🤷♂️ |
Proposed changes
This is needed to fill the excision region with (non-constraint-satisfying) data so it can be imported into moving puncture codes.
Upgrade instructions
Code review checklist
make docto generate the documentation locally intoBUILD_DIR/docs/html.Then open
index.html.code review guide.
bugfixornew featureif appropriate.Further comments