Skip to content

Add InterfaceDerivativeMatrix for conforming meshes - #498

Draft
PaulineVidal wants to merge 53 commits into
develfrom
pvidal_interface_approx_derivative_matrix
Draft

Add InterfaceDerivativeMatrix for conforming meshes#498
PaulineVidal wants to merge 53 commits into
develfrom
pvidal_interface_approx_derivative_matrix

Conversation

@PaulineVidal

@PaulineVidal PaulineVidal commented Jan 12, 2026

Copy link
Copy Markdown
Collaborator

Operator to compute the coupled interface derivatives along a direction with the approximation formula.

Waiting for #499, better to wait for #500.


Please complete the checklist to ensure that all tasks are completed before marking your pull request as ready for review.

All Submissions

  • Have you ensured that all lines changed in this PR are justified by a comment found in the description ?
  • Have you updated the CHANGELOG.md ?
  • Have you linked any issues that should be closed when this PR is merged (using closing keywords) ?
  • Have you checked that the AUTHORS file is up to date ?
  • Have you checked that the copyright information in the LICENCE file is up to date (including dates) ?
  • Do you follow the conventions specified in our coding standards ?

New Feature Submissions

  • Have you added tests for the new functionalities ?
  • Have you documented the new functionalities:
    • API documentation describing the available methods, when each should be used and how to use them ?
    • User-friendly documentation in README files (which may link to the API documentation).
    • If the new functionality is non-trivial to use, provide a tutorial or example ? (optional)

Changes to Existing Features

  • Have you checked that existing tests cover all code after the changes ?
  • Have you checked that existing tests are still passing ?
  • Have you checked that the existing documentation is still accurate (API and README files) ?

Changes to the CI

  • Have you made the same changes to both the GitHub CI and the GitLab CI (for the private fork) ?

PaulineVidal and others added 16 commits December 10, 2025 11:30
Add `SingleInterfaceDerivativesCalculator` operator and tests to compute an interface derivative with the different cases:

- Interpolation type: Hermite boundary conditions ; additional interpolation points
- Mesh type: uniform ; non-uniform
- Formula type: exact ; approximation. 
- Patch connections: E|E ; E|W ; W|E ; W|W ; E|S

See merge request gysela-developpers/gyselalibxx!869

--------------------------------------------

Co-authored-by: Emily Bourne <emily.bourne@epfl.ch>
Add a `SingleInterfaceDerivativeCalculatorCollection` to store a collection of `SingleInterfaceDerivativeCalculator`.

Change the template parameters of `SingleInterfaceDerivativeCalculator`. The boundary conditions are now inputs in the constructor.  

Move useful operators in the test files `single_interface_derivatives_calculator_test.cpp` in a `interface_derivatives_utils.hpp` file.

See merge request gysela-developpers/gyselalibxx!874

--------------------------------------------

Co-authored-by: Pauline Vidal <142578924+PaulineVidal@users.noreply.github.com>
Co-authored-by: Emily Bourne <emily.bourne@epfl.ch>
Add `SingleInterfaceDerivativesCalculator` operator and tests to compute an interface derivative with the different cases:

- Interpolation type: Hermite boundary conditions ; additional interpolation points
- Mesh type: uniform ; non-uniform
- Formula type: exact ; approximation. 
- Patch connections: E|E ; E|W ; W|E ; W|W ; E|S

See merge request gysela-developpers/gyselalibxx!869

--------------------------------------------

Co-authored-by: Emily Bourne <emily.bourne@epfl.ch>
Add a `SingleInterfaceDerivativeCalculatorCollection` to store a collection of `SingleInterfaceDerivativeCalculator`.

Change the template parameters of `SingleInterfaceDerivativeCalculator`. The boundary conditions are now inputs in the constructor.  

Move useful operators in the test files `single_interface_derivatives_calculator_test.cpp` in a `interface_derivatives_utils.hpp` file.

See merge request gysela-developpers/gyselalibxx!874

--------------------------------------------

Co-authored-by: Pauline Vidal <142578924+PaulineVidal@users.noreply.github.com>
Co-authored-by: Emily Bourne <emily.bourne@epfl.ch>
Add `SingleInterfaceDerivativesCalculator` operator and tests to compute an interface derivative with the different cases:

- Interpolation type: Hermite boundary conditions ; additional interpolation points
- Mesh type: uniform ; non-uniform
- Formula type: exact ; approximation. 
- Patch connections: E|E ; E|W ; W|E ; W|W ; E|S

See merge request gysela-developpers/gyselalibxx!869

--------------------------------------------

Co-authored-by: Emily Bourne <emily.bourne@epfl.ch>
Add a `SingleInterfaceDerivativeCalculatorCollection` to store a collection of `SingleInterfaceDerivativeCalculator`.

Change the template parameters of `SingleInterfaceDerivativeCalculator`. The boundary conditions are now inputs in the constructor.  

