Add reflecting boundary conditions for Svärd-Kalisch equations #166

JoshuaLampert · 2024-12-06T13:25:06Z

In the case of $\alpha = \gamma = 0$ reflecting boundary conditions for the Svärd-Kalisch equations satisfy an entropy condition. This PR adds this possibility.

I used an LU factorization instead of a Cholesky decomposition for the system matrix for now because it was a bit simpler to implement and I was not 100% sure if the Dirichlet operator matrix $(\eta + D) - P_DD_1\hat{\beta}D_1$ is always guaranteed to be SPD (or rather the [2:(end - 1), (2:end - 1)] block of it). I would need to check that.

Regarding the treatment of the boundary conditions: We apply a Dirichlet operator in the equation for the velocity and explicitly set the boundary conditions v_1 = v_N = 0 enforcing the boundary condition strongly. As I understand, for the mass equation no special treatment of the boundary condition is necessary as we don't have any elliptic operator and a SAT term would simply be zero as the flux vanishes on the boundary (because the velocity vanishes on the boundary). Correct me if I'm wrong, @ranocha.

Convergence study looks a bit odd sometimes (e.g., convergence in $v$ for p = 6), but seems mostly ok:

p = 2

julia> convergence_test("examples/svaerd_kalisch_1d/svaerd_kalisch_1d_basic_reflecting.jl", 4, N = 16, accuracy_order = 2);
[...]
####################################################################################################
l2
η                        v                        D                        
N    error     EOC       N    error     EOC       N    error     EOC       
16   3.82e-01  -         16   2.20e-02  -         16   0.00e+00  -         
32   1.27e-01  1.59      32   4.82e-03  2.19      32   0.00e+00  NaN       
64   4.44e-02  1.52      64   1.15e-03  2.07      64   0.00e+00  NaN       
128  1.58e-02  1.48      128  2.82e-04  2.03      128  0.00e+00  NaN       

mean           1.53      mean           2.10      mean           NaN       
----------------------------------------------------------------------------------------------------
linf
η                        v                        D                        
N    error     EOC       N    error     EOC       N    error     EOC       
16   1.77e+00  -         16   4.06e-02  -         16   0.00e+00  -         
32   8.73e-01  1.02      32   9.30e-03  2.13      32   0.00e+00  NaN       
64   4.33e-01  1.01      64   2.24e-03  2.05      64   0.00e+00  NaN       
128  2.15e-01  1.01      128  5.50e-04  2.03      128  0.00e+00  NaN       

mean           1.01      mean           2.07      mean           NaN       
----------------------------------------------------------------------------------------------------

p = 4

julia> convergence_test("examples/svaerd_kalisch_1d/svaerd_kalisch_1d_basic_reflecting.jl", 4, N = 16, accuracy_order = 4);
[...]
####################################################################################################
l2
η                        v                        D                        
N    error     EOC       N    error     EOC       N    error     EOC       
16   1.62e-01  -         16   2.73e-03  -         16   0.00e+00  -         
32   1.62e-02  3.32      32   2.20e-04  3.64      32   0.00e+00  NaN       
64   1.74e-03  3.21      64   3.04e-05  2.85      64   0.00e+00  NaN       
128  2.21e-04  2.98      128  4.11e-06  2.89      128  0.00e+00  NaN       

mean           3.17      mean           3.12      mean           NaN       
----------------------------------------------------------------------------------------------------
linf
η                        v                        D                        
N    error     EOC       N    error     EOC       N    error     EOC       
16   5.99e-01  -         16   5.12e-03  -         16   0.00e+00  -         
32   5.68e-02  3.40      32   5.87e-04  3.12      32   0.00e+00  NaN       
64   6.22e-03  3.19      64   7.38e-05  2.99      64   0.00e+00  NaN       
128  1.55e-03  2.00      128  9.30e-06  2.99      128  0.00e+00  NaN       

mean           2.87      mean           3.03      mean           NaN       
----------------------------------------------------------------------------------------------------

p = 6

julia> convergence_test("examples/svaerd_kalisch_1d/svaerd_kalisch_1d_basic_reflecting.jl", 4, N = 16, accuracy_order = 6);
####################################################################################################
l2
η                        v                        D                        
N    error     EOC       N    error     EOC       N    error     EOC       
16   7.65e-01  -         16   1.55e-02  -         16   0.00e+00  -         
32   6.70e-02  3.51      32   4.49e-05  8.43      32   0.00e+00  NaN       
64   3.97e-03  4.08      64   2.28e-05  0.98      64   0.00e+00  NaN       
128  3.06e-04  3.70      128  1.52e-06  3.91      128  0.00e+00  NaN       

