Add higher order derivatives support to the evaluators #973
Conversation
Codecov Report❌ Patch coverage is
❌ Your patch check has failed because the patch coverage (85.00%) is below the target coverage (100.00%). You can increase the patch coverage or adjust the target coverage. Additional details and impacted files@@ Coverage Diff @@
## main #973 +/- ##
==========================================
- Coverage 91.76% 91.63% -0.13%
==========================================
Files 57 57
Lines 2756 2643 -113
Branches 924 865 -59
==========================================
- Hits 2529 2422 -107
+ Misses 136 130 -6
Partials 91 91 ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
| template <class... DDims> | ||
| double error_bound( | ||
| ddc::DiscreteElement<DDims...> const& orders, | ||
| std::array<double, sizeof...(DDims)> const& cell_width, | ||
| std::array<int, sizeof...(DDims)> const& degrees) const | ||
| { | ||
| using ddims = ddc::detail::TypeSeq<DDims...>; | ||
| return (tikhomirov_error_bound( | ||
| cell_width[ddc::type_seq_rank_v<DDims, ddims>], | ||
| degrees[ddc::type_seq_rank_v<DDims, ddims>] - orders.template uid<DDims>(), | ||
| max_norm<DDims>(orders, degrees)) | ||
| + ...); | ||
| } |
There was a problem hiding this comment.
It is outside the scope of this PR, but in the future all the other error_bound_on* functions in this file could be removed in favor of this generic implementation
tpadioleau
left a comment
tpadioleau
left a comment
There was a problem hiding this comment.
Ok looks pretty good, just a few comments/questions.
| if (order1 == 0) { | ||
| jmin1 = ddc::discrete_space<bsplines_type1>() | ||
| .eval_basis(vals1, coord_eval_interest1); | ||
| } else if (order1 == 1) { | ||
| jmin1 = ddc::discrete_space<bsplines_type1>() | ||
| .eval_deriv(vals1, coord_eval_interest1); | ||
| } else { | ||
| std::array<double, 3 * (bsplines_type1::degree() + 1)> derivs1_ptr; | ||
| Kokkos::mdspan< | ||
| double, | ||
| Kokkos::extents<std::size_t, bsplines_type1::degree() + 1, 3>> const | ||
| derivs1(derivs1_ptr.data()); | ||
|
|
||
| jmin1 = ddc::discrete_space<bsplines_type1>() | ||
| .eval_basis_and_n_derivs(derivs1, coord_eval_interest1, order1); | ||
|
|
||
| for (std::size_t i = 0; i < bsplines_type1::degree() + 1; ++i) { | ||
| vals1[i] = DDC_MDSPAN_ACCESS_OP(derivs1, i, order1); | ||
| } | ||
| } |
There was a problem hiding this comment.
I would be in favor of only keeping the cases order1 > 0 && order1 < 3. The case order1 == 0, can be represented by not adding the dimension in the discrete element. I wonder if the case order1 == 1 is necessary, I guess we can consider it an optimization path ?
@EmilyBourne Any opinion ?
There was a problem hiding this comment.
I would be in favor of only keeping the cases order1 > 0 && order1 < 3.
I agree with order1 > 0 as it seems counter-intuitive to use a derivs function to get the values. But I don't see any reason to restrict to <3. Unless you meant <degree()?
I would agree with restricting to order1<degree()+1 as the (d+1)-th derivative of a d-th order polynomial is 0.
There was a problem hiding this comment.
Yes I guess we got confused between the degree and dimensionality. We will fix that.
There was a problem hiding this comment.
I wonder if the case order1 == 1 is necessary, I guess we can consider it an optimization path ?
We could only have a case order > 0, the only advantages I see for order == 1 are that it does a few less operations, and it uses less stack memory (2 arrays in eval_basis_and_n_derivs vs 0 in eval_deriv)
Co-authored-by: Thomas Padioleau <[email protected]>
tretre91
force-pushed
the
nth_order_derivatives
branch
from
3490048 to
96349fb
Compare
- Use the degree of the spline instead of the Evaluator's dimension as the upper limit for the derivation order - Update the tests - Don't split the 3D test into multiple executables anymore
3 participants
This PR adds the ability to compute derivatives higher than the first order using the
SplineEvaluatorclasses.The
deriv,deriv_dim_1,deriv_1_and_2,deriv_1_2_3, etc. methods are deprecated in favor of aderivmethod taking as first argument aDiscreteElementdefined on the dimensions of interest whose components are the desired orders of derivation along each of the dimension.For example, for a 3D evaluator with
continuous_dimension_typesI1,I2andI3will use the second derivative on the first dimension, the first derivative on the third dimension, and the values (no derivation) on the second dimension because the DiscreteElement does not contain a component in
Deriv<I2>.For practical reasons, the 1D evaluator can go up to the first derivative, the 2D evaluator up to the second derivative and the 3D evaluator up to the third derivative.
Note that I still need to create some tests for the higher order derivatives, the current tests only go up to cross-differentiation with the first derivative on all dimensions.
This should resolve #926