Move useful operators in the test files `single_interface_derivatives_calculator_test.cpp` in a `interface_derivatives_utils.hpp` file.

See merge request gysela-developpers/gyselalibxx!874

--------------------------------------------

Co-authored-by: Pauline Vidal <142578924+PaulineVidal@users.noreply.github.com>
Co-authored-by: Emily Bourne <emily.bourne@epfl.ch>
Update the `main_pvidal` branch after merge including DDC changes.

See merge request gysela-developpers/gyselalibxx!882

--------------------------------------------
@PaulineVidal PaulineVidal changed the title Add InterfaceExactDerivativeMatrix Add InterfaceDerivativeMatrix Jan 12, 2026
@github-actions

Copy link
Copy Markdown

This PR is failing tests so it has been put back into draft. Please remove the draft status when the tests pass.

@github-actions
github-actions Bot marked this pull request as draft January 12, 2026 15:12
// Update the derivatives.
// Change the sign of the derivative on Patch1, if the grid directions desagreed.
function_and_derivs_1(idx_slice_deriv_1, idx_slice_1) = interface_deriv * sign_1;
function_and_derivs_2(idx_slice_deriv_2, idx_slice_2) = interface_deriv * sign_2;

@PaulineVidal PaulineVidal Jan 13, 2026

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@EmilyBourne
There is these CI test errors:

 Variable 'function_and_derivs_1(idx_slice_deriv_1,idx_slice_1)' is assigned a value that is never used. [unreadVariable]

see https://github.com/gyselax/gyselalibxx/actions/runs/20954118794/job/60214406603?pr=498

However, these lines modify the values pointed by the Field given as input.

@github-actions
github-actions Bot marked this pull request as draft January 19, 2026 16:44
@github-actions github-actions Bot added the Ready to review Label to be automatically added to a PR when it is ready to be reviewed label Jan 19, 2026
Comment on lines +21 to +23
template <class DimOut, class DimIn>
std::vector<Coord<DimOut>> const convert_dim(std::vector<Coord<DimIn>> const& input_vec)
{

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Doesn't this function do the same thing as fill_in?

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It is true they are similar 🤔
Actually I use fill_in to fill in the vectors of global coordinates. The idea is more like an .append(). I use it in the test files to fill in the global vectors with values of the different patches.
E.g.

fill_in(global_points, patch1_points);
fill_in(global_points, patch2_points);
fill_in(global_points, patch3_points);

Which in python could be equivalent to

global_points = patch1_points
global_points.append(patch2_points)
global_points.append(patch3_points)

It also convert the points to from local coordinates to global coordinates.
Do you think I should rename the operator? I choose fill_in with the idea of filling in a bottle of water.

On the other hand, convert_dim cannot be apply to .append() the vectors.
As it directly returns a vector, I use it to avoid to instantiate a vector.
E.g. (see interface_derivative_matrix_Greville_periodic_test.cpp)

        // Patch 2 ...............................................................................
        ddc::init_discrete_space<BSplinesX<2>>(break_points_x258);
        ddc::init_discrete_space<BSplinesY<2>>(convert_dim<Y<2>, Y<1>>(break_points_y123));

        ddc::init_discrete_space<GridX<2>>(interpolation_points_x258);
        ddc::init_discrete_space<GridY<2>>(convert_dim<Y<2>, Y<1>>(interpolation_points_y123));

        // Patch 3 ...............................................................................
        ddc::init_discrete_space<BSplinesX<3>>(break_points_x369);
        ddc::init_discrete_space<BSplinesY<3>>(convert_dim<Y<3>, Y<1>>(break_points_y123));

        ddc::init_discrete_space<GridX<3>>(interpolation_points_x369);
        ddc::init_discrete_space<GridY<3>>(convert_dim<Y<3>, Y<1>>(interpolation_points_y123));
        
        // ...

Comment on lines +239 to +252
// The patches are conforming between each other.
std::vector<Coord<X<1>>> break_points_x147
= build_random_non_uniform_break_points(x1_min, x1_max, x1_ncells);
std::vector<Coord<X<2>>> break_points_x258
= build_random_non_uniform_break_points(x2_min, x2_max, x2_ncells);
std::vector<Coord<X<3>>> break_points_x369
= build_random_non_uniform_break_points(x3_min, x3_max, x3_ncells);

std::vector<Coord<Y<1>>> break_points_y123
= build_random_non_uniform_break_points(y1_min, y1_max, y1_ncells);
std::vector<Coord<Y<4>>> break_points_y456
= build_random_non_uniform_break_points(y4_min, y4_max, y4_ncells);
std::vector<Coord<Y<7>>> break_points_y789
= build_random_non_uniform_break_points(y7_min, y7_max, y7_ncells);

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If the patches use similar coordinate systems is it possible to have:

using Patch3 = Patch<GridX<3>, GridY<1>, BSplinesX<3>, BSplinesY<1>>

?

@PaulineVidal PaulineVidal Feb 4, 2026

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In the source file, there is

// Sequence of 1D grids in the direction of Grid1D over the patches.
    using Grid1DSeq = collect_grids_on_dim_t<
            find_patch_t<Grid1D, all_patches>,
            Grid1D,
            interface_sorted_collection>;

I am not sure about the result with all_patches = (Patch1, Patch2, ..., Patch9) and Grid1D = GridY<1>.
Does it work?

Maybe it could work. I am not sure. To be tested.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think we should check this at some point

Comment on lines +378 to +393
return get_function_coefficients_with_sign(function_1, function_2, 1);
}