mean           3.76      mean           4.44      mean           NaN       
----------------------------------------------------------------------------------------------------
linf
η                        v                        D                        
N    error     EOC       N    error     EOC       N    error     EOC       
16   2.63e+00  -         16   4.05e-02  -         16   0.00e+00  -         
32   3.68e-01  2.84      32   1.30e-04  8.28      32   0.00e+00  NaN       
64   2.72e-02  3.76      64   4.13e-05  1.65      64   0.00e+00  NaN       
128  2.37e-03  3.52      128  2.43e-06  4.09      128  0.00e+00  NaN       

mean           3.37      mean           4.68      mean           NaN       
----------------------------------------------------------------------------------------------------

TODO:

Upwind operators

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

src/equations/bbm_bbm_1d.jl

src/equations/svaerd_kalisch_1d.jl

github-actions · 2024-12-06T13:39:15Z

Benchmark Results

	main	`7092265`...	main/7092265877d86b...
bbm_1d/bbm_1d_basic.jl	14.4 ± 0.31 μs	13.8 ± 0.28 μs	1.04
bbm_1d/bbm_1d_fourier.jl	0.529 ± 0.0031 ms	0.217 ± 0.29 ms	2.44
bbm_bbm_1d/bbm_bbm_1d_basic_reflecting.jl	0.122 ± 0.0032 ms	0.121 ± 0.0033 ms	1.01
bbm_bbm_1d/bbm_bbm_1d_dg.jl	0.0345 ± 0.00053 ms	0.0345 ± 0.00055 ms	1
bbm_bbm_1d/bbm_bbm_1d_relaxation.jl	27.4 ± 0.52 μs	27.4 ± 0.49 μs	1
bbm_bbm_1d/bbm_bbm_1d_upwind_relaxation.jl	0.0485 ± 0.00061 ms	0.0485 ± 0.00057 ms	1
hyperbolic_serre_green_naghdi_1d/hyperbolic_serre_green_naghdi_dingemans.jl	4.11 ± 0.012 μs	4.09 ± 0.013 μs	1.01
serre_green_naghdi_1d/serre_green_naghdi_well_balanced.jl	0.197 ± 0.0047 ms	0.201 ± 0.0046 ms	0.98
svaerd_kalisch_1d/svaerd_kalisch_1d_dingemans_relaxation.jl	0.149 ± 0.0049 ms	0.15 ± 0.0039 ms	0.995
time_to_load	2.01 ± 0.036 s	2.06 ± 0.026 s	0.976

Benchmark Plots

A plot of the benchmark results have been uploaded as an artifact to the workflow run for this PR.

codecov · 2024-12-06T13:52:18Z

Codecov Report

Attention: Patch coverage is 99.19355% with 1 line in your changes missing coverage. Please review.

Project coverage is 98.02%. Comparing base (4d1cc6c) to head (7092265).
Report is 1 commits behind head on main.

Files with missing lines	Patch %	Lines
src/equations/svaerd_kalisch_1d.jl	99.09%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #166      +/-   ##
==========================================
+ Coverage   97.97%   98.02%   +0.05%     
==========================================
  Files          19       19              
  Lines        1776     1877     +101     
==========================================
+ Hits         1740     1840     +100     
- Misses         36       37       +1

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

src/equations/svaerd_kalisch_1d.jl

JoshuaLampert · 2024-12-08T20:48:32Z

Some test values change significantly between different operating systems. In this example we would need an absolute tolerance of ~1.0e-4 to let the MacOS test pass using the test values I get locally, while CI is successful for Linux and Windows. Do you know why these tests are particularly sensitive to rounding errors @ranocha? The condition number of the system matrix doesn't seem to be too high:

julia> include("examples/svaerd_kalisch_1d/svaerd_kalisch_1d_basic_reflecting.jl");
[...]
julia> cond(Matrix(DispersiveShallowWater.assemble_system_matrix!(semi.cache, ones(N), equations)))
107031.29687711073

What do you think is the best way to deal with this? I am a bit hesitant to just put such a high absolute tolerance in the test.
Edit: Could it be that using a Cholesky decomposition helps stabilizing the method?

src/equations/svaerd_kalisch_1d.jl

Co-authored-by: Hendrik Ranocha <[email protected]>

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

ranocha · 2024-12-10T12:24:14Z

Some test values change significantly between different operating systems. In this example we would need an absolute tolerance of ~1.0e-4 to let the MacOS test pass using the test values I get locally, while CI is successful for Linux and Windows.

It also fails on Linux with Julia 1.11, see https://github.com/JoshuaLampert/DispersiveShallowWater.jl/actions/runs/12222702421/job/34093312682?pr=166

Do you know why these tests are particularly sensitive to rounding errors @ranocha? The condition number of the system matrix doesn't seem to be too high:

This is for the central SBP operators, isn't it? We know that these methods are not particularly stable for elliptic solves.

What do you think is the best way to deal with this? I am a bit hesitant to just put such a high absolute tolerance in the test.

I agree. However, if upwind is fine, we may need to accept it 🤷 Does it help if you set stricter tolerances of the time integration methods?

Edit: Could it be that using a Cholesky decomposition helps stabilizing the method?

I am not sure

JoshuaLampert · 2024-12-10T12:36:11Z

It also fails on Linux with Julia 1.11, see JoshuaLampert/DispersiveShallowWater.jl/actions/runs/12222702421/job/34093312682?pr=166

Yes, but it's not as severe. A tolerance of ~1e-7 would be enough. But yes, it is still quite high.

This is for the central SBP operators, isn't it? We know that these methods are not particularly stable for elliptic solves.

The tests fail for both central and upwind operators. For the central operators we would a tolerance of 1e-4 and for the upwind operators 1e-5. The condition number I posted is for central operators.

I agree. However, if upwind is fine, we may need to accept it 🤷 Does it help if you set stricter tolerances of the time integration methods?

Good idea. I'll try.

Edit: Could it be that using a Cholesky decomposition helps stabilizing the method?

I am not sure

Maybe it's worth a try. A Cholesky decomposition should be desirable anyway. And I think the matrix is SPD. I'll give it a try.

JoshuaLampert · 2024-12-10T15:22:40Z

I just found a bug in the upwind discretization, which was introduced in #146. I accidentally pushed the fix 4d1cc6c to the main branch (sorry, I thought I was on another branch). I recognized this bug, when running the following elixir:

using OrdinaryDiffEq
using DispersiveShallowWater
using SummationByPartsOperators: upwind_operators, periodic_derivative_operator

###############################################################################
# Semidiscretization of the Svärd-Kalisch equations

equations = SvaerdKalischEquations1D(gravity_constant = 1.0, eta0 = 0.0,
                                     alpha = 0.0004040404040404049,
                                     beta = 0.49292929292929294,
                                     gamma = 0.15707070707070708)

initial_condition = initial_condition_manufactured
source_terms = source_terms_manufactured
boundary_conditions = boundary_condition_periodic

# create homogeneous mesh
coordinates_min = 0.0
coordinates_max = 1.0
N = 128
mesh = Mesh1D(coordinates_min, coordinates_max, N)

# create solver with periodic upwind SBP operators of accuracy order 4
accuracy_order = 4
D1 = upwind_operators(periodic_derivative_operator; derivative_order = 1,
                      accuracy_order = accuracy_order, xmin = mesh.xmin, xmax = mesh.xmax,
                      N = mesh.N)
D2 = periodic_derivative_operator(2, accuracy_order, mesh.xmin, mesh.xmax, mesh.N)
solver = Solver(D1, D2)

# semidiscretization holds all the necessary data structures for the spatial discretization
semi = Semidiscretization(mesh, equations, initial_condition, solver,
                          boundary_conditions = boundary_conditions,
                          source_terms = source_terms)

###############################################################################
# Create `ODEProblem` and run the simulation
tspan = (0.0,  1.0)
ode = semidiscretize(semi, tspan)
summary_callback = SummaryCallback()
analysis_callback = AnalysisCallback(semi; interval = 10000,
                                     extra_analysis_errors = (:conservation_error,),
                                     extra_analysis_integrals = (waterheight_total,
                                                                 velocity, entropy))
callbacks = CallbackSet(analysis_callback, summary_callback)

saveat = range(tspan..., length = 100)
sol = solve(ode, Tsit5(), abstol = 1e-12, reltol = 1e-12,
            save_everystep = false, callback = callbacks, saveat = saveat)

Before the patch:

julia> trixi_include(joinpath("elixirs", "svaerd_kalisch_1d_manufactured.jl"), tspan = (0.0, 0.01));
[...]
────────────────────────────────────────────────────────────────────────────────────────────────────
Simulation running 'SvaerdKalischEquations1D-BathymetryVariable' with 'initial_condition_manufactured'
────────────────────────────────────────────────────────────────────────────────────────────────────
 #timesteps:              10983                run time:       2.95456787e+00 s
 Δt:             2.43355447e-07                └── GC time:    4.91423915e-01 s (16.633%)
 sim. time:      1.00000000e-02 (100.000%)
 #DOF:                      128                alloc'd memory:        255.540 MiB

 Variable:       η                v                D             
 L2 error:       9.13215214e-07   1.60475943e-08   0.00000000e+00
 Linf error:     1.61551232e-06   2.80500931e-08   0.00000000e+00
 |∫q - ∫q₀|:     1.12757026e-16   1.40679864e-09   0.00000000e+00

 Integrals:    
 ∫η          :   1.23165367e-16
 ∫v          :  -1.40679864e-09
 ∫U          :   1.51761312e+00
────────────────────────────────────────────────────────────────────────────────────────────────────
──────────────────────────────────────────────────────────────────────────────────────
           DispersiveSWE                     Time                    Allocations      
                                    ───────────────────────   ────────────────────────
         Tot / % measured:               2.95s /  97.2%           5.20GiB / 100.0%    

Section                     ncalls     time    %tot     avg     alloc    %tot      avg
──────────────────────────────────────────────────────────────────────────────────────
rhs!                         65.9k    2.87s  100.0%  43.5μs   5.20GiB  100.0%  82.7KiB
  solving elliptic system    65.9k    1.79s   62.3%  27.1μs   3.61GiB   69.3%  57.3KiB
  source terms               65.9k    604ms   21.1%  9.17μs   54.3MiB    1.0%     864B
  assemble system matrix     65.9k    392ms   13.7%  5.95μs   1.54GiB   29.6%  24.5KiB
  deta hyperbolic            65.9k   38.6ms    1.3%   586ns     0.00B    0.0%    0.00B
  dv hyperbolic              65.9k   25.1ms    0.9%   380ns     0.00B    0.0%    0.00B
  ~rhs!~                     65.9k   21.2ms    0.7%   321ns   1.84KiB    0.0%    0.03B
analyze solution                 2    305μs    0.0%   153μs   47.8KiB    0.0%  23.9KiB
──────────────────────────────────────────────────────────────────────────────────────

After the patch:

julia> trixi_include(joinpath("elixirs", "svaerd_kalisch_1d_manufactured.jl"), tspan = (0.0, 0.01));
[...]
────────────────────────────────────────────────────────────────────────────────────────────────────
Simulation running 'SvaerdKalischEquations1D-BathymetryVariable' with 'initial_condition_manufactured'
────────────────────────────────────────────────────────────────────────────────────────────────────
 #timesteps:               4092                run time:       1.02080916e+00 s
 Δt:             9.19401860e-07                └── GC time:    7.80598360e-02 s (7.647%)
 sim. time:      1.00000000e-02 (100.000%)
 #DOF:                      128                alloc'd memory:         94.419 MiB

 Variable:       η                v                D             
 L2 error:       9.24533471e-07   1.60355836e-08   0.00000000e+00
 Linf error:     1.63662047e-06   2.79071288e-08   0.00000000e+00
 |∫q - ∫q₀|:     2.39391840e-16   1.26550453e-09   0.00000000e+00

 Integrals:    
 ∫η          :   2.49800181e-16
 ∫v          :  -1.26550453e-09
 ∫U          :   1.51761312e+00
────────────────────────────────────────────────────────────────────────────────────────────────────
──────────────────────────────────────────────────────────────────────────────────────
           DispersiveSWE                     Time                    Allocations      
                                    ───────────────────────   ────────────────────────
         Tot / % measured:               1.02s /  97.0%           2.04GiB /  99.9%    

Section                     ncalls     time    %tot     avg     alloc    %tot      avg
──────────────────────────────────────────────────────────────────────────────────────
rhs!                         25.8k    989ms   99.9%  38.3μs   2.04GiB  100.0%  82.7KiB
  solving elliptic system    25.8k    613ms   61.9%  23.7μs   1.41GiB   69.3%  57.3KiB
  source terms               25.8k    234ms   23.6%  9.05μs   21.3MiB    1.0%     864B
  assemble system matrix     25.8k    110ms   11.1%  4.24μs    618MiB   29.6%  24.5KiB
  deta hyperbolic            25.8k   14.6ms    1.5%   564ns     0.00B    0.0%    0.00B
  dv hyperbolic              25.8k   9.76ms    1.0%   378ns     0.00B    0.0%    0.00B
  ~rhs!~                     25.8k   8.87ms    0.9%   343ns   1.84KiB    0.0%    0.07B
analyze solution                 5    777μs    0.1%   155μs    116KiB    0.0%  23.3KiB
──────────────────────────────────────────────────────────────────────────────────────

Especially, note the difference in the number of timesteps needed.

JoshuaLampert · 2024-12-10T16:05:58Z

I agree. However, if upwind is fine, we may need to accept it 🤷 Does it help if you set stricter tolerances of the time integration methods?

Good idea. I'll try.

Seems like, it indeed really helps. With abstol and reltol equal to 1e-12, we only need test tolerances of 1e-11, which is acceptable.

ranocha · 2024-12-11T08:18:42Z

Sounds good 👍

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

JoshuaLampert · 2024-12-11T09:32:31Z

I have now switched from an LU decomposition to a Cholesky decomposition in 8859d06, but the errors are too high. So it seems there is some bug somewhere. However, the system_matrix seems to be symmetric and positive definite. Do you have an idea, where things go wrong, @ranocha?

src/equations/svaerd_kalisch_1d.jl

ranocha · 2024-12-11T10:06:10Z

Did you perform a convergence test with the new implementation?

JoshuaLampert · 2024-12-11T10:09:14Z

julia> convergence_test("examples/svaerd_kalisch_1d/svaerd_kalisch_1d_basic_reflecting.jl", 4, N = 16, accuracy_order = 6)
[...]
####################################################################################################
l2
η                        v                        D                        
N    error     EOC       N    error     EOC       N    error     EOC       
16   2.42e+00  -         16   1.39e-01  -         16   0.00e+00  -         
32   6.55e-01  1.89      32   2.30e-02  2.60      32   0.00e+00  NaN       
64   1.79e-01  1.87      64   4.38e-03  2.39      64   0.00e+00  NaN       
128  5.63e-02  1.67      128  9.58e-04  2.19      128  0.00e+00  NaN       

mean           1.81      mean           2.39      mean           NaN       
----------------------------------------------------------------------------------------------------
linf
η                        v                        D                        
N    error     EOC       N    error     EOC       N    error     EOC       
16   6.34e+00  -         16   3.53e-01  -         16   0.00e+00  -         
32   2.35e+00  1.43      32   4.31e-02  3.03      32   0.00e+00  NaN       
64   8.58e-01  1.45      64   6.95e-03  2.63      64   0.00e+00  NaN       
128  3.63e-01  1.24      128  1.36e-03  2.35      128  0.00e+00  NaN       

mean           1.38      mean           2.67      mean           NaN       
----------------------------------------------------------------------------------------------------

JoshuaLampert · 2024-12-11T11:26:58Z

Do you mean something like this:

julia> include("examples/svaerd_kalisch_1d/svaerd_kalisch_1d_basic_reflecting.jl");
[...]
julia> t = 0.1
0.1

julia> q = DispersiveShallowWater.compute_coefficients(initial_condition, t, semi)
([1.2214027581601699, 1.2213796755201909, 1.2213104284727068, 1.2211950196350434, 1.2210334533693017, 1.2208257357821906, 1.2205718747247982, 1.2202718797922947, 1.2199257623235684, 1.219533535400799  …  -1.219533535400799, -1.2199257623235684, -1.2202718797922947, -1.2205718747247982, -1.2208257357821906, -1.2210334533693017, -1.2211950196350434, -1.2213104284727068, -1.2213796755201909, -1.2214027581601699], [0.0, 1.3296421832613442e-5, 5.3184682202485554e-5, 0.0001196617657826978, 0.00021272264721833685, 0.00033236029141141103, 0.00047856565391971273, 0.0006513276814696079, 0.000850633312582725, 0.0010764674783165142  …  0.06004296379054353, 0.05348356952863884, 0.04689559306581195, 0.04027927587157583, 0.033634861490834986, 0.0269625955349242, 0.02026272567253702, 0.013535501620532575, 0.006781175134633048, 0.0], [7.0, 6.999848813690903, 6.9993952776209145, 6.9986394603584765, 6.997581476172813, 6.996221485016647, 6.994559692502023, 6.992596349869215, 6.9903317539487535, 6.987766247116534  …  6.987766247116534, 6.9903317539487535, 6.992596349869215, 6.994559692502023, 6.996221485016647, 6.997581476172813, 6.9986394603584765, 6.9993952776209145, 6.999848813690903, 7.0])

