conservation-laws / ryujin

High-performance high-order finite element solver for hyperbolic conservation equations
https://conservation-laws.org
Other
103 stars 25 forks source link

`SparsityPatternSIMD`: incorrect exchange of SIMD vectorized to serial vector-valued entry #46

Closed tamiko closed 6 months ago

tamiko commented 6 months ago

While debugging #44 I noticed that we have a rather nasty bug in our ghost row exchange: When exchanging the c_ij matrix via SparseMatrixSIMD::update_ghost_rows() we end up with an incorrect vector-valued matrix entry.

This is hard to spot because (a) scalar, non vector-valued matrices are exchanged correctly; (b) we never used the "ghost row" range of the c_ij matrix for symmetrizing d_ij due to another bug fixed by #43. But now we do.

Example: when symmetrizing d_ij I encounter the following:

rank 0: boundary symmetrization beteen i=4 and j=27
    U_i=0.9999999745486675 0 0 3.00094432415411 - U_j=0.9999999995620856 0 0 3.000154764003423
    c_ij=0.2083333333333333 -0.1041666666666667 - c_ji=-0.2083333333333333 -0.1041666666666667
    norm_ij=0.2329237476562281 - norm_ji=0.2329237476562281
    d_ij=0.3019511981394248 - d_ji=0.3019511981394248
rank 1: boundary symmetrization beteen i=11 and j=24
    U_i=0.9999999995620856 0 0 3.000154764003423 - U_j=0.9999999745486675 0 0 3.00094432415411
    c_ij=-0.2083333333333333 -0.1041666666666667 - c_ji=0.2083333333333333 0.2083333333333333
                                                                           ~~~~~~~~~~~~~~~~~~
    norm_ij=0.2329237476562281 - norm_ji=0.2946278254943948
    d_ij=0.3019511981394248 - d_ji=0.3819414113349516

One can see the c_ji matrix entry on rank 1 is wrong. It had been incorrectly imported from a vectorized range on rank0 to a serial non-vectorized range on rank1.

Looking at the exchange code in question: https://github.com/conservation-laws/ryujin/blob/development/source/sparse_matrix_simd.h#L592-L631 we do not slice correctly.

tamiko commented 6 months ago

@kronbichler ping