Add 10 and 12 order boundary optimized operators of Mattsson, Almquist, van der Weide by xlxs4 · Pull Request #152 · ranocha/SummationByPartsOperators.jl

xlxs4 · 2022-11-01T11:27:43Z

Closes #37

10th order accurate
12th order accurate
10th order minimal
12th order minimal

pending tests

ranocha · 2022-11-01T12:26:59Z

Feel free to ask here if you have any questions along the way. And please ping me when you're ready for a review 🙂

xlxs4 · 2022-11-02T14:57:34Z

Hi @ranocha, thanks!

From a quick look, I can see it's split between /dev/MattssonAlmquistVanDerWeide2018... and /src/MattssonAlmquistVanDerWeide2018....

$x_k$, the non-equidistant grid points in MATLAB (D1_accurate_8th.m):

%%%% Non-equidistant grid points %%%%%
x0 = 0.0000000000000e+00;
x1 = 3.8118550247622e-01;
x2 = 1.1899550868338e+00;
x3 = 2.2476300175641e+00;
x4 = 3.3192851303204e+00;
x5 = 4.3192851303204e+00;
x6 = 5.3192851303204e+00;
x7 = 6.3192851303204e+00;
x8 = 7.3192851303204e+00;

Reside in /src/:

elseif accuracy_order == 8
    xstart = SVector(
      T(0.0000000000000e+00),
      T(3.8118550247622e-01),
      T(1.1899550868338e+00),
      T(2.2476300175641e+00),
      T(3.3192851303204e+00),
      T(4.3192851303204e+00),
      T(5.3192851303204e+00),
      T(6.3192851303204e+00),
      T(7.3192851303204e+00),
    )

The interior stencil of the $Q$ matrix in MATLAB:

switch order
	case 2
		d = [-1/2,0,1/2];
	case 4
		d = [1/12,-2/3,0,2/3,-1/12];
	case 6
		d = [-1/60,3/20,-3/4,0,3/4,-3/20,1/60];
	case 8
		d = [1/280,-4/105,1/5,-4/5,0,4/5,-1/5,4/105,-1/280];
	case 10
		d = [-1/1260,5/504,-5/84,5/21,-5/6,0,5/6,-5/21,5/84,-5/504,1/1260];
	case 12
		d = [1/5544,-1/385,1/56,-5/63,15/56,-6/7,0,6/7,-15/56,5/63,-1/56,1/385,-1/5544];
end

For, e.g. order 8, the relevant part is

	case 8
		d = [1/280,-4/105,1/5,-4/5,0,4/5,-1/5,4/105,-1/280];

which can be found in src/:

upper_coef = SVector(T(4//5), T(-1//5), T(4//105), T(-1//280))
lower_coef = -upper_coef

For the norm matrix, we have:

P0 = 1.0758368078310e-01;
P1 = 6.1909685107891e-01;
P2 = 9.6971176519117e-01;
P3 = 1.1023441350947e+00;
P4 = 1.0244688965833e+00;
P5 = 9.9533550116831e-01;
P6 = 1.0008236941028e+00;
P7 = 9.9992060631812e-01;

which can be found in dev/:

h = [
  1.0758368078310e-01,
  6.1909685107891e-01,
  9.6971176519117e-01,
  1.1023441350947e+00,
  1.0244688965833e+00,
  9.9533550116831e-01,
  1.0008236941028e+00,
  9.9992060631812e-01,
  1, 1, 1, 1]

question: what about the trailing ones?

The boundaries, $Q_k$ are:

...
Q0_0 = -5.0000000000000e-01;
Q0_1 = 6.7284756079369e-01;
Q0_2 = -2.5969732837062e-01;
Q0_3 = 1.3519390385721e-01;
Q0_4 = -6.9678474730984e-02;
Q0_5 = 2.6434024071371e-02;
Q0_6 = -5.5992311465618e-03;
Q0_7 = 4.9954552590464e-04;
Q0_8 = 0.0000000000000e+00;
Q0_9 = 0.0000000000000e+00;
Q0_10 = 0.0000000000000e+00;
Q0_11 = 0.0000000000000e+00;

which also reside in dev/:

q1 = [0,
  6.7284756079369e-01,
  -2.5969732837062e-01,
  1.3519390385721e-01,
  -6.9678474730984e-02,
  2.6434024071371e-02,
  -5.5992311465618e-03,
  4.9954552590464e-04,
  0, 0, 0, 0]'

question: in this example, why is the original Q0_0 = -5.0000000000000e-01;, while it's 0 in the Julia implementation?

question: why D1 = H \ (Q - 1//2*e*e') to construct the difference operator $D1$?

question: I understand that dev/ is used first to get the values used in src/, e.g. for the DerivativeCoefficientRows. How?

ranocha · 2022-11-02T15:21:24Z

dev is just my set of notes while translating the stuff, the real coefficients live in src. You can basically ignore the stuff in dev.
For the derivative operator, I use the coefficients of D1, not of Q. Thus, I usually translated Q coefficients from a paper to D1 coefficients for this package.
The SBP property is M * D1 + D1' * M' = eR * eR' - eL * eL' with eL = [1, 0, ..., 0] and eR = [0, ..., 0, 1]. This allows to translate Q to D1 and vice versa.
You do not need to put anything in dev in his PR. (but you may do so if you like). We just need the coefficients in src.

ranocha · 2022-11-03T05:26:14Z

When you have something that works, you should also add some tests, e.g., to https://github.com/ranocha/SummationByPartsOperators.jl/blob/main/test/SBP_operators_test.jl. I think I set up tests only up to order 8, so there need to be more tests for the higher-order operators.

xlxs4 · 2022-11-03T11:18:31Z

Makes perfect sense. Thanks!

xlxs4 · 2022-11-03T11:31:48Z

What about the weights?

ranocha · 2022-11-04T10:47:54Z

What about the weights?

What exactly do you mean? The trailing ones in dev? I just used them to have a boundary part big enough to reach the inner stencil.

xlxs4 · 2022-11-04T16:10:28Z

No, I mean the weights as in, e.g., here

ranocha · 2022-11-07T09:26:41Z

These weights should be the entries on the diagonal of the SBP mass/norm matrix M.

xlxs4 · 2022-11-07T22:56:31Z

@ranocha I think you can give this a quick review

ranocha

Looks good so far (without checking the coefficients) - thanks a lot! Could you please make sure that the new operators are also tested?

xlxs4 · 2022-11-17T12:42:14Z

Will do soon enough, sorry for the delay

xlxs4 · 2022-12-04T19:47:43Z

Hi! I'm getting most tests to pass, although in a few cases I'm getting some values that are way off.
For example, testing the accurate 10th order operator in the @test all(i->abs(res[i]) < <a_number>*eps(T), eachindex(res)) check, every value is alright, save for the last value in res:

using Test
using LinearAlgebra
using SummationByPartsOperators

source = MattssonAlmquistVanDerWeide2018Accurate()
for T in (Float32, Float64)
    xmin = -one(T)
    xmax = 2*one(T)
    N = 101
    der_order = 1
	acc_order = 10
    D = try
        derivative_operator(source, der_order, acc_order, xmin, xmax, N)
    catch
        nothing
    end
    if D !== nothing
        D = derivative_operator(source, der_order, acc_order, xmin, xmax, N)
        x1 = grid(D)
        x0 = fill(one(eltype(x1)), length(x1))
        x2 = x1 .* x1
        x3 = x2 .* x1
        x4 = x2 .* x2
        x5 = x2 .* x3
        x6 = x3 .* x3
        x7 = x4 .* x3
        res = fill(zero(eltype(x0)), length(x0))
        inner_indices = length(D.coefficients.left_boundary)+1:length(res)-length(D.coefficients.left_boundary)-1
        mul!(res, D, x0)
        for i in eachindex(res)
            res[i] >= 6000 * eps(T) && println("value: $(res[i]) at index $i")
        end
        println(length(res))
    end
end

Testing Running tests...
value: 1647.6987 at index 101
101
value: 1647.6989635229024 at index 101
101

Any ideas?

ranocha · 2022-12-05T07:08:15Z

src/SBP_coefficients/MattssonAlmquistVanDerWeide2018Accurate.jl

+        DerivativeCoefficientRow{T,1,10}(SVector(T(-55),
+                                                T(6.7548747038001995),
+                                                T(-2.6691978151545994),
+                                                T(1.4438714982129999),
+                                                T(-0.7727367375076),
+                                                T(0.25570078343005),
+                                                T(0.042808774693299),
+                                                T(-0.082902108933389),
+                                                T(0.032031176427907995),
+                                                T(-0.0044502749689555995) )),


There seems to be something wrong with these coefficients. They are intended to be the coefficients of the derivative operator D, not the skew-symmetric part of it. Thus, each row needs to sum to zero. This is not the case here:

using StaticArrays T = Float64 SVector(T(-55), T(6.7548747038001995), T(-2.6691978151545994), T(1.4438714982129999), T(-0.7727367375076), T(0.25570078343005), T(0.042808774693299), T(-0.082902108933389), T(0.032031176427907995), T(-0.0044502749689555995)) |> sum

This yields

-50.000000000000085

Note that the first row is basically also used as last row, explaining your observation described above.

I understand what you say and the problem does go away with, e.g. changing that -55 to ~-5. Still, I don't understand how that -55 came to be in the first place. I've cross-checked the coefficients from MATLAB.

Code

using LinearAlgebra # order 10 h = [ 1.0000000000000e-01, 5.8980851260667e-01, 9.5666820955973e-01, 1.1500297411596e+00, 1.1232986993248e+00, 1.0123020150951e+00, 9.9877122702527e-01, 1.0000873322761e+00, 1.0000045540888e+00, 9.9999861455083e-01, 1, 1, 1, 1, 1] e = [1//1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] H = Diagonal(h) q1 = [-5, 6.7548747038002e-01, -2.6691978151546e-01, 1.4438714982130e-01, -7.7273673750760e-02, 2.5570078343005e-02, 4.2808774693299e-03, -8.2902108933389e-03, 3.2031176427908e-03, -4.4502749689556e-04, 0, 0, 0, 0, 0]' q2 = [-6.7548747038002e-01, 0, 9.5146052715180e-01, -4.2442349882626e-01, 2.1538865145190e-01, -7.1939778160350e-02, -8.2539187832840e-03, 1.9930661669090e-02, -7.7433256989613e-03, 1.0681515760869e-03, 0, 0, 0, 0, 0]' q3 = [2.6691978151546e-01, -9.5146052715180e-01, 0, 9.6073770842387e-01, -3.9378595264609e-01, 1.3302097358959e-01, 8.1200458151489e-05, -2.3849770528789e-02, 9.6600442856829e-03, -1.3234579460680e-03, 0, 0, 0, 0, 0]' q4 = [-1.4438714982130e-01, 4.2442349882626e-01, -9.6073770842387e-01, 0, 9.1551097634196e-01, -2.8541713079648e-01, 4.1398809121293e-02, 1.7256059167927e-02, -9.4349194803610e-03, 1.3875650645663e-03, 0, 0, 0, 0, 0]' q5 = [7.7273673750760e-02, -2.1538865145190e-01, 3.9378595264609e-01, -9.1551097634196e-01, 0, 8.3519401865051e-01, -2.0586492924974e-01, 3.1230261235901e-02, -2.0969453466651e-04, -5.0965470499782e-04, 0, 0, 0, 0, 0]' q6 = [-2.5570078343005e-02, 7.1939778160350e-02, -1.3302097358959e-01, 2.8541713079648e-01, -8.3519401865051e-01, 0, 8.1046389580138e-01, -2.1879194972141e-01, 5.2977237804899e-02, -9.0146730522360e-03, 7.9365079365079e-04, 0, 0, 0, 0]' q7 = [-4.2808774693299e-03, 8.2539187832840e-03, -8.1200458151489e-05, -4.1398809121293e-02, 2.0586492924974e-01, -8.1046389580138e-01, 0, 8.2787884456005e-01, -2.3582460382545e-01, 5.9178678209520e-02, -9.9206349206349e-03, 7.9365079365079e-04, 0, 0, 0]' q8 = [8.2902108933389e-03, -1.9930661669090e-02, 2.3849770528789e-02, -1.7256059167927e-02, -3.1230261235901e-02, 2.1879194972141e-01, -8.2787884456005e-01, 0, 8.3301028859275e-01, -2.3804321850015e-01, 5.9523809523809e-02, -9.9206349206349e-03, 7.9365079365079e-04, 0, 0]' q9 = [-3.2031176427908e-03, 7.7433256989613e-03, -9.6600442856829e-03, 9.4349194803610e-03, 2.0969453466651e-04, -5.2977237804899e-02, 2.3582460382545e-01, -8.3301028859275e-01, 0, 8.3333655748509e-01, -2.3809523809524e-01, 5.9523809523809e-02, -9.9206349206349e-03, 7.9365079365079e-04, 0]' q10 = [4.4502749689556e-04, -1.0681515760869e-03, 1.3234579460680e-03, -1.3875650645663e-03, 5.0965470499782e-04, 9.0146730522360e-03, -5.9178678209520e-02, 2.3804321850015e-01, -8.3333655748509e-01, 0, 8.3333333333333e-01, -2.3809523809524e-01, 5.9523809523809e-02, -9.9206349206349e-03, 7.9365079365079e-04]' q11 = [0 0 0 0 -1//1260 5//504 -5//84 5//21 -5//6 0 5//6 -5//21 5//84 -5//504 1//1260] q12 = [0 0 0 0 0 -1//1260 5//504 -5//84 5//21 -5//6 0 5//6 -5//21 5//84 -5//504] q13 = [0 0 0 0 0 0 -1//1260 5//504 -5//84 5//21 -5//6 0 5//6 -5//21 5//84] q14 = [0 0 0 0 0 0 0 -1//1260 5//504 -5//84 5//21 -5//6 0 5//6 -5//21] q15 = [0 0 0 0 0 0 0 0 -1//1260 5//504 -5//84 5//21 -5//6 0 5//6] Q = vcat(q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15) D1 = H \ (Q - 1//2*e*e'); display(D1)

Results

[-55.0, 6.7548747038001995, -2.6691978151545994, 1.4438714982129999, -0.7727367375076, 0.25570078343005, 0.042808774693299, -0.082902108933389, 0.032031176427907995, -0.0044502749689555995, 0.0, 0.0, 0.0, 0.0, 0.0] [-1.145265719198745, 0.0, 1.6131685230292827, -0.7195954106367713, 0.3651840332042441, -0.12197141381091767, -0.01399423475053903, 0.03379174976808331, -0.013128541778312972, 0.0018110141736784769, 0.0, 0.0, 0.0, 0.0, 0.0] [0.27900977459917853, -0.9945564383180177, 0.0, 1.004253824704819, -0.4116222831605484, 0.13904608960593332, 8.487839079429473e-5, -0.024930033516808243, 0.010097590982069497, -0.0013834032874125409, 0.0, 0.0, 0.0, 0.0, 0.0] [-0.12555079634350264, 0.36905436758383703, -0.8354024892044419, 0.0, 0.7960759131488463, -0.2481823909255492, 0.035998033476551373, 0.01500488078728063, -0.008204065636465682, 0.0012065471134400296, 0.0, 0.0, 0.0, 0.0, 0.0] [0.06879174149957458, -0.1917465511011161, 0.3505621014987283, -0.8150200626888124, 0.0, 0.7435190828161148, -0.18326819871996888, 0.027802276682660717, -0.00018667744812003666, -0.0004537125390638899, 0.0, 0.0, 0.0, 0.0, 0.0] [-0.02525933759067232, 0.07106552895046016, -0.1314044342558119, 0.28194859492566227, -0.8250443110814595, 0.0, 0.8006147214131956, -0.2161330773414056, 0.052333431144975136, -0.008905122105668357, 0.0007840059407332415, 0.0, 0.0, 0.0, 0.0] [-0.004286144167448658, 0.008264073453404727, -8.130035783403134e-5, -0.04144974144338818, 0.2061182017256204, -0.8114609971447189, 0.0, 0.8288973712486651, -0.23611473523103743, 0.05925148483279516, -0.009932840126144218, 0.0007946272100915355, 0.0, 0.0, 0.0] [0.008289486953575544, -0.01992892123103868, 0.023847687855928824, -0.01725455228860255, -0.031227534064274177, 0.2187728437910129, -0.8278065503298393, 0.0, 0.8329375462609859, -0.23802243145944782, 0.05951861162798522, -0.009919768604664269, 0.0007935814883731396, 0.0, 0.0] [-0.003203103055575049, 0.007743290435329053, -0.009660000293183757, 0.00943487651309554, 0.00020969357970332722, -0.05297699654295238, 0.23582352986415386, -0.8330064950072007, 0.0, 0.8333327624136899, -0.23809415379332086, 0.059523538448329215, -0.009920589741388267, 0.0007936471793110595, 0.0] [0.0004450281134593904, -0.0010681530559586649, 0.0013234597796542534, -0.0013875669869698303, 0.0005096554110994863, 0.0090146855416246, -0.0591787601986842, 0.23804354829738644, -0.8333377120321315, 0.0, 0.8333344878759047, -0.23809556796454703, 0.05952389199113576, -0.009920648665189359, 0.0007936518932151467] [0.0, 0.0, 0.0, 0.0, -0.0007936507936507937, 0.00992063492063492, -0.05952380952380952, 0.23809523809523808, -0.8333333333333334, 0.0, 0.8333333333333334, -0.23809523809523808, 0.05952380952380952, -0.00992063492063492, 0.0007936507936507937] [0.0, 0.0, 0.0, 0.0, 0.0, -0.0007936507936507937, 0.00992063492063492, -0.05952380952380952, 0.23809523809523808, -0.8333333333333334, 0.0, 0.8333333333333334, -0.23809523809523808, 0.05952380952380952, -0.00992063492063492] [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.0007936507936507937, 0.00992063492063492, -0.05952380952380952, 0.23809523809523808, -0.8333333333333334, 0.0, 0.8333333333333334, -0.23809523809523808, 0.05952380952380952] [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.0007936507936507937, 0.00992063492063492, -0.05952380952380952, 0.23809523809523808, -0.8333333333333334, 0.0, 0.8333333333333334, -0.23809523809523808] [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.0007936507936507937, 0.00992063492063492, -0.05952380952380952, 0.23809523809523808, -0.8333333333333334, 0.0, 0.8333333333333334]

Your q1 in the code should have the entries of Q0_ in the MATLAB code, i.e., q1[i] = Q0_i, correct? They use

Q0_0 = -5.0000000000000e-01;

in https://bitbucket.org/martinalmquist/optimized_sbp_operators/src/6bad240b56cddc50a4845ed8a39a907c7aae7a5e/Accurate/D1_accurate_10th.m#lines-88
but you used

q1 = [-5,

so there is a factor of ten difference?

xlxs4 · 2022-12-09T23:31:42Z

No checks error on the inner_indices, if that holds any value

ranocha · 2022-12-12T07:19:17Z

test/SBP_operators_test.jl

+        mul!(res, D, x0)
+        @test all(i->abs(res[i]) < 16000*eps(T), eachindex(res))
+        mul!(res, D, x1)
+        @test all(i->isapprox(res[i], x0[i], rtol=28000*eps(T)), eachindex(res))
+        mul!(res, D, x2)
+        @test all(i->isapprox(res[i], 2*x1[i], rtol=18000*eps(T)), eachindex(res))
+        mul!(res, D, x3)
+        @test all(i->isapprox(res[i], 3*x2[i], rtol=18000*eps(T)), eachindex(res))


Some tests in these lines fail in CI

I know, I'll go over the coefficients again to check if I've missed anything. I just pushed it just in case anyone else wanted to take a look

xlxs4 added 4 commits October 29, 2022 01:28

Add norm matrix and boundaries for order 10 accurate mattson

5d7d852

Bump minor version

044305f

Add acc order 10 grid points vector

c8a04cd

Remove commented out code

07d2004

xlxs4 added 3 commits November 2, 2022 20:42

Correct boundaries

fe44ec5

Working difference operator for order 10

3864143

Add 10th order coeffs

fb5dbf9

xlxs4 added 4 commits November 3, 2022 16:05

Add notes on 12th order

87e173d

Add xvector and coeffs for 12th order

1a64de6

Correct dims 12th order coeffs

9a1068f

Add 10th order minimal

6617d10

xlxs4 added 2 commits November 4, 2022 17:33

Add minimal 12th order notes

99e394c

Add minimal order 12th coeffs

1e0bcc2

xlxs4 added 2 commits November 7, 2022 23:52

Add norm matrix weights for accurate 10 and 12 order

40553e9

Add norm matrix weights for minimal 10 and 12 order

3db8bd2

ranocha requested changes Nov 8, 2022

View reviewed changes

ranocha reviewed Dec 5, 2022

View reviewed changes

Accurate VanDerWeide coeff fix

23d343e

Add VanDerWeide 1st derivative 10th 12th order tests

6b19bd2

Merge branch 'main' into boundary-optimized-operators

6bac4b4

ranocha reviewed Dec 12, 2022

View reviewed changes

Conversation

xlxs4 commented Nov 1, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

ranocha commented Nov 1, 2022

Uh oh!

xlxs4 commented Nov 2, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

ranocha commented Nov 2, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

ranocha commented Nov 3, 2022

Uh oh!

xlxs4 commented Nov 3, 2022

Uh oh!

xlxs4 commented Nov 3, 2022

Uh oh!

ranocha commented Nov 4, 2022

Uh oh!

xlxs4 commented Nov 4, 2022

Uh oh!

ranocha commented Nov 7, 2022

Uh oh!

xlxs4 commented Nov 7, 2022

Uh oh!

ranocha left a comment

Choose a reason for hiding this comment

Uh oh!

xlxs4 commented Nov 17, 2022

Uh oh!

xlxs4 commented Dec 4, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

ranocha Dec 5, 2022

Choose a reason for hiding this comment

Uh oh!

ranocha Dec 5, 2022

Choose a reason for hiding this comment

Uh oh!

xlxs4 Dec 8, 2022

Choose a reason for hiding this comment

Uh oh!

ranocha Dec 9, 2022

Choose a reason for hiding this comment

Uh oh!

xlxs4 commented Dec 9, 2022

Uh oh!

ranocha Dec 12, 2022

Choose a reason for hiding this comment

Uh oh!

xlxs4 Dec 12, 2022

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

xlxs4 commented Nov 1, 2022 •

edited

Loading

xlxs4 commented Nov 2, 2022 •

edited

Loading

ranocha commented Nov 2, 2022 •

edited

Loading

xlxs4 commented Dec 4, 2022 •

edited

Loading