julia> dq = similar(q)
([9.7969281e-316, 9.70988004e-316, 9.7969281e-316, 9.70988004e-316, 6.0e-323, 6.518773660173e-311, 6.4e-323, 7.0620045826603e-311, 7.0e-323, 7.6052355051474e-311  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [9.7969281e-316, 9.8260313e-316, 9.7969281e-316, 9.8260313e-316, 2.172928870602e-311, 2.7161597930893e-311, 5.432314405525e-312, 1.0864623630396e-311, 1.6296932855267e-311, 5.180654e-318  …  4.243991582e-313, 0.0, 3.75e-321, 2.7585945283e-313, 0.0, 3.814e-321, 2.9707941074e-313, 0.0, 4.076e-321, 1.9097962119e-313], [1.5993478e-316, 9.8260349e-316, 9.8260349e-316, 0.0, 9.8260349e-316, 3.78818866e-316, 3.78818905e-316, 3.78818905e-316, 9.82603764e-316, 9.82603764e-316  …  2.53e-321, 2.100775833056e-312, 2.0825436e-316, 6.50493647688795e-310, 3.5e-323, 2.54e-321, 2.100775833056e-312, 2.7170591e-316, 6.50493647688795e-310, 3.5e-323])

julia> DispersiveShallowWater.rhs!(dq, q, t, mesh, equations, initial_condition, boundary_conditions, source_terms, solver, semi.cache)

julia> dq_lu = similar(q)
([8.487983164e-314, 4.14173848e-316, 4.243991582e-314, 7.2792943613299e-310, NaN, 4.14174006e-316, 4.243991582e-314, 8.487983164e-314, 8.4879831856e-314, 4.14174164e-316  …  NaN, 1.586e-321, 7.42420905e-316, 7.4196439e-316, NaN, NaN, NaN, NaN, NaN, NaN], [2.5e-323, 6.7765372e-316, 1.3190567232077176e-84, 4.0e-323, 0.0, 5.0e-324, 0.0, 1.69759663287e-313, 2.121995791e-314, 0.0  …  0.0, 6.96e-321, 6.50494391036804e-310, 4.92504477e-316, 0.0, 0.0, 3.28093665e-316, 6.94252987e-316, 4.03878433e-316, 9.8453521e-316], [6.50494391039334e-310, 4.1419547e-316, 4.1419547e-316, 4.1419547e-316, 1.497e-321, 5.0e-324, 0.0, 0.0, 1.7381080806e-314, 6.5049397079524e-310  …  4.0387843e-316, 3.8132271e-316, NaN, NaN, NaN, 3.00831513e-316, NaN, 4.14055865e-316, NaN, 4.38453854e-316])

julia> DispersiveShallowWater.rhs_lu!(dq_lu, q, t, mesh, equations, initial_condition, boundary_conditions, source_terms, solver, semi.cache)

julia> dq .- dq_lu
([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 2.5885250037231663e-5, -1.62061519686071e-5, 3.576633358458717e-5, -1.0175667095673287e-5, 3.6540107521403726e-5, -1.0095704447866147e-5, 3.653070993839909e-5, -1.011872634096723e-5, 3.650777580948075e-5  …  -3.877105139398623e-5, 1.1622132469489566e-5, -3.8796828582446374e-5, 1.1602498280804951e-5, -3.880869896399214e-5, 1.1689689387396618e-5, -3.7999106218711276e-5, 1.8091870581942934e-5, -2.7143511910356424e-5, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])

There are differences of the order of 1e-5 in the velocity.

JoshuaLampert · 2024-12-11T14:33:06Z

Ohh, I found the bug (f017619 fixes this). scale_by_mass_matrix! was used in solve_system_matrix! for a smaller vector than the size of the mass matrix. Why doesn't this trigger the boundscheck in https://github.com/ranocha/SummationByPartsOperators.jl/blob/00164c715f6df4f1ea7297d28e2a39ade3334223/src/matrix_operators.jl#L55?
Let me clean up this PR. Then this should be fine.

ranocha · 2024-12-11T15:09:48Z

Ohh, I found the bug (f017619 fixes this).

🥳

scale_by_mass_matrix! was used in solve_system_matrix! for a smaller vector than the size of the mass matrix. Why doesn't this trigger the boundscheck in https://github.com/ranocha/SummationByPartsOperators.jl/blob/00164c715f6df4f1ea7297d28e2a39ade3334223/src/matrix_operators.jl#L55?

That's a good question...

JoshuaLampert · 2024-12-11T15:27:44Z

scale_by_mass_matrix! was used in solve_system_matrix! for a smaller vector than the size of the mass matrix. Why doesn't this trigger the boundscheck in ranocha/SummationByPartsOperators.jl@00164c7/src/matrix_operators.jl#L55?

That's a good question...

I opened an issue: ranocha/SummationByPartsOperators.jl#310

JoshuaLampert · 2024-12-11T15:33:49Z

This is ready for a review @ranocha.

ranocha

Thanks! This looks already quite good to me 👍

examples/svaerd_kalisch_1d/svaerd_kalisch_1d_basic_reflecting.jl

src/equations/equations.jl

src/equations/svaerd_kalisch_1d.jl

Co-authored-by: Hendrik Ranocha <[email protected]>

src/equations/equations.jl

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

ranocha

Thanks!

add reflecting boundary conditions for Svärd-Kalisch equations

3db9929

JoshuaLampert added the enhancement New feature or request label Dec 6, 2024

Apply suggestions from code review

94dad02

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

JoshuaLampert commented Dec 6, 2024

View reviewed changes

src/equations/bbm_bbm_1d.jl Show resolved Hide resolved

src/equations/svaerd_kalisch_1d.jl Show resolved Hide resolved

add back accidently removed FourierDerivativeOperator

0d5ca48

JoshuaLampert added 2 commits December 8, 2024 15:38

also implement upwind operators

c7f3d4e

reformat

9e9d472

JoshuaLampert commented Dec 8, 2024

View reviewed changes

src/equations/svaerd_kalisch_1d.jl Show resolved Hide resolved

JoshuaLampert commented Dec 8, 2024

View reviewed changes

src/equations/svaerd_kalisch_1d.jl Show resolved Hide resolved

ranocha reviewed Dec 10, 2024

View reviewed changes

src/equations/svaerd_kalisch_1d.jl Outdated Show resolved Hide resolved

src/equations/svaerd_kalisch_1d.jl Outdated Show resolved Hide resolved

src/equations/svaerd_kalisch_1d.jl Outdated Show resolved Hide resolved

src/equations/svaerd_kalisch_1d.jl Outdated Show resolved Hide resolved

JoshuaLampert and others added 2 commits December 10, 2024 12:38

Apply suggestions from code review

874028b

Co-authored-by: Hendrik Ranocha <[email protected]>

Update src/equations/svaerd_kalisch_1d.jl

25ad669

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

Merge branch 'main' into reflecting-bc-SK

d8149c6

JoshuaLampert added 2 commits December 10, 2024 16:40

set stricter time integration tolerances

8fe7401

lower tolerance

75ec37f

JoshuaLampert and others added 2 commits December 11, 2024 10:29

use Cholesky decomposition

8859d06

Update src/equations/svaerd_kalisch_1d.jl

40468f3

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

ranocha reviewed Dec 11, 2024

View reviewed changes

src/equations/svaerd_kalisch_1d.jl Show resolved Hide resolved

ranocha reviewed Dec 11, 2024

View reviewed changes

src/equations/svaerd_kalisch_1d.jl Show resolved Hide resolved

fix scale_by_mass_matrix!

f017619

JoshuaLampert added 8 commits December 11, 2024 15:37

remove LU decomposition again

ca345a5

scale_by_mass_matrix! outside of solve_system_matrix!

3fa969f

try CI without tolerances

82becee

reorder imports

60e800e

fix variable name

08c04e6

refactor to reduce code redundancy

94991f3

fix

320d883

pass boundary_conditions

a2d8497

format

9bfebe3

adjust tolerances

2d50833

JoshuaLampert marked this pull request as ready for review December 11, 2024 15:33

JoshuaLampert requested a review from ranocha December 11, 2024 15:33

ranocha requested changes Dec 11, 2024

View reviewed changes

examples/svaerd_kalisch_1d/svaerd_kalisch_1d_basic_reflecting.jl Outdated Show resolved Hide resolved

src/equations/equations.jl Outdated Show resolved Hide resolved

src/equations/svaerd_kalisch_1d.jl Show resolved Hide resolved

Apply suggestions from code review

209cc76

Co-authored-by: Hendrik Ranocha <[email protected]>

github-actions bot reviewed Dec 11, 2024

View reviewed changes

src/equations/equations.jl Outdated Show resolved Hide resolved

JoshuaLampert and others added 4 commits December 11, 2024 17:11

Update src/equations/equations.jl

2308121

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

begin + 1instead of 2

3302cd0

set tolerances in tests

7358093

add test for ArgumentError

7092265

JoshuaLampert requested a review from ranocha December 11, 2024 16:45

JoshuaLampert mentioned this pull request Dec 11, 2024

Update changelog #167

Merged

ranocha approved these changes Dec 11, 2024

View reviewed changes

JoshuaLampert merged commit 1a509a9 into main Dec 11, 2024
18 checks passed

JoshuaLampert deleted the reflecting-bc-SK branch December 11, 2024 19:00

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add reflecting boundary conditions for Svärd-Kalisch equations #166

Add reflecting boundary conditions for Svärd-Kalisch equations #166

JoshuaLampert commented Dec 6, 2024 •

edited

Loading

github-actions bot commented Dec 6, 2024 •

edited

Loading

codecov bot commented Dec 6, 2024 •

edited

Loading

JoshuaLampert commented Dec 8, 2024 •

edited

Loading

ranocha commented Dec 10, 2024

JoshuaLampert commented Dec 10, 2024

JoshuaLampert commented Dec 10, 2024 •

edited

Loading

JoshuaLampert commented Dec 10, 2024

ranocha commented Dec 11, 2024

JoshuaLampert commented Dec 11, 2024 •

edited

Loading

ranocha commented Dec 11, 2024

JoshuaLampert commented Dec 11, 2024

JoshuaLampert commented Dec 11, 2024

JoshuaLampert commented Dec 11, 2024 •

edited

Loading

ranocha commented Dec 11, 2024

JoshuaLampert commented Dec 11, 2024

JoshuaLampert commented Dec 11, 2024

ranocha left a comment

ranocha left a comment

Add reflecting boundary conditions for Svärd-Kalisch equations #166

Add reflecting boundary conditions for Svärd-Kalisch equations #166

Conversation

JoshuaLampert commented Dec 6, 2024 • edited Loading

github-actions bot commented Dec 6, 2024 • edited Loading

Benchmark Results

Benchmark Plots

codecov bot commented Dec 6, 2024 • edited Loading

Codecov Report

JoshuaLampert commented Dec 8, 2024 • edited Loading

ranocha commented Dec 10, 2024

JoshuaLampert commented Dec 10, 2024

JoshuaLampert commented Dec 10, 2024 • edited Loading

JoshuaLampert commented Dec 10, 2024

ranocha commented Dec 11, 2024

JoshuaLampert commented Dec 11, 2024 • edited Loading

ranocha commented Dec 11, 2024

JoshuaLampert commented Dec 11, 2024

JoshuaLampert commented Dec 11, 2024

JoshuaLampert commented Dec 11, 2024 • edited Loading

ranocha commented Dec 11, 2024

JoshuaLampert commented Dec 11, 2024

JoshuaLampert commented Dec 11, 2024

ranocha left a comment

Choose a reason for hiding this comment

ranocha left a comment

Choose a reason for hiding this comment

JoshuaLampert commented Dec 6, 2024 •

edited

Loading

github-actions bot commented Dec 6, 2024 •

edited

Loading

codecov bot commented Dec 6, 2024 •

edited

Loading

JoshuaLampert commented Dec 8, 2024 •

edited

Loading

JoshuaLampert commented Dec 10, 2024 •

edited

Loading

JoshuaLampert commented Dec 11, 2024 •

edited

Loading

JoshuaLampert commented Dec 11, 2024 •

edited

Loading