/**
* @brief Get the linear combination of the function values (c).
* See @ref get_function_coefficients.
* @param function_1 Function values at the interpolation points on patch 1.
* @param function_2 Function values at the interpolation points on patch 2.
* @return the linear combination of the function values (c).
*/
template <class Layout1, class Layout2>
inline double get_function_coefficients(
DConstField<IdxRange1DPerp_2, Kokkos::HostSpace, Layout2> const& function_2,
DConstField<IdxRange1DPerp_1, Kokkos::HostSpace, Layout1> const& function_1) const
{
return get_function_coefficients_with_sign(function_1, function_2, 1);

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why is sign hard-coded to 1 in these 2 functions?

@PaulineVidal PaulineVidal Feb 9, 2026

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

To get an interface derivatives, I compute
$\partial_x f = \sum_{k \text{ on Patch1}} \omega_k f_k + \sum_{k \text{ on Patch2}} \omega_k f_k $

To get an interface cross-derivatives, I compute
$\partial_{xy} f = \sum_{k \text{ on Patch1}} \omega_k \partial_y f_k + \sum_{k \text{ on Patch2}} \omega_k \partial_y f_k $

The computed derivatives correspond to the derivatives of an equivalent global spline on the two patches.
In the case of cross-derivatives, we have to be careful that the sign of the $\partial_y f$ agree and correspond to a same equivalent global spline.
For example,

    ^^^^^^^^|^^^^^^^^
    Patch1       Patch2 

can be represented by

    ^^^^^^^^^^^^^^^^
    Global

However, we have a problem with

    ^^^^^^^^|vvvvvvvv
    Patch1       Patch2 

because GridY1 and GridY2 do not have the same orientation. By flipping the sign of the derivative on Patch2, this correspond to orientate GridY2 to the same orientation of GridY1 and bring us to the situation1 the problem.

Of course, the change of sign is not needed for the computation of the first derivative, because the function values are not oriented.

Tell me if this part is not enough https://github.com/gyselax/gyselalibxx/tree/pvidal_interface_approx_derivative_matrix/src/multipatch/interface_derivatives
image

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Of course, the change of sign is not needed for the computation of the first derivative, because the function values are not oriented.

Do you mean 0-th derivative?

@PaulineVidal PaulineVidal Feb 9, 2026

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The first derivatives are computed from the "0th derivatives" = function values -> no need of sign flip.
The cross-derivatives are computed from the 1rst derivatives -> need to flip sign if the grid orientations desagree.

@PaulineVidal PaulineVidal Feb 9, 2026

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe the name of the operator is confusing?

The exact value of the derivative (of a global spline) is
deriv = a x previous_deriv + b x next_deriv + c
with
$c = \sum_{k \text{ on Patch1}} \omega_k f_k + \sum_{k \text{ on Patch2}} \omega_k f_k $
or
$c \sum_{k \text{ on Patch1}} \omega_k \partial_y f_k + \sum_{k \text{ on Patch2}} \omega_k \partial_y f_k $

The coefficient c I called it get_function_coefficient in the first case and get_derivatives_coefficient in the second case, because it is a coefficient computed from the function values or the derivatives in the second case.

Comment on lines +111 to +115
using Interface_1_2 = Interface<EastEdge<1>, WestEdge<2>, true>;
using Interface_2_3 = Interface<EastEdge<2>, SouthEdge<3>, false>;

using OutsideInterface1 = Interface<OutsideEdge, WestEdge<1>, true>;
using OutsideInterface3 = Interface<NorthEdge<3>, OutsideEdge, true>;

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't understand why you have:
using Interface_2_3 = Interface<EastEdge<2>, SouthEdge<3>, false>;
image
Doesn't this interface connect $x_2$ to $y_3$? Those grids are aligned in the same direction

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The boolean refers to the orientation of the parallel grid to the interface.
Here, $y_2$ and $x_3$ dont have the same orientation.

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I added a comment here: 81f26af

@EmilyBourne EmilyBourne removed the Ready to review Label to be automatically added to a PR when it is ready to be reviewed label Feb 16, 2026
EmilyBourne and others added 3 commits February 17, 2026 11:31
Reorganise the InterfaceDerivative tests in order to:
- reduce code duplication (reduce overall code quantity)
- group decisions (i.e. all choices related to a given patch
configuration should be found together)
@EmilyBourne

Copy link
Copy Markdown
Member

Is this PR ready for review again?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants