diff --git a/dev/MattssonAlmquistVanDerWeide2018Accurate_dev.jl b/dev/MattssonAlmquistVanDerWeide2018Accurate_dev.jl index d34d6664..205aa9f1 100644 --- a/dev/MattssonAlmquistVanDerWeide2018Accurate_dev.jl +++ b/dev/MattssonAlmquistVanDerWeide2018Accurate_dev.jl @@ -152,3 +152,362 @@ q12 = [0 0 0 0 0 0 0 1//280 -4//105 1//5 -4//5 0] Q = vcat(q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12) D1 = H \ (Q - 1//2*e*e'); display(D1) + + + +# 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.0000000000000e-01, + 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) + + + +# order 12 +h = [ + 1.0000000000011e-01, + 5.9616216757547e-01, + 9.9065699844442e-01, + 1.2512548713913e+00, + 1.3316678989403e+00, + 1.2093375037721e+00, + 1.0236491595704e+00, + 9.9685258909811e-01, + 1.0004766563923e+00, + 9.9993617879146e-01, + 1.0000063122914e+00, + 9.9999966373260e-01, + 1, 1, 1, 1, 1, 1] +e = [1//1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] +H = Diagonal(h) + +q1 = [-5, + 6.7597560728423e-01, + -2.6859785384416e-01, + 1.4850302678903e-01, + -8.7976689586154e-02, + 4.1833336322613e-02, + -2.2216684976993e-03, + -1.5910034062022e-02, + 1.1296706376589e-02, + -3.1823678285130e-03, + 2.4843594063649e-04, + 3.1501105449828e-05, + 0, 0, 0, 0, 0, 0]' +q2 = [-6.7597560728423e-01, + 0, + 9.5424013647146e-01, + -4.3389334603464e-01, + 2.4285669347653e-01, + -1.1443465137214e-01, + 8.5942765682435e-03, + 4.0290424215772e-02, + -2.9396383714543e-02, + 8.5601827834256e-03, + -7.8128092862319e-04, + -6.0444181254875e-05, + 0, 0, 0, 0, 0, 0]' +q3 = [2.6859785384416e-01, + -9.5424013647146e-01, + 0, + 9.7065114311923e-01, + -4.3205328628292e-01, + 1.9549970932735e-01, + -2.4406885385172e-02, + -5.5737279079895e-02, + 4.3772338637753e-02, + -1.3727655130726e-02, + 1.6271304373071e-03, + 1.7066984372933e-05, + 0, 0, 0, 0, 0, 0]' +q4 = [-1.4850302678903e-01, + 4.3389334603464e-01, + -9.7065114311923e-01, + 0, + 9.5375878629204e-01, + -3.6113954384951e-01, + 6.9749289223875e-02, + 6.5161366516465e-02, + -6.0325702283960e-02, + 2.1188913621662e-02, + -3.2632650250470e-03, + 1.3097937809499e-04, + 0, 0, 0, 0, 0, 0]' +q5 = [8.7976689586154e-02, + -2.4285669347653e-01, + 4.3205328628292e-01, + -9.5375878629204e-01, + 0, + 8.8676146394834e-01, + -2.1292503103800e-01, + -4.6037018833218e-02, + 7.4338719466734e-02, + -3.1217656663809e-02, + 6.1239492854797e-03, + -4.5892226603067e-04, + 0, 0, 0, 0, 0, 0]' +q6 = [-4.1833336322613e-02, + 1.1443465137214e-01, + -1.9549970932735e-01, + 3.6113954384951e-01, + -8.8676146394834e-01, + 0, + 7.7461223007026e-01, + -1.0609547334165e-01, + -4.4853791547749e-02, + 3.2436468405486e-02, + -8.4387621360184e-03, + 8.5964292632428e-04, + 0, 0, 0, 0, 0, 0]' +q7 = [2.2216684976993e-03, + -8.5942765682435e-03, + 2.4406885385172e-02, + -6.9749289223875e-02, + 2.1292503103800e-01, + -7.7461223007026e-01, + 0, + 7.4758103262966e-01, + -1.5730779067906e-01, + 2.6517620342970e-02, + -4.3175367549700e-03, + 1.1092605832824e-03, + -1.8037518037522e-04, + 0, 0, 0, 0, 0]' +q8 = [1.5910034062022e-02, + -4.0290424215772e-02, + 5.5737279079895e-02, + -6.5161366516465e-02, + 4.6037018833218e-02, + 1.0609547334165e-01, + -7.4758103262966e-01, + 0, + 8.0975719267918e-01, + -2.3568822398349e-01, + 6.9373143801571e-02, + -1.6606121869177e-02, + 2.5974025974031e-03, + -1.8037518037522e-04, + 0, 0, 0, 0]' +q9 = [-1.1296706376589e-02, + 2.9396383714543e-02, + -4.3772338637753e-02, + 6.0325702283960e-02, + -7.4338719466734e-02, + 4.4853791547749e-02, + 1.5730779067906e-01, + -8.0975719267918e-01, + 0, + 8.4765775072084e-01, + -2.6369594097148e-01, + 7.8759594625702e-02, + -1.7857142857146e-02, + 2.5974025974031e-03, + -1.8037518037522e-04, + 0, 0, 0]' +q10 = [3.1823678285130e-03, + -8.5601827834256e-03, + 1.3727655130726e-02, + -2.1188913621662e-02, + 3.1217656663809e-02, + -3.2436468405486e-02, + -2.6517620342970e-02, + 2.3568822398349e-01, + -8.4765775072084e-01, + 0, + 8.5631774953989e-01, + -2.6769768119702e-01, + 7.9365079365093e-02, + -1.7857142857146e-02, + 2.5974025974031e-03, + -1.8037518037522e-04, + 0, 0]' +q11 = [-2.4843594063649e-04, + 7.8128092862319e-04, + -1.6271304373071e-03, + 3.2632650250470e-03, + -6.1239492854797e-03, + 8.4387621360184e-03, + 4.3175367549700e-03, + -6.9373143801571e-02, + 2.6369594097148e-01, + -8.5631774953989e-01, + 0, + 8.5712580212095e-01, + -2.6785714285718e-01, + 7.9365079365093e-02, + -1.7857142857146e-02, + 2.5974025974031e-03, + -1.8037518037522e-04, + 0]' +q12 = [-3.1501105449828e-05, + 6.0444181254875e-05, + -1.7066984372933e-05, + -1.3097937809499e-04, + 4.5892226603067e-04, + -8.5964292632428e-04, + -1.1092605832824e-03, + 1.6606121869177e-02, + -7.8759594625702e-02, + 2.6769768119702e-01, + -8.5712580212095e-01, + 0, + 8.5714285714289e-01, + -2.6785714285718e-01, + 7.9365079365093e-02, + -1.7857142857146e-02, + 2.5974025974031e-03, + -1.8037518037522e-04]' +q13 = [0 0 0 0 0 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] +q14 = [0 0 0 0 0 0 1//5544 -1//385 1//56 -5//63 15//56 -6//7 0 6//7 -15//56 5//63 -1//56 1//385] +q15 = [0 0 0 0 0 0 0 1//5544 -1//385 1//56 -5//63 15//56 -6//7 0 6//7 -15//56 5//63 -1//56] +q16 = [0 0 0 0 0 0 0 0 1//5544 -1//385 1//56 -5//63 15//56 -6//7 0 6//7 -15//56 5//63] +q17 = [0 0 0 0 0 0 0 0 0 1//5544 -1//385 1//56 -5//63 15//56 -6//7 0 6//7 -15//56] +q18 = [0 0 0 0 0 0 0 0 0 0 1//5544 -1//385 1//56 -5//63 15//56 -6//7 0 6//7] + +Q = vcat(q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18) + +D1 = H \ (Q - 1//2*e*e'); display(D1) diff --git a/dev/MattssonAlmquistVanDerWeide2018Minimal_dev.jl b/dev/MattssonAlmquistVanDerWeide2018Minimal_dev.jl index 39faac1f..e73533dc 100644 --- a/dev/MattssonAlmquistVanDerWeide2018Minimal_dev.jl +++ b/dev/MattssonAlmquistVanDerWeide2018Minimal_dev.jl @@ -125,3 +125,274 @@ q10 = [0 0 0 0 0 1//280 -4//105 1//5 -4//5 0] Q = vcat(q1, q2, q3, q4, q5, q6, q7, q8, q9, q10) D1 = H \ (Q - 1//2*e*e'); display(D1) + + + +# order 10 +h = [ + 1.6717213975289e-01, + 9.3675739171278e-01, + 1.3035532379753e+00, + 1.1188461804303e+00, + 9.6664345922660e-01, + 1.0083235564392e+00, + 9.9858767377362e-01, + 1.0001163606893e+00, + 1, 1, 1, 1, 1 + ] +e = [1//1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] +H = Diagonal(h) + +q1 = [-5.0000000000000e-01, + 6.7349296966214e-01, + -2.5186401896559e-01, + 8.3431385420901e-02, + 2.5480326895984e-02, + -4.5992420658252e-02, + 1.7526412909003e-02, + -2.0746552641799e-03, + 0, 0, 0, 0, 0]' +q2 = [-6.7349296966214e-01, + 0, + 9.1982892384044e-01, + -2.7262271754043e-01, + -5.0992113348238e-02, + 1.1814647281129e-01, + -4.6693123378079e-02, + 5.8255272771571e-03, + 0, 0, 0, 0, 0]' +q3 = [2.5186401896559e-01, + -9.1982892384044e-01, + 0, + 7.8566746772741e-01, + -2.4097806629929e-02, + -1.5312168858669e-01, + 6.9451518963875e-02, + -9.9345865998262e-03, + 0, 0, 0, 0, 0]' +q4 = [-8.3431385420901e-02, + 2.7262271754043e-01, + -7.8566746772741e-01, + 0, + 6.2047871210535e-01, + 1.4776775176509e-02, + -4.6889652372990e-02, + 7.3166499053672e-03, + 7.9365079365079e-04, + 0, 0, 0, 0]' +q5 = [-2.5480326895984e-02, + 5.0992113348238e-02, + 2.4097806629929e-02, + -6.2047871210535e-01, + 0, + 6.9425006383507e-01, + -1.5686345740485e-01, + 4.2609496719925e-02, + -9.9206349206349e-03, + 7.9365079365079e-04, + 0, 0, 0]' +q6 = [4.5992420658252e-02, + -1.1814647281129e-01, + 1.5312168858669e-01, + -1.4776775176509e-02, + -6.9425006383507e-01, + 0, + 8.0719535654891e-01, + -2.2953297936781e-01, + 5.9523809523809e-02, + -9.9206349206349e-03, + 7.9365079365079e-04, + 0, 0]' +q7 = [-1.7526412909003e-02, + 4.6693123378079e-02, + -6.9451518963875e-02, + 4.6889652372990e-02, + 1.5686345740485e-01, + -8.0719535654891e-01, + 0, + 8.3142546796428e-01, + -2.3809523809524e-01, + 5.9523809523809e-02, + -9.9206349206349e-03, + 7.9365079365079e-04, + 0]' +q8 = [2.0746552641799e-03, + -5.8255272771571e-03, + 9.9345865998262e-03, + -7.3166499053672e-03, + -4.2609496719925e-02, + 2.2953297936781e-01, + -8.3142546796428e-01, + 0, + 8.3333333333333e-01, + -2.3809523809524e-01, + 5.9523809523809e-02, + -9.9206349206349e-03, + 7.9365079365079e-04]' +q9 = [0 0 -1//1260 5//504 -5//84 5//21 -5//6 0 5//6 -5//21 5//84 -5//504 1//1260] +q10 = [0 0 0 -1//1260 5//504 -5//84 5//21 -5//6 0 5//6 -5//21 5//84 -5//504] +q11 = [0 0 0 0 -1//1260 5//504 -5//84 5//21 -5//6 0 5//6 -5//21 5//84] +q12 = [0 0 0 0 0 -1//1260 5//504 -5//84 5//21 -5//6 0 5//6 -5//21] +q13 = [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) + +D1 = H \ (Q - 1//2*e*e'); display(D1) + + + +# order 12 +h = [ + 1.3013597111750e-01, + 7.6146045079020e-01, + 1.1984222247012e+00, + 1.3340123109301e+00, + 1.0951811473364e+00, + 9.7569096377130e-01, + 1.0061945410831e+00, + 9.9874339446564e-01, + 1.0001702615573e+00, + 9.9998873424721e-01, + 1, 1, 1, 1, 1, 1 +] +e = [1//1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] +H = Diagonal(h) + +q1 = [-5.0000000000000e-01, + 6.7603132599815e-01, + -2.6781065957921e-01, + 1.4050310470012e-01, + -5.4072653004710e-02, + -1.1876984028213e-02, + 2.6300694680362e-02, + -9.8077210531438e-03, + 4.2848959311712e-04, + 3.0440269352791e-04, + 0, 0, 0, 0, 0, 0]' +q2 = [-6.7603132599815e-01, + 0, + 9.5204118058043e-01, + -4.1306598236120e-01, + 1.5442577883533e-01, + 2.6535212157067e-02, + -6.7869317213141e-02, + 2.6431850942376e-02, + -1.8383496124689e-03, + -6.2904733024363e-04, + 0, 0, 0, 0, 0, 0]' +q3 = [2.6781065957921e-01, + -9.5204118058043e-01, + 0, + 9.4424869445124e-01, + -3.0369922793820e-01, + -1.7036409572828e-02, + 9.7546158402857e-02, + -4.2534720340735e-02, + 5.3471186513813e-03, + 3.5890734751923e-04, + 0, 0, 0, 0, 0, 0]' +q4 = [-1.4050310470012e-01, + 4.1306598236120e-01, + -9.4424869445124e-01, + 0, + 8.1369662782755e-01, + -8.4027084126181e-02, + -1.0721180825279e-01, + 6.1098180874949e-02, + -1.2618762739267e-02, + 7.4866320589496e-04, + 0, 0, 0, 0, 0, 0]' +q5 = [5.4072653004710e-02, + -1.5442577883533e-01, + 3.0369922793820e-01, + -8.1369662782755e-01, + 0, + 6.8140317057259e-01, + -5.0090848997730e-02, + -3.2156238350691e-02, + 1.2270208460707e-02, + -8.9539078453821e-04, + -1.8037518037522e-04, + 0, 0, 0, 0, 0]' +q6 = [1.1876984028213e-02, + -2.6535212157067e-02, + 1.7036409572828e-02, + 8.4027084126181e-02, + -6.8140317057259e-01, + 0, + 7.3535220394540e-01, + -1.7565390898074e-01, + 4.5853976429252e-02, + -1.2971393808506e-02, + 2.5974025974031e-03, + -1.8037518037522e-04, + 0, 0, 0, 0]' +q7 = [-2.6300694680362e-02, + 6.7869317213141e-02, + -9.7546158402857e-02, + 1.0721180825279e-01, + 5.0090848997730e-02, + -7.3535220394540e-01, + 0, + 8.2185236816776e-01, + -2.4842386107781e-01, + 7.6038690915127e-02, + -1.7857142857146e-02, + 2.5974025974031e-03, + -1.8037518037522e-04, + 0, 0, 0]' +q8 = [9.8077210531438e-03, + -2.6431850942376e-02, + 4.2534720340735e-02, + -6.1098180874949e-02, + 3.2156238350691e-02, + 1.7565390898074e-01, + -8.2185236816776e-01, + 0, + 8.5207110387533e-01, + -2.6676625654053e-01, + 7.9365079365093e-02, + -1.7857142857146e-02, + 2.5974025974031e-03, + -1.8037518037522e-04, + 0, 0]' +q9 = [-4.2848959311712e-04, + 1.8383496124689e-03, + -5.3471186513813e-03, + 1.2618762739267e-02, + -1.2270208460707e-02, + -4.5853976429252e-02, + 2.4842386107781e-01, + -8.5207110387533e-01, + 0, + 8.5702210251244e-01, + -2.6785714285718e-01, + 7.9365079365093e-02, + -1.7857142857146e-02, + 2.5974025974031e-03, + -1.8037518037522e-04, + 0]' +q10 = [-3.0440269352791e-04, + 6.2904733024363e-04, + -3.5890734751923e-04, + -7.4866320589496e-04, + 8.9539078453821e-04, + 1.2971393808506e-02, + -7.6038690915127e-02, + 2.6676625654053e-01, + -8.5702210251244e-01, + 0, + 8.5714285714289e-01, + -2.6785714285718e-01, + 7.9365079365093e-02, + -1.7857142857146e-02, + 2.5974025974031e-03, + -1.8037518037522e-04]' +q11 = [0 0 0 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] +q12 = [0 0 0 0 1//5544 -1//385 1//56 -5//63 15//56 -6//7 0 6/7 -15//56 5//63 -1//56 1//385] +q13 = [0 0 0 0 0 1//5544 -1//385 1//56 -5//63 15//56 -6//7 0 6/7 -15//56 5//63 -1//56] +q14 = [0 0 0 0 0 0 1//5544 -1//385 1//56 -5//63 15//56 -6//7 0 6/7 -15//56 5//63] +q15 = [0 0 0 0 0 0 0 1//5544 -1//385 1//56 -5//63 15//56 -6//7 0 6/7 -15//56] +q16 = [0 0 0 0 0 0 0 0 1//5544 -1//385 1//56 -5//63 15//56 -6//7 0 6/7] +Q = vcat(q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16) + +D1 = H \ (Q - 1//2*e*e'); display(D1) diff --git a/src/SBP_coefficients/MattssonAlmquistVanDerWeide2018Accurate.jl b/src/SBP_coefficients/MattssonAlmquistVanDerWeide2018Accurate.jl index a8d7c989..8abd6c19 100644 --- a/src/SBP_coefficients/MattssonAlmquistVanDerWeide2018Accurate.jl +++ b/src/SBP_coefficients/MattssonAlmquistVanDerWeide2018Accurate.jl @@ -55,6 +55,36 @@ function construct_grid(::MattssonAlmquistVanDerWeide2018Accurate, accuracy_orde T(6.3192851303204e+00), T(7.3192851303204e+00), ) + elseif accuracy_order == 10 + xstart = SVector( + T(0.0000000000000e+00), + T(3.5902433622052e-01), + T(1.1436659188355e+00), + T(2.2144895894456e+00), + T(3.3682742337736e+00), + T(4.4309689056870e+00), + T(5.4309689056870e+00), + T(6.4309689056870e+00), + T(7.4309689056870e+00), + T(8.4309689056870e+00), + T(9.4309689056870e+00), + ) + elseif accuracy_order == 12 + xstart = SVector( + T(0.0000000000000e+00), + T(3.6098032343909e-01), + T(1.1634317168086e+00), + T(2.2975905356987e+00), + T(3.6057529790929e+00), + T(4.8918275675510e+00), + T(6.0000000000000e+00), + T(7.0000000000000e+00), + T(8.0000000000000e+00), + T(9.0000000000000e+00), + T(1.0000000000000e+01), + T(1.1000000000000e+01), + T(1.2000000000000e+01) + ) else throw(ArgumentError("Accuracy order $accuracy_order not implemented/derived.")) end @@ -281,79 +311,356 @@ function first_derivative_coefficients(source::MattssonAlmquistVanDerWeide2018Ac left_boundary_derivatives, right_boundary_derivatives, lower_coef, central_coef, upper_coef, left_weights, right_weights, mode, 1, order, source) - # elseif order == 8 - # left_boundary = ( - # # d1 - # DerivativeCoefficientRow{T,1,6}(SVector(T(), - # T(), - # T(), - # T(), - # T(), - # T() )), - # # d2 - # DerivativeCoefficientRow{T,1,6}(SVector(T(), - # T(), - # T(), - # T(), - # T(), - # T() )), - # # d3 - # DerivativeCoefficientRow{T,1,7}(SVector(T(), - # T(), - # T(), - # T(), - # T(), - # T(), - # T() )), - # # d4 - # DerivativeCoefficientRow{T,1,8}(SVector(T(), - # T(), - # T(), - # T(), - # T(), - # T(), - # T(), - # T() )), - # # d5 - # DerivativeCoefficientRow{T,1,9}(SVector(T(), - # T(), - # T(), - # T(), - # T(), - # T(), - # T(), - # T(), - # T() )), - # # d6 - # DerivativeCoefficientRow{T,1,10}(SVector(T(), - # T(), - # T(), - # T(), - # T(), - # T(), - # T(), - # T(), - # T(), - # T() )), - # ) - # right_boundary = .- left_boundary - # upper_coef = SVector(T(4//5), T(-1//5), T(4//105), T(-1//280)) - # central_coef = zero(T) - # lower_coef = -upper_coef - # left_weights = SVector( T(), - # T(), - # T(), - # T(), - # T(), - # T() ) - # right_weights = left_weights - # left_boundary_derivatives = Tuple{}() - # right_boundary_derivatives = left_boundary_derivatives + elseif order == 10 + left_boundary = ( + # d1 + DerivativeCoefficientRow{T,1,10}(SVector(T(-5), + 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) )), + # d2 + DerivativeCoefficientRow{T,1,10}(SVector(T(-1.145265719198745), + T(0), + T(1.6131685230292827), + T(-0.7195954106367713), + T(0.3651840332042441), + T(-0.12197141381091767), + T(-0.01399423475053903), + T(0.0337917497680833), + T(-0.013128541778312972), + T(0.0018110141736784769) )), + # d3 + DerivativeCoefficientRow{T,1,10}(SVector(T(0.27900977459917853), + T(-0.9945564383180177), + T(0), + T(1.004253824704819), + T(-0.4116222831605484), + T(0.13904608960593332), + T(8.487839079429473e-5), + T(-0.024930033516808243), + T(0.010097590982069497), + T(-0.0013834032874125409) )), + # d4 + DerivativeCoefficientRow{T,1,10}(SVector(T(-0.12555079634350264), + T(0.36905436758383703), + T(-0.8354024892044419), + T(0), + T(0.7960759131488463), + T(-0.2481823909255492), + T(0.035998033476551373), + T(0.01500488078728063), + T(-0.008204065636465682), + T(0.0012065471134400296) )), + # d5 + DerivativeCoefficientRow{T,1,10}(SVector(T(0.06879174149957458), + T(-0.1917465511011161), + T(0.3505621014987283), + T(-0.8150200626888124), + T(0), + T(0.7435190828161148), + T(-0.18326819871996888), + T(0.027802276682660717), + T(-0.00018667744812003666), + T(-0.0004537125390638899) )), + # d6 + DerivativeCoefficientRow{T,1,11}(SVector(T(-0.02525933759067232), + T(0.07106552895046016), + T(-0.1314044342558119), + T(0.28194859492566227), + T(-0.8250443110814595), + T(0), + T(0.8006147214131956), + T(-0.2161330773414056), + T(0.052333431144975136), + T(-0.008905122105668357), + T(0.0007840059407332415) )), + # d7 + DerivativeCoefficientRow{T,1,12}(SVector(T(-0.004286144167448658), + T(0.008264073453404727), + T(-8.130035783403134e-5), + T(-0.04144974144338818), + T(0.2061182017256204), + T(-0.8114609971447189), + T(0), + T(0.8288973712486651), + T(-0.23611473523103743), + T(0.05925148483279516), + T(-0.009932840126144218), + T(0.0007946272100915355) )), + # d8 + DerivativeCoefficientRow{T,1,13}(SVector(T(0.008289486953575544), + T(-0.01992892123103868), + T(0.023847687855928824), + T(-0.01725455228860255), + T(-0.031227534064274177), + T(0.2187728437910129), + T(-0.8278065503298393), + T(0), + T(0.8329375462609859), + T(-0.23802243145944782), + T(0.05951861162798522), + T(-0.009919768604664269), + T(0.0007935814883731396) )), + # d9 + DerivativeCoefficientRow{T,1,14}(SVector(T(-0.003203103055575049), + T(0.007743290435329053), + T(-0.009660000293183757), + T(0.00943487651309554), + T(0.00020969357970332722), + T(-0.05297699654295238), + T(0.23582352986415386), + T(-0.8330064950072007), + T(0), + T(0.8333327624136899), + T(-0.23809415379332086), + T(0.059523538448329215), + T(-0.009920589741388267), + T(0.0007936471793110595) )), + # d10 + DerivativeCoefficientRow{T,1,15}(SVector(T(0.0004450281134593904), + T(-0.0010681530559586649), + T(0.0013234597796542534), + T(-0.0013875669869698303), + T(0.0005096554110994863), + T(0.0090146855416246), + T(-0.0591787601986842), + T(0.23804354829738644), + T(-0.8333377120321315), + T(0), + T(0.8333344878759047), + T(-0.23809556796454703), + T(0.05952389199113576), + T(-0.009920648665189359), + T(0.0007936518932151467) )), + ) + right_boundary = .- left_boundary + upper_coef = SVector(T(5//6), T(-5//21), T(5//84), T(-5//504), T(1//1260)) + central_coef = zero(T) + lower_coef = -upper_coef + left_weights = SVector( T(0.1), + T(0.58980851260667), + T(0.95666820955973), + T(1.1500297411596), + T(1.1232986993248), + T(1.0123020150951), + T(0.99877122702527), + T(1.0000873322761), + T(1.0000045540888), + T(0.99999861455083) ) + right_weights = left_weights + left_boundary_derivatives = Tuple{}() + right_boundary_derivatives = left_boundary_derivatives - # DerivativeCoefficients(left_boundary, right_boundary, - # left_boundary_derivatives, right_boundary_derivatives, - # lower_coef, central_coef, upper_coef, - # left_weights, right_weights, mode, 1, order, source) + DerivativeCoefficients(left_boundary, right_boundary, + left_boundary_derivatives, right_boundary_derivatives, + lower_coef, central_coef, upper_coef, + left_weights, right_weights, mode, 1, order, source) + elseif order == 12 + left_boundary = ( + # d1 + DerivativeCoefficientRow{T,1,12}(SVector(T(-5), + T(6.759756072834865), + T(-2.6859785384386456), + T(1.4850302678886664), + T(-0.8797668958605722), + T(0.41833336322566983), + T(-0.022216684976968562), + T(-0.15910034062004497), + T(0.11296706376576574), + T(-0.031823678285095), + T(0.0024843594063621672), + T(0.00031501105449793353) )), + # d2 + DerivativeCoefficientRow{T,1,12}(SVector(T(-1.133878739795504), + T(0), + T(1.6006385315463008), + T(-0.7278109374153), + T(0.4073668318541633), + T(-0.1919522197081608), + T(0.01441600463041044), + T(0.06758299403605063), + T(-0.04930937471945772), + T(0.014358815854147503), + T(-0.0013105174583630138), + T(-0.00010138882428701446) )), + # d3 + DerivativeCoefficientRow{T,1,12}(SVector(T(0.2711310314931667), + T(-0.9632396863595133), + T(0), + T(0.9798054671227233), + T(-0.4361280311564467), + T(0.19734348986009648), + T(-0.024637069564437473), + T(-0.0562629438518241), + T(0.044185160662556816), + T(-0.013857122245420829), + T(0.0016424760940084238), + T(1.7227945090714997e-5) )), + # d4 + DerivativeCoefficientRow{T,1,12}(SVector(T(-0.11868327563344944), + T(0.3467665588803532), + T(-0.7757421491913473), + T(0), + T(0.7622418166744341), + T(-0.288621888399084), + T(0.05574347066982373), + T(0.05207681345049272), + T(-0.0482121617771465), + T(0.01693413077233541), + T(-0.002607993862528182), + T(0.00010467841611625528) )), + # d5 + DerivativeCoefficientRow{T,1,12}(SVector(T(0.06606503742874865), + T(-0.18237031445286608), + T(0.3244452213849524), + T(-0.7162136949092276), + T(0), + T(0.6659028611067349), + T(-0.15989349236956088), + T(-0.03457094585658544), + T(0.055823767717079045), + T(-0.023442523987137517), + T(0.0045987060965823005), + T(-0.00034462215871980254) )), + # d6 + DerivativeCoefficientRow{T,1,12}(SVector(T(-0.034591944922016164), + T(0.09462590138418898), + T(-0.16165851858356986), + T(0.298625935872379), + T(-0.733262187918924), + T(0), + T(0.6405260960270657), + T(-0.08773024322054244), + T(-0.03708955639583118), + T(0.02682168402477549), + T(-0.00697800416318577), + T(0.000710837895660168) )), + # d7 + DerivativeCoefficientRow{T,1,13}(SVector(T(0.002170341739577727), + T(-0.00839572473429305), + T(0.023843018046744612), + T(-0.06813788549696756), + T(0.20800586709547958), + T(-0.7567165203313853), + T(0), + T(0.7303098191800412), + T(-0.15367354059576246), + T(0.02590498912156464), + T(-0.0042177895762469655), + T(0.001083633560299047), + T(-0.00017620800905157678) )), + # d8 + DerivativeCoefficientRow{T,1,14}(SVector(T(0.015960267582207317), + T(-0.04041763512118102), + T(0.05591326108740171), + T(-0.0653671036511215), + T(0.04618237373980181), + T(0.1064304537119561), + T(-0.749941406387904), + T(0), + T(0.8123138782352943), + T(-0.23643237381439314), + T(0.06959217898439275), + T(-0.01665855318106881), + T(0.002605603502272154), + T(-0.00018094468765779322) )), + # d9 + DerivativeCoefficientRow{T,1,15}(SVector(T(-0.011291324294686406), + T(0.02938237841605001), + T(-0.043751484213130194), + T(0.060296961351895356), + T(-0.07430330232271287), + T(0.04483242188727414), + T(0.1572328446385835), + T(-0.8093714006273262), + T(0), + T(0.847253901732678), + T(-0.26357030849911345), + T(0.07872207124722692), + T(-0.017848635191088338), + T(0.0025961651187038034), + T(-0.00018028924435443554) )), + # d10 + DerivativeCoefficientRow{T,1,16}(SVector(T(0.003182570944036913), + T(-0.008560729139505266), + T(0.013728531302185184), + T(-0.021190266010048045), + T(0.031219649139547282), + T(-0.03243853867222733), + T(-0.02651931283756494), + T(0.23570326685083723), + T(-0.847711852715774), + T(0), + T(0.8563724042616903), + T(-0.2677147670769988), + T(0.07937014486366019), + T(-0.017858282594323618), + T(0.002597568377356208), + T(-0.00018038669287196362) )), + # d11 + DerivativeCoefficientRow{T,1,13}(SVector(T(-0.00024843437244633735), + T(0.0007812759969814333), + T(0.004317509501591903), + T(-0.06937270590083615), + T(0.2636942764563665), + T(-0.8563123442468437), + T(0), + T(0.8571203917272725), + T(-0.2678554520755134), + T(0.07936457839274735), + T(-0.017857030138368226), + T(0.002597386201944515), + T(-0.00018037404180170716) )), + # d12 + DerivativeCoefficientRow{T,1,18}(SVector(T(-3.150111604262639e-5), + T(6.0444201580289515e-5), + T(-1.7066990112005393e-5), + T(-0.00013097942213909975), + T(0.0004589224203513191), + T(-0.0008596432153942689), + T(-0.0011092609562906976), + T(0.016606127453276304), + T(-0.07875962110999503), + T(0.2676977712150535), + T(-0.8571260903445119), + T(0), + T(0.8571431453721869), + T(-0.2678572329288353), + T(0.07936510605299087), + T(-0.01785714886192302), + T(0.002597403470825212), + T(-0.00018037524102953333) )), + ) + right_boundary = .- left_boundary + upper_coef = SVector(T(6//7), T(-15//56), T(5//63), T(-1//56), T(1//385), T(-1//5544)) + central_coef = zero(T) + lower_coef = -upper_coef + left_weights = SVector( T(0.10000000000011), + T(0.59616216757547), + T(0.99065699844442), + T(1.2512548713913), + T(1.3316678989403), + T(1.2093375037721), + T(1.0236491595704), + T(0.99685258909811), + T(1.0004766563923), + T(0.99993617879146), + T(1.0000063122914), + T(0.9999996637326) ) + right_weights = left_weights + left_boundary_derivatives = Tuple{}() + right_boundary_derivatives = left_boundary_derivatives + + DerivativeCoefficients(left_boundary, right_boundary, + left_boundary_derivatives, right_boundary_derivatives, + lower_coef, central_coef, upper_coef, + left_weights, right_weights, mode, 1, order, source) else throw(ArgumentError("Order $order not implemented/derived.")) end diff --git a/src/SBP_coefficients/MattssonAlmquistVanDerWeide2018Minimal.jl b/src/SBP_coefficients/MattssonAlmquistVanDerWeide2018Minimal.jl index 9ff9a00b..0b12e328 100644 --- a/src/SBP_coefficients/MattssonAlmquistVanDerWeide2018Minimal.jl +++ b/src/SBP_coefficients/MattssonAlmquistVanDerWeide2018Minimal.jl @@ -100,6 +100,32 @@ function construct_grid(::MattssonAlmquistVanDerWeide2018Minimal, accuracy_order T(4.4051531374839e+00), T(5.4051531374839e+00), ) + elseif accuracy_order == 10 + xstart = SVector( + T(0.0000000000000e+00), + T(5.8556160757529e-01), + T(1.7473267488572e+00), + T(3.0000000000000e+00), + T(4.0000000000000e+00), + T(5.0000000000000e+00), + T(6.0000000000000e+00), + T(7.0000000000000e+00), + T(8.0000000000000e+00), + ) + elseif accuracy_order == 12 + xstart = SVector( + T(0.0000000000000e+00), + T(4.6552112904489e-01), + T(1.4647984306493e+00), + T(2.7620429464763e+00), + T(4.0000000000000e+00), + T(5.0000000000000e+00), + T(6.0000000000000e+00), + T(7.0000000000000e+00), + T(8.0000000000000e+00), + T(9.0000000000000e+00), + T(1.0000000000000e+01) + ) else throw(ArgumentError("Accuracy order $accuracy_order not implemented/derived.")) end @@ -272,79 +298,269 @@ function first_derivative_coefficients(source::MattssonAlmquistVanDerWeide2018Mi left_boundary_derivatives, right_boundary_derivatives, lower_coef, central_coef, upper_coef, left_weights, right_weights, mode, 1, order, source) - # elseif order == 8 - # left_boundary = ( - # # d1 - # DerivativeCoefficientRow{T,1,6}(SVector(T(), - # T(), - # T(), - # T(), - # T(), - # T() )), - # # d2 - # DerivativeCoefficientRow{T,1,6}(SVector(T(), - # T(), - # T(), - # T(), - # T(), - # T() )), - # # d3 - # DerivativeCoefficientRow{T,1,7}(SVector(T(), - # T(), - # T(), - # T(), - # T(), - # T(), - # T() )), - # # d4 - # DerivativeCoefficientRow{T,1,8}(SVector(T(), - # T(), - # T(), - # T(), - # T(), - # T(), - # T(), - # T() )), - # # d5 - # DerivativeCoefficientRow{T,1,9}(SVector(T(), - # T(), - # T(), - # T(), - # T(), - # T(), - # T(), - # T(), - # T() )), - # # d6 - # DerivativeCoefficientRow{T,1,10}(SVector(T(), - # T(), - # T(), - # T(), - # T(), - # T(), - # T(), - # T(), - # T(), - # T() )), - # ) - # right_boundary = .- left_boundary - # upper_coef = SVector(T(4//5), T(-1//5), T(4//105), T(-1//280)) - # central_coef = zero(T) - # lower_coef = -upper_coef - # left_weights = SVector( T(), - # T(), - # T(), - # T(), - # T(), - # T() ) - # right_weights = left_weights - # left_boundary_derivatives = Tuple{}() - # right_boundary_derivatives = left_boundary_derivatives + elseif order == 10 + left_boundary = ( + # d1 + DerivativeCoefficientRow{T,1,5}(SVector(T(0.4990746995535701), + T(0.15241969704789585), + T(-0.2751201290253085), + T(0.10484051310768731), + T(-0.012410293169942119) )), + # d2 + DerivativeCoefficientRow{T,1,8}(SVector(T(-0.7189620019231621), + T(0), + T(0.981928653008664), + T(-0.29102809324190426), + T(-0.054434706146276914), + T(0.1261228081640961), + T(-0.0498454816488874), + T(0.006218821787470101) )), + # d3 + DerivativeCoefficientRow{T,1,8}(SVector(T(0.19321345045852464), + T(-0.7056320348443407), + T(0), + T(0.60271222136483), + T(-0.018486246612649363), + T(-0.11746485231744051), + T(0.05327862103411152), + T(-0.0076211590830434816) )), + # d4 + DerivativeCoefficientRow{T,1,8}(SVector(T(-0.07456912923348757), + T(0.24366416251747966), + T(-0.7022122267291899), + T(0), + T(0.5545701660854921), + T(0.013207155223809193), + T(-0.04190893546685442), + T(0.006539460055673846), + T(0.0007093475470824396) )), + # d5 + DerivativeCoefficientRow{T,1,9}(SVector(T(-0.026359591690994846), + T(0.05275172853187895), + T(0.024929363975844176), + T(-0.6418899400631001), + T(0), + T(0.7182069637035886), + T(-0.1622764380264411), + T(0.04407984796588429), + T(-0.010262972170289429), + T(0.0008210377736231522) )), + # d6 + DerivativeCoefficientRow{T,1,10}(SVector(T(0.04561276027376562), + T(-0.11717119178343229), + T(0.1518576925123369), + T(-0.014654795161874228), + T(-0.6885191359475415), + T(0), + T(0.8005320825781802), + T(-0.22763821979760565), + T(0.059032449597837176), + T(-0.009838741599639595), + T(0.0007870993279711656) )), + # d7 + DerivativeCoefficientRow{T,1,10}(SVector(T(-0.01755120093038144), + T(0.04675916256969975), + T(-0.06954974589403921), + T(0.04695596952022851), + T(0.15708531311233767), + T(-0.8083369920825813), + T(0), + T(0.8326013727190911), + T(-0.2384319818363953), + T(0.059607995459097826), + T(-0.009934665909849703), + T(0.0007947732727879743) )), + # d8 + DerivativeCoefficientRow{T,1,10}(SVector(T(0.0020744138839504703), + T(-0.005824849493654949), + T(0.009933430738978299), + T(-0.0073157986339953684), + T(-0.042604539226373306), + T(0.22950627385958502), + T(-0.8313287339797591), + T(0), + T(0.8332363773741089), + T(-0.23806753639260542), + T(0.059516884098150355), + T(-0.009919480683025124), + T(0.0007935584546420079) )), + ) + right_boundary = .- left_boundary + upper_coef = SVector(T(5//6), T(-5//21), T(5//84), T(-5//504), T(1//1260)) + central_coef = zero(T) + lower_coef = -upper_coef + left_weights = SVector( T(0.16717213975289), + T(0.93675739171278), + T(1.3035532379753), + T(1.1188461804303), + T(0.9666434592266), + T(1.0083235564392), + T(0.99858767377362), + T(1.0001163606893) ) + right_weights = left_weights + left_boundary_derivatives = Tuple{}() + right_boundary_derivatives = left_boundary_derivatives - # DerivativeCoefficients(left_boundary, right_boundary, - # left_boundary_derivatives, right_boundary_derivatives, - # lower_coef, central_coef, upper_coef, - # left_weights, right_weights, mode, 1, order, source) + DerivativeCoefficients(left_boundary, right_boundary, + left_boundary_derivatives, right_boundary_derivatives, + lower_coef, central_coef, upper_coef, + left_weights, right_weights, mode, 1, order, source) + elseif order == 12 + left_boundary = ( + # d1 + DerivativeCoefficientRow{T,1,10}(SVector(T(-7.684270470438172), + T(5.194807555458745), + T(-2.0579295430730933), + T(1.0796638584520148), + T(-0.41550889074234293), + T(-0.09126595764586296), + T(0.202101651484316), + T(-0.07536518127096767), + T(0.003292629927279952), + T(0.0023391126289983594) )), + # d2 + DerivativeCoefficientRow{T,1,8}(SVector(T(-0.8878088485050056), + T(0), + T(1.250283162562752), + T(-0.5424654450963853), + T(0.2028021004571883), + T(0.034847787734123656), + T(-0.0891304560108277), + T(0.03471204698149003), + T(-0.0024142417515724767), + T(-0.0008261063717634195) )), + # d3 + DerivativeCoefficientRow{T,1,8}(SVector(T(0.22346937002606296), + T(-0.7944121537113519), + T(0), + T(0.7879098659795528), + T(0.2534158843840873), + T(-0.014215698959584674), + T(0.08139548515731004), + T(-0.035492265967731106), + T(0.004461798639218733), + T(0.00029948322062261124) )), + # d4 + DerivativeCoefficientRow{T,1,8}(SVector(T(-0.10532369420351033), + T(0.3096418068834778), + T(-0.7078260723043035), + T(0), + T(0.6099618580432923), + T(-0.06298823739309845), + T(-0.08036793017152877), + T(0.04580031261656808), + T(-0.00945925508773524), + T(0.0005612116168350629) )), + # d5 + DerivativeCoefficientRow{T,1,9}(SVector(T(0.049373250385308945), + T(-0.14100478191293775), + T(0.2773050181486685), + T(-0.7429790311918253), + T(0), + T(0.6221830719327455), + T(-0.045737501160932524), + T(-0.029361570393079246), + T(0.011203816364579948), + T(-0.0008175732267816134), + T(-0.00016469894575332322) )), + # d6 + DerivativeCoefficientRow{T,1,10}(SVector(T(0.012172895383088677), + T(-0.027196328696641285), + T(0.017460866406897767), + T(0.08612059273501356), + T(-0.6983801181665034), + T(0), + T(0.7536732748892866), + T(-0.18003027136972943), + T(0.046996413958795336), + T(-0.013294572041918046), + T(0.002662116073478288), + T(-0.00018486917176933042) )), + # d7 + DerivativeCoefficientRow{T,1,10}(SVector(T(-0.02613877695266672), + T(0.06745148621069272), + T(-0.09694562474753163), + T(0.10655176894259812), + T(0.04978246944553149), + T(-0.7308250779753221), + T(0), + T(0.8167927121559334), + T(-0.2468944631824365), + T(0.07557056593974017), + T(-0.017747207053940086), + T(0.0025814119351186032), + T(-0.00017926471771657437) )), + # d8 + DerivativeCoefficientRow{T,1,10}(SVector(T(0.009820060996139302), + T(-0.02646510714247867), + T(0.04258823695499128), + T(-0.061175053786101385), + T(0.03219669689820139), + T(0.1758749143714943), + T(-0.8228864117869611), + T(0), + T(0.85314316830222), + T(-0.2671018982641268), + T(0.07946493544275793), + T(-0.01787961047462061), + T(0.002600670614490316), + T(-0.00018060212600627668) )), + # d9 + DerivativeCoefficientRow{T,1,10}(SVector(T(-0.0004284166502310784), + T(0.0018380366654838602), + T(-0.005346208397613871), + T(0.012616614614814827), + T(-0.01226811967154658), + T(-0.0458461705888513), + T(0.24838157124468524), + T(-0.8519260536187364), + T(0), + T(0.8568762094346083), + T(-0.2678115448464915), + T(0.07935156884340752), + T(-0.017854102989766765), + T(0.0025969604348752115), + T(-0.00018034447464411663) )), + # d10 + DerivativeCoefficientRow{T,1,10}(SVector(T(-0.0003044061228920383), + T(0.0006290544170151834), + T(-0.0003589113909262337), + T(-0.0007486716402445799), + T(0.0008954008719030808), + T(0.012971539942668298), + T(-0.07603954755787204), + T(0.2667692618970865), + T(-0.8570317576203546), + T(0), + T(0.8571525136112116), + T(-0.2678601605035305), + T(0.07936597348253019), + T(-0.017857344033569367), + T(0.002597431859428317), + T(-0.00018037721246030451) )), + ) + right_boundary = .- left_boundary + upper_coef = SVector(T(6//7), T(-15//56), T(5//63), T(-1//56), T(1//385), T(-1//5544)) + central_coef = zero(T) + lower_coef = -upper_coef + left_weights = SVector( T(0.1301359711175), + T(0.7614604507902), + T(1.1984222247012), + T(1.3340123109301), + T(1.0951811473364), + T(0.9756909637713), + T(1.0061945410831), + T(0.99874339446564), + T(1.0001702615573), + T(0.99998873424721) ) + right_weights = left_weights + left_boundary_derivatives = Tuple{}() + right_boundary_derivatives = left_boundary_derivatives + + DerivativeCoefficients(left_boundary, right_boundary, + left_boundary_derivatives, right_boundary_derivatives, + lower_coef, central_coef, upper_coef, + left_weights, right_weights, mode, 1, order, source) else throw(ArgumentError("Order $order not implemented/derived.")) end diff --git a/test/SBP_operators_test.jl b/test/SBP_operators_test.jl index 3b8eb6d4..f690bd4e 100644 --- a/test/SBP_operators_test.jl +++ b/test/SBP_operators_test.jl @@ -232,6 +232,114 @@ for source in D_test_list, T in (Float32,Float64) @test derivative_right(D, x1, Val{0}()) ≈ x1[end] end + 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) + for compact in (true, false) + show(IOContext(devnull, :compact=>false), D) + show(IOContext(devnull, :compact=>false), D.coefficients) + end + 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 + @test derivative_order(D) == der_order + @test accuracy_order(D) == acc_order + @test issymmetric(D) == false + @test SummationByPartsOperators.xmin(D) ≈ xmin + @test SummationByPartsOperators.xmax(D) ≈ xmax + # interior and boundary + 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)) + mul!(res, D, x4) + @test all(i->res[i] ≈ 4*x3[i], inner_indices) + # only interior + mul!(res, D, x5) + @test all(i->res[i] ≈ 5*x4[i], inner_indices) + mul!(res, D, x6) + @test all(i->isapprox(res[i], 6*x5[i], rtol=8000*eps(T)), inner_indices) + # integration + k=0; @test integrate(x0, D) ≈ (xmax^(k+1)-(xmin)^(k+1))/(k+1) + k=1; @test integrate(x1, D) ≈ (xmax^(k+1)-(xmin)^(k+1))/(k+1) + k=2; @test integrate(x2, D) ≈ (xmax^(k+1)-(xmin)^(k+1))/(k+1) + k=3; @test integrate(x3, D) ≈ (xmax^(k+1)-(xmin)^(k+1))/(k+1) + k=4; @test integrate(x4, D) ≈ (xmax^(k+1)-(xmin)^(k+1))/(k+1) + k=5; @test integrate(x5, D) ≈ (xmax^(k+1)-(xmin)^(k+1))/(k+1) + k=6; @test integrate(x6, D) ≈ (xmax^(k+1)-(xmin)^(k+1))/(k+1) + k=7; @test integrate(x7, D) ≈ (xmax^(k+1)-(xmin)^(k+1))/(k+1) + end + + acc_order = 12 + 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) + for compact in (true, false) + show(IOContext(devnull, :compact=>false), D) + show(IOContext(devnull, :compact=>false), D.coefficients) + end + 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 + @test derivative_order(D) == der_order + @test accuracy_order(D) == acc_order + @test issymmetric(D) == false + @test SummationByPartsOperators.xmin(D) ≈ xmin + @test SummationByPartsOperators.xmax(D) ≈ xmax + # interior and boundary + 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)) + mul!(res, D, x4) + @test all(i->res[i] ≈ 4*x3[i], inner_indices) + # only interior + mul!(res, D, x5) + @test all(i->res[i] ≈ 5*x4[i], inner_indices) + mul!(res, D, x6) + @test all(i->isapprox(res[i], 6*x5[i], rtol=8000*eps(T)), inner_indices) + # integration + k=0; @test integrate(x0, D) ≈ (xmax^(k+1)-(xmin)^(k+1))/(k+1) + k=1; @test integrate(x1, D) ≈ (xmax^(k+1)-(xmin)^(k+1))/(k+1) + k=2; @test integrate(x2, D) ≈ (xmax^(k+1)-(xmin)^(k+1))/(k+1) + k=3; @test integrate(x3, D) ≈ (xmax^(k+1)-(xmin)^(k+1))/(k+1) + k=4; @test integrate(x4, D) ≈ (xmax^(k+1)-(xmin)^(k+1))/(k+1) + k=5; @test integrate(x5, D) ≈ (xmax^(k+1)-(xmin)^(k+1))/(k+1) + k=6; @test integrate(x6, D) ≈ (xmax^(k+1)-(xmin)^(k+1))/(k+1) + k=7; @test integrate(x7, D) ≈ (xmax^(k+1)-(xmin)^(k+1))/(k+1) + end + @test_throws Union{MethodError,ArgumentError} derivative_operator(source, der_order, 16, xmin, xmax, N) end