From bb92b2b319f35249e688db2f72904cbe55aab103 Mon Sep 17 00:00:00 2001 From: Utkarsh Date: Sun, 23 Feb 2020 16:39:07 +0530 Subject: [PATCH 1/4] Add Butcher Table for PFRK87 --- src/alg_utils.jl | 1 + src/algorithms.jl | 5 + src/tableaus/high_order_rk_tableaus.jl | 158 +++++++++++++++++++++++++ 3 files changed, 164 insertions(+) diff --git a/src/alg_utils.jl b/src/alg_utils.jl index 52da7dfb7e..26a6e0cf06 100644 --- a/src/alg_utils.jl +++ b/src/alg_utils.jl @@ -290,6 +290,7 @@ alg_order(alg::Hairer42) = 4 alg_order(alg::Feagin10) = 10 alg_order(alg::Feagin12) = 12 alg_order(alg::Feagin14) = 14 +alg_order(alg::PFRK87) = 8 alg_order(alg::Rosenbrock23) = 2 alg_order(alg::Rosenbrock32) = 3 diff --git a/src/algorithms.jl b/src/algorithms.jl index 64f958c0d4..c92f07514c 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -403,6 +403,11 @@ struct FRK65{T} <: OrdinaryDiffEqAdaptiveAlgorithm FRK65(omega=0.0) = new{typeof(omega)}(omega) end +struct PFRK87{T} <: OrdinaryDiffEqAdaptiveAlgorithm + omega::T + PFRK87(omega=0.0) = new{typeof(omega)}(omega) +end + ################################################################################ diff --git a/src/tableaus/high_order_rk_tableaus.jl b/src/tableaus/high_order_rk_tableaus.jl index 07b54bf5a6..8219a46efa 100644 --- a/src/tableaus/high_order_rk_tableaus.jl +++ b/src/tableaus/high_order_rk_tableaus.jl @@ -1020,3 +1020,161 @@ function DP8Interp_polyweights(T::Type) return d401,d406,d407,d408,d409,d410,d411,d412,d413,d414,d415,d416,d501,d506,d507,d508,d509,d510,d511,d512,d513,d514,d515,d516,d601,d606,d607,d608,d609,d610,d611,d612,d613,d614,d615,d616,d701,d706,d707,d708,d709,d710,d711,d712,d713,d714,d715,d716 end +struct PFRK87ConstantCache{T1,T2} <: OrdinaryDiffEqConstantCache + + α0201::T1 + α0301::T1 + α0401::T1 + α0501::T1 + α0601::T1 + α0701::T1 + + α0302::T1 + + α0403::T1 + α0503::T1 + + α0504::T1 + α0604::T1 + α0704::T1 + + α0605::T1 + α0705::T1 + + α0706::T1 + + α1108::T1 + α1208::T1 + α1308::T1 + + α1009::T1 + α1109::T1 + α1209::T1 + α1309::T1 + + α1110::T1 + α1210::T1 + α1310::T1 + + α1211::T1 + α1311::T1 + + β1::T1 + β6::T1 + β7::T1 + β8::T1 + β9::T1 + β10::T1 + β11::T1 + β12::T1 + β13::T1 + + β1tilde::T1 + β6tilde::T1 + β7tilde::T1 + β8tilde::T1 + β9tilde::T1 + β10tilde::T1 + β11tilde::T1 + β12tilde::T1 + + c2::T2 + c3::T2 + c4::T2 + c5::T2 + c6::T2 + c7::T2 + c8::T2 + c9::T2 + + d1::T1 + d2::T1 + d3::T1 + d4::T1 + d5::T1 + d6::T1 + d7::T1 + d8::T1 + d9::T1 + d10::T1 + d11::T1 + d12::T1 + d13::T1 +end + + +function PFRK87ConstantCache(T1::Type, T2::Type) + + # elements of Butcher Table + α0201 = convert(T1,1//18) + α0301 = convert(T1,1//48) + α0401 = convert(T1, 1//32) + α0501 = convert(T1, 5//16) + α0601 = convert(T1, 3//80) + α0701 = convert(T1, 29443841//614563906) + + + α0302 = convert(T1, 1//16) + + α0403 = convert(T1, 3//32) + α0503 = convert(T1, -75//64) + + α0504 = convert(T1, 75//64) + α0604 = convert(T1, 3//16) + α0704 = convert(T1, 77736538//692538347) + + α0605 = convert(T1, 3//20) + α0705 = convert(T1, -28693883//1125000000) + + α0706 = convert(T1, 23124283//1800000000) + + α1108 = convert(T1, 15336726248//1032824649) + α1208 = convert(T1, 5232866602//850066563) + α1308 = convert(T1, -13158990841//6184727034) + + α1009 = convert(T1, 123872331//1001029789) + α1109 = convert(T1, -45442868181//3398467696) + α1209 = convert(T1, -4093664535//808688257) + α1309 = convert(T1, 3936647629//1978049680) + + α1110 = convert(T1, 3065993473//597172653) + α1210 = convert(T1, 3962137247//1805957418) + α1310 = convert(T1, -160528059//685178525) + + α1211 = convert(T1, 65686358//487910083) + α1311 = convert(T1, 248638103//1413531060) + + β1 = convert(T1, 14005451//335480064) + β6 = convert(T1, -59238493//1068277825) + β7 = convert(T1, 181606767//758867731) + β8 = convert(T1, 561292985//797845732) + β9 = convert(T1, -1041891430//1371343529) + β10 = convert(T1, 760417239//1151165299) + β11 = convert(T1, 118820643//751138087) + β12 = convert(T1, -528747749//220607170) + β13 = convert(T1, 1//4) + + β1tilde = convert(T1, 13451932//455176623) + β6tilde = convert(T1, -808719846//976000145) + β7tilde = convert(T1, 1757004468//5645159321) + β8tilde = convert(T1, 656045339//265891186) + β9tilde = convert(T1, -3867574721//1518517206) + β10tilde = convert(T1, 465885868//322736535) + β11tilde = convert(T1, 53011238//667516719) + β12tilde = convert(T1, 2//45) + + c2 = convert(T2, 1//18) + c3 = convert(T2, 1//12) + c4 = convert(T2, 1//8) + c5 = convert(T2, 5//16) + c6 = convert(T2, 3//8) + c7 = convert(T2, 59//400) + c8 = convert(T2, 93//200) + c9 = convert(T2, 5490023248//9719169821) + c10 = convert(T2, 5490023248//9719169821) + c11 = convert(T2, 1201146811//1299019798) + c12 = convert(T2, 1//1) + c13 = convert(T2, 1//1) + + PFRK87ConstantCache(α0201, α0301, α0401, α0501, α0601, α0701, α0801, α0901, α0302, α0403, α0503, α63, α73, α83, α54, α64, α74, α84, α94, α65, α75, α85, α95, α76, α86, α96, α87, α97, α98, β1, β7, β8, β1tilde, β4tilde, β5tilde, β6tilde, β7tilde, β8tilde, β9tilde, c2, c3, c4, c5, c6, c7, c8, c9, d1, d2, d3, d4, d5, d6, d7, d8, d9, d10, d11, d12, d13, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, e11, f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11) +end \ No newline at end of file From 9e15224d52a8a2ae7b95cfa0df20bf9249e8579d Mon Sep 17 00:00:00 2001 From: Utkarsh Date: Sun, 23 Feb 2020 21:35:21 +0530 Subject: [PATCH 2/4] Add perform step for PFRK87 --- src/OrdinaryDiffEq.jl | 2 +- src/caches/high_order_rk_caches.jl | 47 +++++ .../high_order_rk_perform_step.jl | 179 ++++++++++++++++++ src/tableaus/high_order_rk_tableaus.jl | 2 +- 4 files changed, 228 insertions(+), 2 deletions(-) diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl index 4521e6a117..b7ba5f8021 100644 --- a/src/OrdinaryDiffEq.jl +++ b/src/OrdinaryDiffEq.jl @@ -176,7 +176,7 @@ module OrdinaryDiffEq export FunctionMap, Euler, Heun, Ralston, Midpoint, RK4, ExplicitRK, OwrenZen3, OwrenZen4, OwrenZen5, BS3, BS5, RK46NL, DP5, Tsit5, DP8, Vern6, Vern7, Vern8, TanYam7, TsitPap8, - Vern9, Feagin10, Feagin12, Feagin14, CompositeAlgorithm, Anas5, RKO65, FRK65 + Vern9, Feagin10, Feagin12, Feagin14, CompositeAlgorithm, Anas5, RKO65, FRK65, PFRK87 export SSPRK22, SSPRK33, SSPRK53, SSPRK53_2N1, SSPRK53_2N2, SSPRK63, SSPRK73, SSPRK83, SSPRK432, SSPRKMSVS32, SSPRKMSVS43, SSPRK932, SSPRK54, SSPRK104 diff --git a/src/caches/high_order_rk_caches.jl b/src/caches/high_order_rk_caches.jl index 1cad1ba6bf..0910d5d8c1 100644 --- a/src/caches/high_order_rk_caches.jl +++ b/src/caches/high_order_rk_caches.jl @@ -123,3 +123,50 @@ function alg_cache(alg::TsitPap8,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNo end alg_cache(alg::TsitPap8,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false}) = TsitPap8ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits)) + +@cache struct PFRK87Cache{uType,rateType,uNoUnitsType,TabType} <: OrdinaryDiffEqMutableCache + u::uType + uprev::uType + fsalfirst::rateType + k2::rateType + k3::rateType + k4::rateType + k5::rateType + k6::rateType + k7::rateType + k8::rateType + k9::rateType + k10::rateType + k11::rateType + k12::rateType + k13::rateType + utilde::uType + tmp::uType + atmp::uNoUnitsType + k::rateType + tab::TabType +end + +function alg_cache(alg::PFRK87,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true}) + tab = PFRK87ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits)) + k1 = zero(rate_prototype) + k2 = zero(rate_prototype) + k3 = zero(rate_prototype) + k4 = zero(rate_prototype) + k5 = zero(rate_prototype) + k6 = zero(rate_prototype) + k7 = zero(rate_prototype) + k8 = zero(rate_prototype) + k9 = zero(rate_prototype) + k10 = zero(rate_prototype) + k11 = zero(rate_prototype) + k12 = zero(rate_prototype) + k13 = zero(rate_prototype) + utilde = similar(u) + k = zero(rate_prototype) + tmp = similar(u) + atmp = similar(u,uEltypeNoUnits) + PFRK87Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,utilde,tmp,atmp,k,tab) +end + +alg_cache(alg::PFRK87,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false}) = PFRK87ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits)) diff --git a/src/perform_step/high_order_rk_perform_step.jl b/src/perform_step/high_order_rk_perform_step.jl index 094707b9d5..1764c3f25c 100644 --- a/src/perform_step/high_order_rk_perform_step.jl +++ b/src/perform_step/high_order_rk_perform_step.jl @@ -601,3 +601,182 @@ end integrator.destats.nf += 1 end =# + +function initialize!(integrator, cache::PFRK87ConstantCache) + integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal + integrator.destats.nf += 1 + integrator.kshortsize = 2 + integrator.k = typeof(integrator.k)(undef, integrator.kshortsize) + + # Avoid undefined entries if k is an array of arrays + integrator.fsallast = zero(integrator.fsalfirst) + integrator.k[1] = integrator.fsalfirst + integrator.k[2] = integrator.fsallast +end + +@muladd function perform_step!(integrator, cache::PFRK87ConstantCache, repeat_step=false) + @unpack t,dt,uprev,u,f,p = integrator + @unpack α0201, α0301, α0401, α0501, α0601, α0701, α0302, α0403, α0503, α0504, α0604, α0704, α0605, α0705, α0706, α1108, α1208, α1308, α1009, α1109, α1209, α1309, α1110, α1210, α1310, α1211, α1311, β1, β6, β7, β8, β9, β10, β11, β12, β13, β1tilde, β6tilde, β7tilde, β8tilde, β9tilde, β10tilde, β11tilde, β12tilde, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12, c13 = cache + alg = unwrap_alg(integrator, true) + ν = alg.omega*dt + νsq = ν^2 + + c_ν = -0.19781108078634084 - (0.164050909125528499*νsq) + (0.042578310088756321*(νsq^2)) - (0.002300513610963998*(νsq^3)) + (0.000033467244551879287*(νsq^4)) - (7.8661142036921924*(10^(-8))*(νsq^5)) + d_ν = 1-(0.296457092123567400*νsq) + (0.0015793885907465726*(νsq^2)) - (0.00018913011771688527*(νsq^3)) + (0.000017089234650765179*(νsq^4)) - (1.2705211682518626*(10^(-7))*(νsq^5)) + + α0807 = c_ν/d_ν + α0801 = 0.026876256 + 0.0576576*α0807 + α0804 = 0.22464336 - 0.944944*α0807 + α0805 = 0.000369024 - 0.2061696*α0807 + α0806 = 0.21311136 + 0.093456*α0807 + α0901 = 0.07239997637512857 + 0.01913119863380767*α0807 + α0904 = -0.688400520601143 - 0.3135390887207368*α0807 + α0905 = -0.688400520601143 - 0.3135390887207368*α0807 + α0906 = -0.17301267570583073 - 0.06840852844816077*α0807 + α0907 = 0.1440060555560846 + 0.031009360422930017*α0807 + α0908 = 0.9982362892760762 + 0.33180705811215994*α0807 + α1001 = 0.16261514523236525 - 0.12125171966747463*α0807 + α1004 = -2.1255544052061124 + 1.9871809612169453*α0807 + α1005 = -0.216403903283323 + 0.43356675517460624*α0807 + + α1006 = -0.060417230254934076 - 0.1965343807796979*α0807 + α1007 = -0.060417230254934076 - 0.1965343807796979*α0807 + α1008 = 2.4846281621788395 - 2.102961615944379*α0807 + α1101 = -1.0320124180911034 + 1.061943768952537*α0807 + α1104 = 13.666683232895137 - 17.40407843561103*α0807 + α1105 = 0.25990355211486116 - 3.797253476860588*α0807 + α1106 = -5.759316475814002 + 1.7212824826428488*α0807 + α1107 = -12.822511612651839 + 18.41810566087623*α0807 + α1201 = 0.2478349764611783 - 0.06383934946543009*α0807 + α1204 = -4.593782880309185 + 1.046256005127882*α0807 + α1205 = -0.39566692537411896 + 0.22827403748244698*α0807 + α1206 = -3.0673550479691665 - 0.10347586863902129*α0807 + α1207 = 5.386688702227177 - 1.1072148245058775*α0807 + α1301 = 0.7332242174431163 - 0.5164807626867616*α0807 + α1304 = -10.196728938160977 + 8.464545832921925*α0807 + α1305 = -0.43865244706547707 + 1.846809999910238*α0807 + α1306 = 0.5693856884667226 - 0.8371528845746959*α0807 + α1307 = 10.52865228002416 - 8.957722185570706*α0807 + + k1 = integrator.fsalfirst + k2 = f(uprev+dt*α0201*k1, p, t + c2*dt) + k3 = f(uprev+dt*(α0301*k1+α0302*k2), p, t + c3*dt) + k4 = f(uprev+dt*(α0401*k1 +α0403*k3), p, t + c4*dt) + k5 = f(uprev+dt*(α0501*k1 +α0503*k3+α0504*k4), p, t + c5*dt) + k6 = f(uprev+dt*(α0601*k1 +α0604*k4+α0605*k5), p, t + c6*dt) + k7 = f(uprev+dt*(α0701*k1 +α0704*k4+α0705*k5+α0706*k6), p, t + c7*dt) + k8 = f(uprev+dt*(α0801*k1 +α0804*k4+α0805*k5+α0806*k6+α0807*k7), p, t + c8*dt) + k9 = f(uprev+dt*(α0901*k1 +α0904*k4+α0905*k5+α0906*k6+α0907*k7+α0908*k8), p, t + c9*dt) + k10 =f(uprev+dt*(α1001*k1 +α1004*k4+α1005*k5+α1006*k6+α1007*k7+α1008*k8+α1009*k9), p, t + c10*dt) + k11= f(uprev+dt*(α1101*k1 +α1104*k4+α1105*k5+α1106*k6+α1107*k7+α1108*k8+α1109*k9+α1110*k10), p, t + c11*dt) + k12= f(uprev+dt*(α1201*k1 +α1204*k4+α1205*k5+α1206*k6+α1207*k7+α1208*k8+α1209*k9+α1210*k10+α1211*k11), p, t+c12*dt) + k13= f(uprev+dt*(α1301*k1 +α1304*k4+α1305*k5+α1306*k6+α1307*k7+α1308*k8+α1309*k9+α1310*k10), p, t+c13*dt) + integrator.destats.nf += 12 + u = uprev + dt*(β1*k1+β6*k6+β7*k7+β8*k8+β9*k9+β10*k10+β11*k11+β12*k12+β13*k12) + if integrator.opts.adaptive + utilde = dt*(β1tilde*k1 + β6tilde*k6 + β7tilde*k7 + β8tilde*k8 + β9tilde*k9 + β10tilde*k10 + β11tilde*k11 + β12tilde*k12) + atmp = calculate_residuals(utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t) + integrator.EEst = integrator.opts.internalnorm(atmp,t) + end + integrator.fsallast = f(u, p, t+dt) + integrator.destats.nf += 1 + integrator.k[1] = integrator.fsalfirst + integrator.k[2] = integrator.fsallast + integrator.u = u +end + + +function initialize!(integrator, cache::PFRK87Cache) + integrator.fsalfirst = cache.fsalfirst + integrator.fsallast = cache.k + integrator.kshortsize = 2 + resize!(integrator.k, integrator.kshortsize) + integrator.k[1] = integrator.fsalfirst + integrator.k[2] = integrator.fsallast + integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) # Pre-start fsal + integrator.destats.nf += 1 +end + +@muladd function perform_step!(integrator, cache::PFRK87Cache, repeat_step=false) + @unpack t,dt,uprev,u,f,p = integrator + @unpack α0201, α0301, α0401, α0501, α0601, α0701, α0302, α0403, α0503, α0504, α0604, α0704, α0605, α0705, α0706, α1108, α1208, α1308, α1009, α1109, α1209, α1309, α1110, α1210, α1310, α1211, α1311, β1, β6, β7, β8, β9, β10, β11, β12, β13, β1tilde, β6tilde, β7tilde, β8tilde, β9tilde, β10tilde, β11tilde, β12tilde, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12, c13 = cache.tab + @unpack k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,utilde,tmp,atmp,k = cache + + alg = unwrap_alg(integrator, true) + ν = alg.omega*dt + νsq = ν^2 + + c_ν = -0.19781108078634084 - (0.164050909125528499*νsq) + (0.042578310088756321*(νsq^2)) - (0.002300513610963998*(νsq^3)) + (0.000033467244551879287*(νsq^4)) - (7.8661142036921924*(10^(-8))*(νsq^5)) + d_ν = 1-(0.296457092123567400*νsq) + (0.0015793885907465726*(νsq^2)) - (0.00018913011771688527*(νsq^3)) + (0.000017089234650765179*(νsq^4)) - (1.2705211682518626*(10^(-7))*(νsq^5)) + + α0807 = c_ν/d_ν + α0801 = 0.026876256 + 0.0576576*α0807 + α0804 = 0.22464336 - 0.944944*α0807 + α0805 = 0.000369024 - 0.2061696*α0807 + α0806 = 0.21311136 + 0.093456*α0807 + α0901 = 0.07239997637512857 + 0.01913119863380767*α0807 + α0904 = -0.688400520601143 - 0.3135390887207368*α0807 + α0905 = -0.688400520601143 - 0.3135390887207368*α0807 + α0906 = -0.17301267570583073 - 0.06840852844816077*α0807 + α0907 = 0.1440060555560846 + 0.031009360422930017*α0807 + α0908 = 0.9982362892760762 + 0.33180705811215994*α0807 + α1001 = 0.16261514523236525 - 0.12125171966747463*α0807 + α1004 = -2.1255544052061124 + 1.9871809612169453*α0807 + α1005 = -0.216403903283323 + 0.43356675517460624*α0807 + + α1006 = -0.060417230254934076 - 0.1965343807796979*α0807 + α1007 = -0.060417230254934076 - 0.1965343807796979*α0807 + α1008 = 2.4846281621788395 - 2.102961615944379*α0807 + α1101 = -1.0320124180911034 + 1.061943768952537*α0807 + α1104 = 13.666683232895137 - 17.40407843561103*α0807 + α1105 = 0.25990355211486116 - 3.797253476860588*α0807 + α1106 = -5.759316475814002 + 1.7212824826428488*α0807 + α1107 = -12.822511612651839 + 18.41810566087623*α0807 + α1201 = 0.2478349764611783 - 0.06383934946543009*α0807 + α1204 = -4.593782880309185 + 1.046256005127882*α0807 + α1205 = -0.39566692537411896 + 0.22827403748244698*α0807 + α1206 = -3.0673550479691665 - 0.10347586863902129*α0807 + α1207 = 5.386688702227177 - 1.1072148245058775*α0807 + α1301 = 0.7332242174431163 - 0.5164807626867616*α0807 + α1304 = -10.196728938160977 + 8.464545832921925*α0807 + α1305 = -0.43865244706547707 + 1.846809999910238*α0807 + α1306 = 0.5693856884667226 - 0.8371528845746959*α0807 + α1307 = 10.52865228002416 - 8.957722185570706*α0807 + + k1 = cache.fsalfirst + f(k1, uprev, p, t) + @.. tmp = uprev+dt*α0201*k1 + f(k2, tmp, p, t + c2*dt) + @.. tmp = uprev+dt*(α0301*k1+α0302*k2) + f(k3, tmp, p, t + c3*dt) + @.. tmp = uprev+dt*(α0401*k1+α0403*k3) + f(k4, tmp, p, t + c4*dt) + @.. tmp = uprev+dt*(α0501*k1+α0503*k3+α0504*k4) + f(k5, tmp, p, t + c5*dt) + @.. tmp = uprev+dt*(α0601*k1+α0604*k4+α0605*k5) + f(k6, tmp, p, t + c6*dt) + @.. tmp = uprev+dt*(α0701*k1+α0704*k4+α0705*k5+α0706*k6) + f(k7, tmp, p, t + c7*dt) + @.. tmp = uprev+dt*(α0801*k1+α0804*k4+α0805*k5+α0806*k6+α0807*k7) + f(k8, tmp, p, t + c8*dt) + @.. tmp = uprev+dt*(α0901*k1+α0904*k4+α0905*k5+α0906*k6+α0907*k7+α0908*k8) + f(k9, tmp, p, t + c9*dt) + @.. tmp = uprev+dt*(α1001*k1+α1004*k4+α1005*k5+α1006*k6+α1007*k7+α1008*k8+α1009*k9) + f(k10, tmp, p, t + c10*dt) + @.. tmp = uprev+dt*(α1101*k1+α1104*k4+α1105*k5+α1106*k6+α1107*k7+α1108*k8+α1109*k9+α1110*k10) + f(k11, tmp, p, t + c11*dt) + @.. tmp = uprev+dt*(α1201*k1+α1204*k4+α1205*k5+α1206*k6+α1207*k7+α1208*k8+α1209*k9+α1210*k10+α1211*k11) + f(k12, tmp, p, t+c12*dt) + @.. tmp = uprev+dt*(α1301*k1+α1304*k4+α1305*k5+α1306*k6+α1307*k7+α1308*k8+α1309*k9+α1310*k10) + f(k13, tmp, p, t+c13*dt) + @.. u = uprev + dt*(β1*k1+β6*k6+β7*k7+β8*k8+β9*k9+β10*k10+β11*k11+β12*k12++β13*k13) + integrator.destats.nf += 13 + if integrator.opts.adaptive + @.. utilde = dt*(β1tilde*k1 + β6tilde*k6 + β7tilde*k7 + β8tilde*k8 + β9tilde*k9 + β10tilde*k10 + β11tilde*k11 + βtilde12*k12) + calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t) + integrator.EEst = integrator.opts.internalnorm(atmp,t) + end + f(k, u, p, t+dt) + integrator.destats.nf += 1 + return nothing +end \ No newline at end of file diff --git a/src/tableaus/high_order_rk_tableaus.jl b/src/tableaus/high_order_rk_tableaus.jl index 8219a46efa..b76d444319 100644 --- a/src/tableaus/high_order_rk_tableaus.jl +++ b/src/tableaus/high_order_rk_tableaus.jl @@ -1176,5 +1176,5 @@ function PFRK87ConstantCache(T1::Type, T2::Type) c12 = convert(T2, 1//1) c13 = convert(T2, 1//1) - PFRK87ConstantCache(α0201, α0301, α0401, α0501, α0601, α0701, α0801, α0901, α0302, α0403, α0503, α63, α73, α83, α54, α64, α74, α84, α94, α65, α75, α85, α95, α76, α86, α96, α87, α97, α98, β1, β7, β8, β1tilde, β4tilde, β5tilde, β6tilde, β7tilde, β8tilde, β9tilde, c2, c3, c4, c5, c6, c7, c8, c9, d1, d2, d3, d4, d5, d6, d7, d8, d9, d10, d11, d12, d13, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, e11, f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11) + PFRK87ConstantCache(α0201, α0301, α0401, α0501, α0601, α0701, α0302, α0403, α0503, α0504, α0604, α0704, α0605, α0705, α0706, α1108, α1208, α1308, α1009, α1109, α1209, α1309, α1110, α1210, α1310, α1211, α1311, β1, β6, β7, β8, β9, β10, β11, β12, β13, β1tilde, β6tilde, β7tilde, β8tilde, β9tilde, β10tilde, β11tilde, β12tilde, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12, c13) end \ No newline at end of file From 01e1cbf919a511b816ee3ce8f2d02e2243803303 Mon Sep 17 00:00:00 2001 From: Utkarsh Date: Sun, 23 Feb 2020 21:35:38 +0530 Subject: [PATCH 3/4] Minor Fixes --- .../high_order_rk_perform_step.jl | 15 ++++++------- src/tableaus/high_order_rk_tableaus.jl | 22 +++++-------------- 2 files changed, 13 insertions(+), 24 deletions(-) diff --git a/src/perform_step/high_order_rk_perform_step.jl b/src/perform_step/high_order_rk_perform_step.jl index 1764c3f25c..2aa93a1d35 100644 --- a/src/perform_step/high_order_rk_perform_step.jl +++ b/src/perform_step/high_order_rk_perform_step.jl @@ -621,8 +621,8 @@ end ν = alg.omega*dt νsq = ν^2 - c_ν = -0.19781108078634084 - (0.164050909125528499*νsq) + (0.042578310088756321*(νsq^2)) - (0.002300513610963998*(νsq^3)) + (0.000033467244551879287*(νsq^4)) - (7.8661142036921924*(10^(-8))*(νsq^5)) - d_ν = 1-(0.296457092123567400*νsq) + (0.0015793885907465726*(νsq^2)) - (0.00018913011771688527*(νsq^3)) + (0.000017089234650765179*(νsq^4)) - (1.2705211682518626*(10^(-7))*(νsq^5)) + c_ν = -0.19781108078634084 - (0.164050909125528499*νsq) + (0.042578310088756321*(νsq^2)) - (0.002300513610963998*(νsq^3)) + (0.000033467244551879287*(νsq^4)) - (7.8661142036921924*(1/100000000)*(νsq^5)) + d_ν = 1-(0.296457092123567400*νsq) + (0.0015793885907465726*(νsq^2)) - (0.00018913011771688527*(νsq^3)) + (0.000017089234650765179*(νsq^4)) - (1.2705211682518626*(1/10000000)*(νsq^5)) α0807 = c_ν/d_ν α0801 = 0.026876256 + 0.0576576*α0807 @@ -638,7 +638,6 @@ end α1001 = 0.16261514523236525 - 0.12125171966747463*α0807 α1004 = -2.1255544052061124 + 1.9871809612169453*α0807 α1005 = -0.216403903283323 + 0.43356675517460624*α0807 - α1006 = -0.060417230254934076 - 0.1965343807796979*α0807 α1007 = -0.060417230254934076 - 0.1965343807796979*α0807 α1008 = 2.4846281621788395 - 2.102961615944379*α0807 @@ -670,9 +669,9 @@ end k10 =f(uprev+dt*(α1001*k1 +α1004*k4+α1005*k5+α1006*k6+α1007*k7+α1008*k8+α1009*k9), p, t + c10*dt) k11= f(uprev+dt*(α1101*k1 +α1104*k4+α1105*k5+α1106*k6+α1107*k7+α1108*k8+α1109*k9+α1110*k10), p, t + c11*dt) k12= f(uprev+dt*(α1201*k1 +α1204*k4+α1205*k5+α1206*k6+α1207*k7+α1208*k8+α1209*k9+α1210*k10+α1211*k11), p, t+c12*dt) - k13= f(uprev+dt*(α1301*k1 +α1304*k4+α1305*k5+α1306*k6+α1307*k7+α1308*k8+α1309*k9+α1310*k10), p, t+c13*dt) + k13= f(uprev+dt*(α1301*k1 +α1304*k4+α1305*k5+α1306*k6+α1307*k7+α1308*k8+α1309*k9+α1310*k10+α1311*k11), p, t+c13*dt) integrator.destats.nf += 12 - u = uprev + dt*(β1*k1+β6*k6+β7*k7+β8*k8+β9*k9+β10*k10+β11*k11+β12*k12+β13*k12) + u = uprev + dt*(β1*k1+β6*k6+β7*k7+β8*k8+β9*k9+β10*k10+β11*k11+β12*k12+β13*k13) if integrator.opts.adaptive utilde = dt*(β1tilde*k1 + β6tilde*k6 + β7tilde*k7 + β8tilde*k8 + β9tilde*k9 + β10tilde*k10 + β11tilde*k11 + β12tilde*k12) atmp = calculate_residuals(utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t) @@ -767,12 +766,12 @@ end f(k11, tmp, p, t + c11*dt) @.. tmp = uprev+dt*(α1201*k1+α1204*k4+α1205*k5+α1206*k6+α1207*k7+α1208*k8+α1209*k9+α1210*k10+α1211*k11) f(k12, tmp, p, t+c12*dt) - @.. tmp = uprev+dt*(α1301*k1+α1304*k4+α1305*k5+α1306*k6+α1307*k7+α1308*k8+α1309*k9+α1310*k10) + @.. tmp = uprev+dt*(α1301*k1+α1304*k4+α1305*k5+α1306*k6+α1307*k7+α1308*k8+α1309*k9+α1310*k10+α1311*k11) f(k13, tmp, p, t+c13*dt) - @.. u = uprev + dt*(β1*k1+β6*k6+β7*k7+β8*k8+β9*k9+β10*k10+β11*k11+β12*k12++β13*k13) + @.. u = uprev + dt*(β1*k1+β6*k6+β7*k7+β8*k8+β9*k9+β10*k10+β11*k11+β12*k12+β13*k13) integrator.destats.nf += 13 if integrator.opts.adaptive - @.. utilde = dt*(β1tilde*k1 + β6tilde*k6 + β7tilde*k7 + β8tilde*k8 + β9tilde*k9 + β10tilde*k10 + β11tilde*k11 + βtilde12*k12) + @.. utilde = dt*(β1tilde*k1 + β6tilde*k6 + β7tilde*k7 + β8tilde*k8 + β9tilde*k9 + β10tilde*k10 + β11tilde*k11 + β12tilde*k12) calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t) integrator.EEst = integrator.opts.internalnorm(atmp,t) end diff --git a/src/tableaus/high_order_rk_tableaus.jl b/src/tableaus/high_order_rk_tableaus.jl index b76d444319..8b6e16e2b2 100644 --- a/src/tableaus/high_order_rk_tableaus.jl +++ b/src/tableaus/high_order_rk_tableaus.jl @@ -1086,20 +1086,10 @@ struct PFRK87ConstantCache{T1,T2} <: OrdinaryDiffEqConstantCache c7::T2 c8::T2 c9::T2 - - d1::T1 - d2::T1 - d3::T1 - d4::T1 - d5::T1 - d6::T1 - d7::T1 - d8::T1 - d9::T1 - d10::T1 - d11::T1 - d12::T1 - d13::T1 + c10::T2 + c11::T2 + c12::T2 + c13::T2 end @@ -1151,7 +1141,7 @@ function PFRK87ConstantCache(T1::Type, T2::Type) β9 = convert(T1, -1041891430//1371343529) β10 = convert(T1, 760417239//1151165299) β11 = convert(T1, 118820643//751138087) - β12 = convert(T1, -528747749//220607170) + β12 = convert(T1, -528747749//2220607170) β13 = convert(T1, 1//4) β1tilde = convert(T1, 13451932//455176623) @@ -1171,7 +1161,7 @@ function PFRK87ConstantCache(T1::Type, T2::Type) c7 = convert(T2, 59//400) c8 = convert(T2, 93//200) c9 = convert(T2, 5490023248//9719169821) - c10 = convert(T2, 5490023248//9719169821) + c10 = convert(T2, 13//20) c11 = convert(T2, 1201146811//1299019798) c12 = convert(T2, 1//1) c13 = convert(T2, 1//1) From 1f870c84c4d57bce1eaddfbe7374f695bb0d380e Mon Sep 17 00:00:00 2001 From: Utkarsh Date: Wed, 26 Feb 2020 00:46:40 +0530 Subject: [PATCH 4/4] Fix Tableau and convergence tests passing --- .../high_order_rk_perform_step.jl | 25 ++++++++----------- src/tableaus/high_order_rk_tableaus.jl | 6 ++++- test/algconvergence/ode_convergence_tests.jl | 4 +++ 3 files changed, 19 insertions(+), 16 deletions(-) diff --git a/src/perform_step/high_order_rk_perform_step.jl b/src/perform_step/high_order_rk_perform_step.jl index 2aa93a1d35..1d0f7bf9a2 100644 --- a/src/perform_step/high_order_rk_perform_step.jl +++ b/src/perform_step/high_order_rk_perform_step.jl @@ -616,7 +616,7 @@ end @muladd function perform_step!(integrator, cache::PFRK87ConstantCache, repeat_step=false) @unpack t,dt,uprev,u,f,p = integrator - @unpack α0201, α0301, α0401, α0501, α0601, α0701, α0302, α0403, α0503, α0504, α0604, α0704, α0605, α0705, α0706, α1108, α1208, α1308, α1009, α1109, α1209, α1309, α1110, α1210, α1310, α1211, α1311, β1, β6, β7, β8, β9, β10, β11, β12, β13, β1tilde, β6tilde, β7tilde, β8tilde, β9tilde, β10tilde, β11tilde, β12tilde, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12, c13 = cache + @unpack α0201, α0301, α0401, α0501, α0601, α0701, α0302, α0403, α0503, α0504, α0604, α0704, α0605, α0705, α0706, α0908, α1008, α1108, α1208, α1308, α1009, α1109, α1209, α1309, α1110, α1210, α1310, α1211, α1311, β1, β6, β7, β8, β9, β10, β11, β12, β13, β1tilde, β6tilde, β7tilde, β8tilde, β9tilde, β10tilde, β11tilde, β12tilde, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12, c13 = cache alg = unwrap_alg(integrator, true) ν = alg.omega*dt νsq = ν^2 @@ -631,16 +631,14 @@ end α0806 = 0.21311136 + 0.093456*α0807 α0901 = 0.07239997637512857 + 0.01913119863380767*α0807 α0904 = -0.688400520601143 - 0.3135390887207368*α0807 - α0905 = -0.688400520601143 - 0.3135390887207368*α0807 - α0906 = -0.17301267570583073 - 0.06840852844816077*α0807 - α0907 = 0.1440060555560846 + 0.031009360422930017*α0807 - α0908 = 0.9982362892760762 + 0.33180705811215994*α0807 + α0905 = -0.17301267570583073 - 0.06840852844816077*α0807 + α0906 = 0.1440060555560846 + 0.031009360422930017*α0807 + α0907 = 0.9982362892760762 + 0.33180705811215994*α0807 α1001 = 0.16261514523236525 - 0.12125171966747463*α0807 α1004 = -2.1255544052061124 + 1.9871809612169453*α0807 α1005 = -0.216403903283323 + 0.43356675517460624*α0807 α1006 = -0.060417230254934076 - 0.1965343807796979*α0807 - α1007 = -0.060417230254934076 - 0.1965343807796979*α0807 - α1008 = 2.4846281621788395 - 2.102961615944379*α0807 + α1007 = 2.4846281621788395 - 2.102961615944379*α0807 α1101 = -1.0320124180911034 + 1.061943768952537*α0807 α1104 = 13.666683232895137 - 17.40407843561103*α0807 α1105 = 0.25990355211486116 - 3.797253476860588*α0807 @@ -698,7 +696,7 @@ end @muladd function perform_step!(integrator, cache::PFRK87Cache, repeat_step=false) @unpack t,dt,uprev,u,f,p = integrator - @unpack α0201, α0301, α0401, α0501, α0601, α0701, α0302, α0403, α0503, α0504, α0604, α0704, α0605, α0705, α0706, α1108, α1208, α1308, α1009, α1109, α1209, α1309, α1110, α1210, α1310, α1211, α1311, β1, β6, β7, β8, β9, β10, β11, β12, β13, β1tilde, β6tilde, β7tilde, β8tilde, β9tilde, β10tilde, β11tilde, β12tilde, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12, c13 = cache.tab + @unpack α0201, α0301, α0401, α0501, α0601, α0701, α0302, α0403, α0503, α0504, α0604, α0704, α0605, α0705, α0706, α0908, α1008, α1108, α1208, α1308, α1009, α1109, α1209, α1309, α1110, α1210, α1310, α1211, α1311, β1, β6, β7, β8, β9, β10, β11, β12, β13, β1tilde, β6tilde, β7tilde, β8tilde, β9tilde, β10tilde, β11tilde, β12tilde, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12, c13 = cache.tab @unpack k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,utilde,tmp,atmp,k = cache alg = unwrap_alg(integrator, true) @@ -715,17 +713,14 @@ end α0806 = 0.21311136 + 0.093456*α0807 α0901 = 0.07239997637512857 + 0.01913119863380767*α0807 α0904 = -0.688400520601143 - 0.3135390887207368*α0807 - α0905 = -0.688400520601143 - 0.3135390887207368*α0807 - α0906 = -0.17301267570583073 - 0.06840852844816077*α0807 - α0907 = 0.1440060555560846 + 0.031009360422930017*α0807 - α0908 = 0.9982362892760762 + 0.33180705811215994*α0807 + α0905 = -0.17301267570583073 - 0.06840852844816077*α0807 + α0906 = 0.1440060555560846 + 0.031009360422930017*α0807 + α0907 = 0.9982362892760762 + 0.33180705811215994*α0807 α1001 = 0.16261514523236525 - 0.12125171966747463*α0807 α1004 = -2.1255544052061124 + 1.9871809612169453*α0807 α1005 = -0.216403903283323 + 0.43356675517460624*α0807 - α1006 = -0.060417230254934076 - 0.1965343807796979*α0807 - α1007 = -0.060417230254934076 - 0.1965343807796979*α0807 - α1008 = 2.4846281621788395 - 2.102961615944379*α0807 + α1007 = 2.4846281621788395 - 2.102961615944379*α0807 α1101 = -1.0320124180911034 + 1.061943768952537*α0807 α1104 = 13.666683232895137 - 17.40407843561103*α0807 α1105 = 0.25990355211486116 - 3.797253476860588*α0807 diff --git a/src/tableaus/high_order_rk_tableaus.jl b/src/tableaus/high_order_rk_tableaus.jl index 8b6e16e2b2..8705ba6f94 100644 --- a/src/tableaus/high_order_rk_tableaus.jl +++ b/src/tableaus/high_order_rk_tableaus.jl @@ -1043,6 +1043,8 @@ struct PFRK87ConstantCache{T1,T2} <: OrdinaryDiffEqConstantCache α0706::T1 + α0908::T1 + α1008::T1 α1108::T1 α1208::T1 α1308::T1 @@ -1118,6 +1120,8 @@ function PFRK87ConstantCache(T1::Type, T2::Type) α0706 = convert(T1, 23124283//1800000000) + α0908 = convert(T1, 800635310//3783071287) + α1008 = convert(T1, 393006217//1396673457) α1108 = convert(T1, 15336726248//1032824649) α1208 = convert(T1, 5232866602//850066563) α1308 = convert(T1, -13158990841//6184727034) @@ -1166,5 +1170,5 @@ function PFRK87ConstantCache(T1::Type, T2::Type) c12 = convert(T2, 1//1) c13 = convert(T2, 1//1) - PFRK87ConstantCache(α0201, α0301, α0401, α0501, α0601, α0701, α0302, α0403, α0503, α0504, α0604, α0704, α0605, α0705, α0706, α1108, α1208, α1308, α1009, α1109, α1209, α1309, α1110, α1210, α1310, α1211, α1311, β1, β6, β7, β8, β9, β10, β11, β12, β13, β1tilde, β6tilde, β7tilde, β8tilde, β9tilde, β10tilde, β11tilde, β12tilde, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12, c13) + PFRK87ConstantCache(α0201, α0301, α0401, α0501, α0601, α0701, α0302, α0403, α0503, α0504, α0604, α0704, α0605, α0705, α0706, α0908, α1008, α1108, α1208, α1308, α1009, α1109, α1209, α1309, α1110, α1210, α1310, α1211, α1311, β1, β6, β7, β8, β9, β10, β11, β12, β13, β1tilde, β6tilde, β7tilde, β8tilde, β9tilde, β10tilde, β11tilde, β12tilde, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12, c13) end \ No newline at end of file diff --git a/test/algconvergence/ode_convergence_tests.jl b/test/algconvergence/ode_convergence_tests.jl index 2a18b50f73..855c44db41 100644 --- a/test/algconvergence/ode_convergence_tests.jl +++ b/test/algconvergence/ode_convergence_tests.jl @@ -11,6 +11,7 @@ dts1 = 1 .//2 .^(9:-1:5) dts2 = 1 .//2 .^(7:-1:3) dts3 = 1 .//2 .^(12:-1:7) dts4 = 1 .//2 .^(5:-1:3) +dts5 = 1 .//2 .^(3:-1:1) testTol = 0.2 @testset "Explicit Solver Convergence Tests ($(["out-of-place", "in-place"][i]))" for i in 1:2 @@ -40,6 +41,9 @@ testTol = 0.2 sim3 = test_convergence(dts4,prob,FRK65()) @test sim3.𝒪est[:l∞] ≈ 6 atol=0.6 + sim3 = test_convergence(dts5,prob,PFRK87()) + @test sim3.𝒪est[:l∞] ≈ 8.4 atol=0.2 + sim4 = test_convergence(dts,prob,BS3()) @test sim4.𝒪est[:l2] ≈ 3 atol=testTol sim5 = test_convergence(dts, prob, AB3())