From 18caa962df4d02127fd893e884c72bbc72aba557 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Mon, 10 Mar 2025 16:51:04 +0100 Subject: [PATCH 01/38] Added MPDeC for conservative PDS --- src/PositiveIntegrators.jl | 4 + src/mpdec.jl | 616 +++++++++++++++++++++++++++++++++++++ 2 files changed, 620 insertions(+) create mode 100644 src/mpdec.jl diff --git a/src/PositiveIntegrators.jl b/src/PositiveIntegrators.jl index 51306dec..3ac47a5d 100644 --- a/src/PositiveIntegrators.jl +++ b/src/PositiveIntegrators.jl @@ -48,6 +48,7 @@ export ConservativePDSFunction, ConservativePDSProblem export MPE, MPRK22, MPRK43I, MPRK43II export SSPMPRK22, SSPMPRK43 +export MPDeC export prob_pds_linmod, prob_pds_linmod_inplace, prob_pds_nonlinmod, prob_pds_robertson, prob_pds_brusselator, prob_pds_sir, @@ -69,6 +70,9 @@ include("mprk.jl") # modified Patankar-Runge-Kutta based on the SSP formulation of RK methods (SSPMPRK) include("sspmprk.jl") +# MPDeC methods +include("mpdec.jl") + # interpolation for dense output include("interpolation.jl") diff --git a/src/mpdec.jl b/src/mpdec.jl new file mode 100644 index 00000000..41792894 --- /dev/null +++ b/src/mpdec.jl @@ -0,0 +1,616 @@ +""" + MPDeC(K; [nodes = :gausslobatto, linsolve = ..., small_constant = ...]) + +EVERYTHING BELOW IS COPIED FROM MPRK22 + +A family of arbitrary order modified Patankar-Runge-Kutta algorithms for +production-destruction systems. Each member of this family is an adaptive, one-step method which is +Kth order accurate, unconditionally positivity-preserving, and linearly +implicit. + +This method supports adaptive time stepping, using TODO: ???, as lower order approximations to estimate the error. + +These schemes were introduced by Torlo and Öffner (2020) for conservative production-destruction systems. +For nonconservative production–destruction systems we use a straight forward extension +analogous to [`MPE`](@ref). + +TODO: We need a sentences about the nodes choices. Literature for Gauss-Lobatto? + +The MPDeC methods require the special structure of a +[`PDSProblem`](@ref) or a [`ConservativePDSProblem`](@ref). + +You can optionally choose the linear solver to be used by passing an +algorithm from [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl) +as keyword argument `linsolve`. +You can also choose the parameter `small_constant` which is added to all Patankar-weight denominators +to avoid divisions by zero. You can pass a value explicitly, otherwise `small_constant` is set to +`floatmin` of the floating point type used. + +## References + +- Davide Torlo, and Philipp Öffner. + "Arbitrary high-order, conservative and positivity preserving Patankar-type deferred correction schemes." + Applied Numerical Mathematics 153 (2020): 15-34. +""" +struct MPDeC{T, N, F, T2} <: OrdinaryDiffEqAlgorithm + K::T + M::T + nodes::N + linsolve::F + small_constant_function::T2 +end + +function MPDeC(K; nodes = :gausslobatto, linsolve = LUFactorization(), + small_constant = nothing) + if !(isinteger(K)) + throw(ArgumentError("MPDeC requires the parameter K to be an integer.")) + end + if !(typeof(K) <: Integer) + K = Int(K) + end + + if isnothing(small_constant) + small_constant_function = floatmin + elseif small_constant isa Number + small_constant_function = Returns(small_constant) + else # assume small_constant isa Function + small_constant_function = small_constant + end + + if nodes == :lagrange + M = K - 1 + else # :gausslobatto + M = ceil(Integer, K / 2) + end + + MPDeC{typeof(K), typeof(nodes), typeof(linsolve), typeof(small_constant_function)}(K, M, + nodes, + linsolve, + small_constant_function) +end + +alg_order(alg::MPDeC) = alg.K +isfsal(::MPDeC) = false + +function get_constant_parameters(alg::MPDeC) + if alg.nodes == :lagrange + if alg.M == 1 + theta = [0.0 0.5; 0.0 0.5] + elseif alg.M == 2 + theta = [0.0 0.20833333333333337 0.16666666666666652; + 0.0 0.33333333333333337 0.6666666666666667; + 0.0 -0.04166666666666667 0.16666666666666663] + elseif alg.M == 3 + theta = [0.0 0.125 0.11111111111111116 0.125; + 0.0 0.26388888888888895 0.4444444444444451 0.3750000000000009; + 0.0 -0.0694444444444445 0.11111111111111072 0.375; + 0.0 0.013888888888888888 5.551115123125783e-17 0.12499999999999978] + elseif alg.M == 4 + theta = [0.0 0.08715277777777783 0.0805555555555556 0.08437500000000009 0.07777777777777928; + 0.0 0.2243055555555556 0.34444444444444455 0.31874999999999964 0.35555555555555785; + 0.0 -0.09166666666666667 0.06666666666666643 0.22499999999999964 0.13333333333333286; + 0.0 0.036805555555555564 0.011111111111111016 0.1312500000000001 0.35555555555555607; + 0.0 -0.0065972222222222265 -0.0027777777777778234 -0.009374999999999911 0.0777777777777775] + elseif alg.M == 5 + theta = [0.0 0.06597222222222225 0.0622222222222224 0.06374999999999953 0.06222222222222573 0.06597222222222676; + 0.0 0.1981944444444445 0.28666666666666707 0.2737500000000015 0.28444444444444983 0.2604166666666785; + 0.0 -0.1108333333333335 0.031111111111109757 0.14249999999999297 0.10666666666664959 0.1736111111111054; + 0.0 0.06694444444444442 0.031111111111110645 0.1424999999999983 0.2844444444444356 0.17361111111109295; + 0.0 -0.024027777777777766 -0.013333333333332975 -0.026250000000000107 0.06222222222222484 0.2604166666666785; + 0.0 0.0037500000000000033 0.0022222222222222365 0.003750000000000364 0.0 0.06597222222222454] + elseif alg.M == 6 + theta = [0.0 0.05259865520282189 0.05022045855379187 0.05096726190476186 0.050440917107582806 0.05118772045855735 0.04880952380950987; + 0.0 0.179431216931217 0.2486772486772483 0.24107142857142755 0.24550264550265144 0.239748677248647 0.25714285714280294; + 0.0 -0.12803406084656088 0.0014550264550234893 0.08638392857142296 0.06772486772487873 0.08783895502639893 0.0321428571428477; + 0.0 0.1033509700176369 0.058553791887126866 0.16190476190476288 0.26525573192239804 0.22045855379185753 0.3238095238095582; + 0.0 -0.055696097883597806 -0.03558201058201127 -0.05424107142856682 0.030687830687839757 0.16017691798938927 0.03214285714290277; + 0.0 0.017394179894179934 0.011640211640212117 0.01607142857142918 0.008465608465618502 0.07771164021165333 0.2571428571429166; + 0.0 -0.00237819664902998 -0.0016313932980599258 -0.002157738095238018 -0.0014109347442683717 -0.003789131393296341 0.048809523809520694] + elseif alg.M == 7 + theta = [0.0 0.04346064814814822 0.0418367346938779 0.04225127551020513 0.04202569916855836 0.042251275510214015 0.041836734693902144 0.04346064814816941; + 0.0 0.1651655801209374 0.22161753590324995 0.21667729591836782 0.21889644746787695 0.21686626039310397 0.2204081632653896 0.20700231481467313; + 0.0 -0.14384566326530607 -0.02414965986394635 0.04390943877551834 0.032653061224550584 0.04118835034029189 0.02755102040850943 0.0765625000003638; + 0.0 0.1454235166288738 0.09251700680271657 0.18899872448978527 0.26969009826137835 0.24580144557796757 0.2775510204078273 0.17297453703622523; + 0.0 -0.10457648337112646 -0.07282690854119345 -0.09671556122447367 -0.016024187452693184 0.08045753023436752 0.02755102040857693 0.17297453703748644; + 0.0 0.04901147959183677 0.03537414965986363 0.04390943877552189 0.03265306122450795 0.1007121598638605 0.22040816326534696 0.07656249999990905; + 0.0 -0.013405848450491272 -0.009863945578231281 -0.011894132653058165 -0.009674981103547253 -0.014615221088391195 0.04183673469395899 0.2070023148147584; + 0.0 0.001623913454270594 0.0012093726379440242 0.0014349489795916492 0.0012093726379429626 0.0016239134542663791 -1.0658141036401503e-14 0.04346064814814099] + elseif alg.M == 8 + theta = [0.0 0.03685850005511468 0.035688932980599594 0.03594029017857134 0.0358289241622568 0.035914937444895045 0.035803571428594694 0.03605492862659787 0.034885361552085214; + 0.0 0.15387641920194012 0.2012610229276906 0.19782924107142996 0.1990828924162642 0.19819740685622378 0.1992857142857929 0.1969121334875581 0.2076895943564523; + 0.0 -0.15861283344356245 -0.046840828924162636 0.009592633928562577 0.0021516754850381403 0.0065018050044045594 0.0016071428572104196 0.011744309413188603 -0.03273368606460281; + 0.0 0.19274133322310427 0.13237213403880332 0.22303013392857451 0.2888183421517425 0.2741522679677928 0.2878571428567511 0.2618484760821502 0.3702292768929283; + 0.0 -0.17337411816578485 -0.12799823633157104 -0.15669642857142208 -0.08007054673731773 -0.003444664903213379 -0.0321428571432989 0.013233024690180173 -0.16014109348088823; + 0.0 0.10838080081569673 0.08237213403880084 0.09607700892858873 0.08141093474421268 0.14719914296711067 0.23785714285577342 0.1774879436707124 0.37022927689395146; + 0.0 -0.044477995480599525 -0.03434082892416224 -0.03923549107142321 -0.03488536155200972 -0.04232631999563363 0.01410714285702852 0.1258791473759011 -0.03273368606702931; + 0.0 0.010777460868606693 0.008403880070546516 0.009492187499996696 0.008606701940021111 0.009860353284821599 0.006428571428614305 0.053813175154459714 0.2076895943564523; + 0.0 -0.0011695670745149895 -0.0009182098765432678 -0.0010295758928574872 -0.0009435626102302086 -0.0010549286265453262 -0.0008035714285732354 -0.0019731385031036552 0.03488536155197153] + elseif alg.M == 9 + theta = [0.0 0.03188616071428564 0.031009210268469034 0.03117187499999874 0.031111111111105316 0.031149339236719698 0.031111111111073342 0.03117187499992724 0.031009210268393872 0.031886160714066136; + 0.0 0.14467159330295903 0.18532725847540765 0.1828236607142859 0.18359396433471176 0.1831509191895293 0.18357142857178133 0.18292556155574857 0.18461297276007826 0.17568080357159488; + 0.0 -0.1725594013325496 -0.06735057809132261 -0.0193750000000108 -0.02461297276119012 -0.022122403488083364 -0.02428571428754367 -0.021130829907633597 -0.02909660984732909 0.012053571401338559; + 0.0 0.24498946698020774 0.1776641191456143 0.26335317460326735 0.3186204193613946 0.30879507152621954 0.31587301587512684 0.3064180384141082 0.3290926905915512 0.21589285713457684; + 0.0 -0.2646060834313152 -0.20377621007249935 -0.23694196428579062 -0.16401332549457948 -0.10071817435664343 -0.11857142856752034 -0.098733067548892 -0.14234763861895772 0.06448660725436639; + 0.0 0.20683424578679332 0.1632196747011676 0.18305803571432477 0.1652047815015294 0.22849993263764645 0.3014285714274365 0.26826281722094336 0.329092690586549 0.06448660718228894; + 0.0 -0.11319983343131501 -0.09052518126592918 -0.09998015873019028 -0.09290221438379831 -0.10272756221900181 -0.04746031746481094 0.03822873799435911 -0.029096609840053134 0.21589285708614625; + 0.0 0.04115018126592196 0.03318440133255063 0.036339285714291236 0.034175974916735186 0.03666654418941562 0.03142857142916 0.07940414951872299 0.18461297276189725 0.0120535714286234; + 0.0 -0.008932169189692364 -0.007244757985499331 -0.007890625000002531 -0.007470115618263051 -0.00791316076327453 -0.007142857142902415 -0.009646454903545987 0.03100921026901915 0.17568080356954852; + 0.0 0.0008769504458161903 0.0007142857142857367 0.000775049603174427 0.0007368214775622661 0.0007750496031686538 0.0007142857142810044 0.0008769504457859512 5.684341886080802e-14 0.03188616071417982] + elseif alg.M == 10 + theta = [0.0 0.028018959644393687 0.027340374645930202 0.027451045048698663 0.0274153171930962 0.02743427683751065 0.027418831168816382 0.027437790813053198 0.02740206295788994 0.027512733360708808 0.02683414836155862; + 0.0 0.13699028395729754 0.17247367858478835 0.17057771915585285 0.17108139597029037 0.17083711202639051 0.17102597402617903 0.17080197227005556 0.17121393832712783 0.16996083603953593 0.17753594143869122; + 0.0 -0.18583978613015117 -0.08617167708834472 -0.04460141030846643 -0.048462401795699606 -0.04691594453948866 -0.048009740260582134 -0.046778097819711206 -0.04896713164639266 -0.04246829342450553 -0.08104357062261158; + 0.0 0.30192072009780363 0.22804745871412313 0.3094549512987534 0.35692031425357507 0.34993098144275514 0.35402597402539016 0.34980383697848083 0.3569305756250287 0.33648092529256246 0.45494628807136905; + 0.0 -0.3806483247655124 -0.3026606541606571 -0.3400126826297907 -0.2703953823952361 -0.21667333679079093 -0.2287597402537358 -0.2184080650258693 -0.23442039440851659 -0.19077242288039997 -0.43515512242447585; + 0.0 0.35715424082090674 0.29001218534554596 0.31687012987021035 0.29602437069246434 0.3568823152170353 0.41774025973563766 0.3968945005565274 0.4237524451623358 0.35661038970238224 0.7137646304354348; + 0.0 -0.24438269976550997 -0.2007347282347136 -0.21674705762993796 -0.20639538239549893 -0.2184817858628776 -0.16475974026434415 -0.0951424400191172 -0.13249446850386448 -0.05450679792556912 -0.43515512242447585; + 0.0 0.11846536295494622 0.09801571268237676 0.10514245129872465 0.10092031425360393 0.10501530683775329 0.09802597402490676 0.14549133697096295 0.22689882957820373 0.15302556815186108 0.45494628810411086; + 0.0 -0.038575277201579265 -0.03207643899310476 -0.03426547280844172 -0.03303383036714891 -0.03412762608719788 -0.032581168831256946 -0.03644216031966607 0.005128106447045866 0.10479621548802243 -0.08104357064439682; + 0.0 0.007575105385869255 0.006322003099781336 0.006733969155847452 0.006509967398880434 0.0066988293985943415 0.006454545454857907 0.006958222269304315 0.005062262842329801 0.04054565746628214 0.17753594143141527; + 0.0 -0.0006785849984634729 -0.0005679145956923939 -0.000603642451298736 -0.0005846828069051568 -0.0006001284755665637 -0.000581168831175205 -0.0006168966868074222 -0.0005062262839601317 -0.0011848112826555735 0.026834148362013366] + else + error("Theta Lagrange M > 9") + end + else # alg.nodes == :gausslobatto + if alg.M == 1 + theta = [0.0 0.5; 0.0 0.5] + elseif alg.M == 2 + theta = [0.0 0.20833333333333337 0.16666666666666652; + 0.0 0.33333333333333337 0.6666666666666667; + 0.0 -0.04166666666666667 0.16666666666666663] + elseif alg.M == 3 + theta = [0.0 0.11030056647916493 0.07303276685416865 0.08333333333333393; + 0.0 0.1896994335208352 0.45057403089581083 0.41666666666666785; + 0.0 -0.033907364229143935 0.22696723314583123 0.4166666666666661; + 0.0 0.010300566479164913 -0.026967233145831604 0.08333333333333326] + elseif alg.M == 4 + theta = [0.0 0.0677284321861569 0.04062499999999991 0.05370013924241501 0.05000000000000071; + 0.0 0.11974476934341162 0.30318418332304287 0.2615863979968083 0.27222222222222214; + 0.0 -0.021735721866558116 0.17777777777777748 0.3772912774221129 0.35555555555555296; + 0.0 0.010635824225415487 -0.0309619611008205 0.1524774528788102 0.27222222222222037; + 0.0 -0.0037001392424145354 0.009375000000000022 -0.017728432186157494 0.04999999999999982] + elseif alg.M == 5 + theta = [0.0 0.04567980513375505 0.025908385387879762 0.03746264288972734 0.03168990732937349 0.033333333333333215; + 0.0 0.08186781700897068 0.2138408086328255 0.177429781771262 0.19370925858950017 0.18923747814892522; + 0.0 -0.01487460578908985 0.13396073565086075 0.30143326325089315 0.2698015123994857 0.2774291885177327; + 0.0 0.007627676118250971 -0.024004074733154912 0.14346845286688126 0.2923037943068376 0.27742918851774334; + 0.0 -0.004471780440573713 0.011807696377659743 -0.02460333048390262 0.10736966113994595 0.18923747814891745; + 0.0 0.0016434260039545345 -0.004129309556393734 0.007424947945453564 -0.012346471800418257 0.03333333333333588] + elseif alg.M == 6 + theta = [0.0 0.03284626432829264 0.018002223201815104 0.027529761904763195 0.02170954269630787 0.02464779425215191 0.023809523809520172; + 0.0 0.059322894027551365 0.15770113064168867 0.1277882555598353 0.14409488824697547 0.13619592709181916 0.1384130236807124; + 0.0 -0.01076859445118926 0.10235481204686203 0.23748565272164512 0.20629511050420746 0.2193616205757678 0.21587269060495373; + 0.0 0.0055975917805697745 -0.018478259273458753 0.12190476190476146 0.2622877830829893 0.23821193202898883 0.24380952380948884; + 0.0 -0.0034889299708074674 0.009577580100741362 -0.02161296211671493 0.1135178785580706 0.22664128505616077 0.2158726906049253; + 0.0 0.002217096588914541 -0.005681864566224326 0.010624768120945982 -0.019288106960911655 0.07909012965322404 0.13841302368079056; + 0.0 -0.000838270442615105 0.0020999811132187685 -0.0037202380952379155 0.005807300607708399 -0.00903674051876635 0.02380952380952195] + elseif alg.M == 7 + theta = [0.0 0.024737514438875514 0.013258719822130019 0.021034356725210923 0.01577972028613317 0.019036721343624663 0.017385669593011244 0.01785714285716722; + 0.0 0.044892662602755 0.12064973282409763 0.09628017831208524 0.11096097521429726 0.10224716355590147 0.10657833455800869 0.10535211357202456; + 0.0 -0.008140767742389747 0.07999763578665159 0.1890416498443277 0.16114433643933257 0.17539401516496778 0.1687159595472414 0.17056134624175545; + 0.0 0.004254082548545937 -0.014480830962625313 0.10124595564769123 0.22436671754081416 0.19859744812260516 0.2089336024046844 0.20622939732966472; + 0.0 -0.002704205075015816 0.0076319492068338685 -0.018137320211455643 0.10498344168165086 0.22071022829197284 0.20197531478076014 0.20622939732958345; + 0.0 0.0018453866944819865 -0.004832668923134664 0.009417009802429988 -0.01848030360257269 0.09056371045502942 0.1787021139839453 0.17056134624169772; + 0.0 -0.0012262209861064882 0.00310495001592085 -0.005608861642537821 0.00907193525967287 -0.015297619252294226 0.060459450968835426 0.10535211357171193; + 0.0 0.0004714732640502682 -0.001179578486445107 0.002077422571009513 -0.003177213868071016 0.0045984230350075705 -0.006880371581776679 0.017857142857105046] + elseif alg.M == 8 + theta = [0.0 0.019293838201043245 0.01018408040822739 0.01656936984357027 0.011990017361096061 0.01514127656190567 0.013175811859724718 0.01417408466352299 0.013888888888914153; + 0.0 0.03512552097762184 0.09508650844936217 0.07509351709620726 0.08786983372308654 0.07945805621039881 0.08459518591985216 0.0820139081133675 0.08274768078126726; + 0.0 -0.006364102418704803 0.0639319956284379 0.15288206102488378 0.12868222230659399 0.14236568211174472 0.13451403216049584 0.13834341694882824 0.13726935625072656; + 0.0 0.0033337771969983894 -0.011585731333842608 0.0840854786889147 0.18975808031592045 0.16521365191611892 0.17719508127352412 0.1717199563247931 0.17321425548699665; + 0.0 -0.002136847017608247 0.006149936900551798 -0.015130673172334241 0.09287981859410088 0.20089031036043536 0.17960970028764223 0.18789648420624872 0.1857596371874024; + 0.0 0.0014942991627282308 -0.003980825787182701 0.008000603570397224 -0.0165438248294123 0.08912877679762232 0.18479998682039422 0.16988047828840536 0.17321425548630032; + 0.0 -0.001074060699381661 0.0027553240894192185 -0.0050963258617406915 0.008587133943507297 -0.015612704774753183 0.0733373606215082 0.14363345866922828 0.13726935624890757; + 0.0 0.0007337726665392911 -0.0018475051393835457 0.0032896245700402282 -0.005122152942657721 0.007654163684208015 -0.012338827668713748 0.04762215980304063 0.08274768078103989; + 0.0 -0.00028519577496630263 0.0007130770290413452 -0.0012523876730271555 0.0018988715277794554 -0.002680480954676767 0.0037048084806841075 -0.005404949312023177 0.013888888889184159] + elseif alg.M == 9 + theta = [0.0 0.015465148886122208 0.008074564828863144 0.013376018525479871 0.009424872116881033 0.012319010010855891 0.010311035538052238 0.01156730916181914 0.010928591594165482 0.011111111110039928; + 0.0 0.02821797600073163 0.07677451467523239 0.060184542475466785 0.07119971523047752 0.06348355183104104 0.06872200531256567 0.06548282490166457 0.06711913714752882 0.0666529954141879; + 0.0 -0.0051090759506765785 0.05211803626519543 0.12565307727248143 0.10482643149413295 0.11734393918791586 0.10937232662209873 0.11414511256730009 0.11177491077136281 0.1124446710354885; + 0.0 0.0026797026134601862 -0.009449705457448665 0.0703555323098648 0.1607094677113423 0.1383496630545551 0.15043195838711654 0.14368298780891564 0.14692275221750606 0.1460213418404237; + 0.0 -0.0017246233123523063 0.005034230848115648 -0.012689136533024836 0.08096672386450621 0.17827052516622643 0.15700518309495237 0.1670603392740304 0.16255069057842775 0.1637698805952823; + 0.0 0.0012191900098376621 -0.003290458683960145 0.006764697496483674 -0.014500644574205523 0.08280315672743654 0.17645901712467094 0.15873564974754117 0.16549450390448328 0.1637698805971013; + 0.0 -0.0009014103789669575 0.002338354035534229 -0.004410616548205404 0.007671678785357017 -0.014688125871543889 0.07566580953016455 0.15547104729739658 0.1433416392246727 0.1460213418395142; + 0.0 0.0006697602581937584 -0.00170044153617834 0.0030723444103197828 -0.0048992681562261 0.007618239537521276 -0.013208406240892145 0.06032663476707967 0.11755374698441301 0.1124446710346092; + 0.0 -0.000466141727447515 0.0011701705230945658 -0.002069009887745693 0.0031694435943515065 -0.004546719804944366 0.0064684529496542575 -0.010121519251002686 0.03843501942257066 0.0666529954255306; + 0.0 0.0001825195178684848 -0.0004561980512006905 0.0008000755736511378 -0.0012078988997608064 0.00168623899423892 -0.002264907414485151 0.00303654628237382 -0.004354037773737218 0.011111111112313665] + else + error("Theta Gauss-Lobatto M > 9") + end + end + return theta +end + +struct MPDeCConstantCache{KType, T, T2} <: OrdinaryDiffEqConstantCache + K::KType + M::KType + theta::T2 + small_constant::T +end + +# Out-of-place +function alg_cache(alg::MPDeC, u, rate_prototype, ::Type{uEltypeNoUnits}, + ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, + uprev, uprev2, f, t, dt, reltol, p, calck, + ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} + if !(f isa PDSFunction || f isa ConservativePDSFunction) + throw(ArgumentError("MPDEC can only be applied to production-destruction systems")) + end + + theta = get_constant_parameters(alg) + MPDeCConstantCache(alg.K, alg.M, theta, alg.small_constant_function(uEltypeNoUnits)) +end + +function initialize!(integrator, cache::MPDeCConstantCache) +end + +# out-of-place +function build_mpdec_matrix(m, prod, C, p, t, dt, theta, small_constant) + N, M = size(C) + M = M - 1 + + Mmat = zeros(eltype(C), N, N) + # P stores the production matrix in in-place computations. + # In out-of-place computations it is not needed + P = nothing + build_mpdec_matrix!(Mmat, m, prod, P, C, p, t, dt, theta, small_constant) + + if C isa StaticArray + return SMatrix(Mmat) + else + return Mmat + end +end + +# in-place for dense arrays +function build_mpdec_matrix!(Mmat, m, prod, P, C, p, t, dt, theta, small_constant) + N, M = size(C) + M = M - 1 + + fill!(Mmat, zero(eltype(Mmat))) + oneMmat = one(eltype(Mmat)) + for i in 1:N + Mmat[i, i] = oneMmat + end + + σ = add_small_constant(C[:, m], small_constant) + + for r in 1:(M + 1) + th = theta[r, m] + dt_th = dt * th + #TODO: This should be checked earlier and only once + if isinplace(prod, 4) + prod(P, C[:, r], p, t) + else + P = prod(C[:, r], p, t) + end + for i in 1:N + for j in 1:N + if th >= 0 + Mmat[i, j] -= dt_th * P[i, j] / σ[j] + Mmat[i, i] += dt_th * P[j, i] / σ[i] + else + Mmat[i, j] += dt_th * P[j, i] / σ[j] + Mmat[i, i] -= dt_th * P[i, j] / σ[i] + end + end + end + end +end + +@muladd function perform_step!(integrator, cache::MPDeCConstantCache, repeat_step = false) + @unpack alg, t, dt, uprev, f, p = integrator + @unpack K, M, theta, small_constant = cache + + C2 = zeros(length(uprev), M + 1) + for i in 1:(M + 1) + C2[:, i] = uprev + end + + if uprev isa StaticArray + #TODO make C static + C = copy(C2) + else + C = copy(C2) + end + + for _ in 1:K + for m in 2:(M + 1) + Mmat = build_mpdec_matrix(m, f.p, C, p, t, dt, theta, small_constant) + + # solve linear system + linprob = LinearProblem(Mmat, uprev) + sol = solve(linprob, alg.linsolve) + C2[:, m] = sol.u + integrator.stats.nsolve += 1 + end + C = copy(C2) + end + u2 = C[:, M + 1] + + #TODO: Remove this + @assert u2 ≈ MPDeC_check(K, M, uprev, theta, f, p, t, dt) + + integrator.u = u1 +end + +struct MPDeCCache{uType, PType, CType, tabType, F} <: MPRKMutableCache + tmp::uType + P::PType + P2::PType + D::uType + D2::uType + σ::uType + C::CType + C2::CType + tab::tabType + linsolve_rhs::uType + linsolve::F +end + +struct MPDeCConservativeCache{uType, PType, CType, tabType, F} <: MPRKMutableCache + P::PType + P2::PType + σ::uType + C::CType + C2::CType + tab::tabType + linsolve_rhs::uType + linsolve::F +end + +get_tmp_cache(integrator, ::MPDeC, cache::OrdinaryDiffEqMutableCache) = (cache.σ,) + +# In-place +function alg_cache(alg::MPDeC, u, rate_prototype, ::Type{uEltypeNoUnits}, + ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, + uprev, uprev2, f, t, dt, reltol, p, calck, + ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} + theta = get_constant_parameters(alg) + tab = MPDeCConstantCache(alg.K, alg.M, theta, + alg.small_constant_function(uEltypeNoUnits)) + + tmp = zero(u) + P = p_prototype(u, f) # stores evaluation of the production matrix + P2 = p_prototype(u, f) # stores the linear system matrix + σ = zero(u) + C = zeros(eltype(u), length(u), alg.M + 1) + C2 = zeros(eltype(u), length(u), alg.M + 1) + linsolve_rhs = zero(u) + + if f isa ConservativePDSFunction + # The right hand side of the linear system is always uprev. But using + # linsolve_rhs instead of uprev for the rhs we allow `alias_b=true`. uprev must + # not be altered, since it is needed to compute the adaptive time step + # size. + linprob = LinearProblem(P2, _vec(linsolve_rhs)) + linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true, + assumptions = LinearSolve.OperatorAssumptions(true)) + + MPDeCConservativeCache(P, P2, σ, C, C2, + tab, #MPDeCConstantCache + linsolve_rhs, + linsolve) + elseif f isa PDSFunction + linprob = LinearProblem(P2, _vec(linsolve_rhs)) + linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true, + assumptions = LinearSolve.OperatorAssumptions(true)) + + MPDeCCache(tmp, P, P2, + similar(u), # D + similar(u), # D2 + σ, + C, C2, + tab, #MPDeCConstantCache + linsolve_rhs, + linsolve) + else + throw(ArgumentError("MPDeC can only be applied to production-destruction systems")) + end +end + +function initialize!(integrator, cache::Union{MPDeCCache, MPDeCConservativeCache}) +end + +#= COPY OF MPRK22 +@muladd function perform_step!(integrator, cache::MPDeCCache, repeat_step = false) + @unpack t, dt, uprev, u, f, p = integrator + @unpack tmp, P, P2, D, D2, σ, linsolve = cache + @unpack a21, b1, b2, small_constant = cache.tab + + # We use P2 to store the last evaluation of the PDS + # as well as to store the system matrix of the linear system + + f.p(P, uprev, p, t) # evaluate production terms + f.d(D, uprev, p, t) # evaluate nonconservative destruction terms + integrator.stats.nf += 1 + if issparse(P) + # We need to keep the structural nonzeros of the production terms. + # However, this is not guaranteed by broadcasting, see + # https://github.com/JuliaSparse/SparseArrays.jl/issues/190 + nz_P = nonzeros(P) + nz_P2 = nonzeros(P2) + @.. broadcast=false nz_P2=a21 * nz_P + else + @.. broadcast=false P2=a21 * P + end + @.. broadcast=false D2=a21 * D + + # avoid division by zero due to zero Patankar weights + @.. broadcast=false σ=uprev + small_constant + + # tmp holds the right hand side of the linear system + tmp .= uprev + @inbounds for i in eachindex(tmp) + tmp[i] += dt * P2[i, i] + end + + build_mprk_matrix!(P2, P2, σ, dt, D2) + + # Same as linres = P2 \ tmp + linsolve.A = P2 + linres = solve!(linsolve) + + u .= linres + integrator.stats.nsolve += 1 + + if isone(a21) + σ .= u + else + @.. broadcast=false σ=σ^(1 - 1 / a21) * u^(1 / a21) + end + @.. broadcast=false σ=σ + small_constant + + f.p(P2, u, p, t + a21 * dt) # evaluate production terms + f.d(D2, u, p, t + a21 * dt) # evaluate nonconservative destruction terms + integrator.stats.nf += 1 + + if issparse(P) + # We need to keep the structural nonzeros of the production terms. + # However, this is not guaranteed by broadcasting, see + # https://github.com/JuliaSparse/SparseArrays.jl/issues/190 + nz_P = nonzeros(P) + nz_P2 = nonzeros(P2) + @.. broadcast=false nz_P2=b1 * nz_P + b2 * nz_P2 + else + @.. broadcast=false P2=b1 * P + b2 * P2 + end + @.. broadcast=false D2=b1 * D + b2 * D2 + + # tmp holds the right hand side of the linear system + tmp .= uprev + @inbounds for i in eachindex(tmp) + tmp[i] += dt * P2[i, i] + end + + build_mprk_matrix!(P2, P2, σ, dt, D2) + + # Same as linres = P2 \ tmp + linsolve.A = P2 + linres = solve!(linsolve) + + u .= linres + integrator.stats.nsolve += 1 + + # Now σ stores the error estimate + # If a21 = 1, then σ is the MPE approximation, i.e. suited for stiff problems. + # Otherwise, this is not clear. + @.. broadcast=false σ=u - σ + + # Now tmp stores error residuals + calculate_residuals!(tmp, σ, uprev, u, integrator.opts.abstol, + integrator.opts.reltol, integrator.opts.internalnorm, t, + False()) + integrator.EEst = integrator.opts.internalnorm(tmp, t) +end +=# + +@muladd function perform_step!(integrator, cache::MPDeCConservativeCache, + repeat_step = false) + @unpack t, dt, uprev, u, f, p = integrator + @unpack P, P2, σ, C, C2, linsolve_rhs, linsolve = cache + @unpack K, M, theta, small_constant = cache.tab + + # Set right hand side of linear system + linsolve_rhs .= uprev + + # Initializ C matrices + for i in 1:(M + 1) + C2[:, i] = uprev + C[:, i] = uprev + end + + for _ in 1:K + for m in 2:(M + 1) + build_mpdec_matrix!(P2, m, f.p, P, C, p, t, dt, theta, small_constant) + + # Same as linres = P2 \ linsolve_rhs + linsolve.A = P2 + linres = solve!(linsolve) + C2[:, m] .= linres + integrator.stats.nsolve += 1 + end + C .= C2 + end + + u .= C[:, M + 1] +end + +######################################################################################################## +######################################################################################################## +#### The functions below are provided with the original paper and used for validation. #### +#### They will be removed in the future. #### +######################################################################################################## +######################################################################################################## +function patanker_type_dec(prod_p, dest_p, delta_t, m_substep::Int, M_sub, theta, u_p, dim) + #This function builds and solves the system u^{(k+1)}=Mu^{(k)} for a subtimestep m + + mass = Matrix{Float64}(I, dim, dim) + #println(mass) + for i in 1:dim + for r in 1:(M_sub + 1) + if theta[r, m_substep] > 0 + for j in 1:dim + mass[i, j] = mass[i, j] - + delta_t * theta[r, m_substep] * + (prod_p[i, j, r] / u_p[j, m_substep]) + mass[i, i] = mass[i, i] + + delta_t * theta[r, m_substep] * + (dest_p[i, j, r] / u_p[i, m_substep]) + end + elseif theta[r, m_substep] < 0 + for j in 1:dim + mass[i, i] = mass[i, i] - + delta_t * theta[r, m_substep] * + (prod_p[i, j, r] / u_p[i, m_substep]) + mass[i, j] = mass[i, j] + + delta_t * theta[r, m_substep] * + (dest_p[i, j, r] / u_p[j, m_substep]) + end + end + end + end + return mass \ u_p[:, 1] +end + +function MPDeC_check(K, M, uprev, theta, f, p, t, dt) + M_sub = M + K_corr = K + dim = length(uprev) + + #Variables of the DeC procedure: one for each subtimesteps and one for previous correction up, one for current correction ua + u_p = zeros(dim, M_sub + 1) + u_a = zeros(dim, M_sub + 1) + + #Matrices of production and destructions applied to up at different subtimesteps + prod_p = zeros(dim, dim, M_sub + 1) + dest_p = zeros(dim, dim, M_sub + 1) + + #coefficients for time integration + Theta = theta + + # Subtimesteps loop to update variables at iteration 0 + for m in 1:(M_sub + 1) + u_a[:, m] = uprev + u_p[:, m] = uprev + end + + #Loop for iterations K + for k in 2:(K_corr + 1) + # Updating previous iteration variabls + u_p = copy(u_a) + #Loop on subtimesteps to compute the production destruction terms + for r in 1:(M_sub + 1) + #prod_p[:,:,r], dest_p[:,:,r]=prod_dest(u_p[:,r]) + prod_p[:, :, r] = f.p(u_p[:, r], p, t) + dest_p[:, :, r] = prod_p[:, :, r]' + end + # Loop on subtimesteps to compute the new variables + for m in 2:(M_sub + 1) + u_a[:, m] = patanker_type_dec(prod_p, dest_p, dt, m, M_sub, Theta, u_p, + dim) + end + end + u1 = u_a[:, M_sub + 1] + return u1 +end From 0b0de40de46a2555a715d03b7cb6d3d2b35efd53 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Mon, 10 Mar 2025 17:38:10 +0100 Subject: [PATCH 02/38] fix --- src/mpdec.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mpdec.jl b/src/mpdec.jl index 41792894..01d57241 100644 --- a/src/mpdec.jl +++ b/src/mpdec.jl @@ -331,7 +331,7 @@ end #TODO: Remove this @assert u2 ≈ MPDeC_check(K, M, uprev, theta, f, p, t, dt) - integrator.u = u1 + integrator.u = u2 end struct MPDeCCache{uType, PType, CType, tabType, F} <: MPRKMutableCache From 6a14d62a622f7ab8878dfda89cc7cbdbd8233821 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Mon, 10 Mar 2025 18:33:09 +0100 Subject: [PATCH 03/38] first tests --- src/mpdec.jl | 2 +- test/runtests.jl | 30 +++++++++++++++++++++++++++--- 2 files changed, 28 insertions(+), 4 deletions(-) diff --git a/src/mpdec.jl b/src/mpdec.jl index 01d57241..a7b08cf0 100644 --- a/src/mpdec.jl +++ b/src/mpdec.jl @@ -329,7 +329,7 @@ end u2 = C[:, M + 1] #TODO: Remove this - @assert u2 ≈ MPDeC_check(K, M, uprev, theta, f, p, t, dt) + #@assert u2 ≈ MPDeC_check(K, M, uprev, theta, f, p, t, dt) integrator.u = u2 end diff --git a/test/runtests.jl b/test/runtests.jl index 02706689..e3a8fcb5 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1092,6 +1092,15 @@ end @test_throws "SSPMPRK43 can only be applied to production-destruction systems" solve(prob_ip, SSPMPRK43(), dt = 0.1) + @test_throws "MPDeC can only be applied to production-destruction systems" solve(prob_ip, + MPDeC(2), + dt = 0.1) + @test_throws "MPDeC can only be applied to production-destruction systems" solve(prob_oop, + MPDeC(2), + dt = 0.1) + @test_throws "MPDeC requires the parameter K to be an integer." solve(prob_pds_linmod, + MPDeC(2.1), + dt = 0.1) end # Here we check that algorithms which accept input parameters return constants @@ -1121,7 +1130,8 @@ end algs = (MPRK22(0.5f0), MPRK22(1.0f0), MPRK22(2.0f0), MPRK43I(1.0f0, 0.5f0), MPRK43I(0.5f0, 0.75f0), MPRK43II(0.5f0), MPRK43II(2.0f0 / 3.0f0), - SSPMPRK22(0.5f0, 1.0f0), SSPMPRK43()) + SSPMPRK22(0.5f0, 1.0f0), SSPMPRK43(), MPDeC(2), + MPDeC(2, nodes = :lagrange)) for alg in algs sol = solve(prob, alg; dt = 0.1f0, save_everystep = false) @test sol.t isa Vector{Float32} @@ -1185,6 +1195,7 @@ end # Here we check that MPRK22(α) = SSPMPRK22(0,α) @testset "MPRK22(α) = SSPMPRK22(0, α)" begin + #TODO: Add MPDeC(2) and MPDeC(2, nodes=:lagrange) REQUIRES ADAPTIVITY for α in (0.5, 2.0 / 3.0, 1.0, 2.0) # conservative PDS sol1 = solve(prob_pds_linmod, MPRK22(α)) @@ -1204,6 +1215,7 @@ end # Here we check that different linear solvers can be used @testset "Different linear solvers" begin + #TODO: Add MPDeC - REQUIRES PDSProblem # problem data u0 = [0.9, 0.1] tspan = (0.0, 2.0) @@ -1264,6 +1276,7 @@ end # Here we check that in-place and out-of-place implementations # deliver the same results @testset "Different matrix types (conservative)" begin + #TODO: Add MPDeC - REQUIRES spare and tridiagonal matrices prod_1! = (P, u, p, t) -> begin fill!(P, zero(eltype(P))) for i in 1:(length(u) - 1) @@ -1348,6 +1361,8 @@ end end @testset "Different matrix types (conservative, adaptive)" begin + #TODO: SSPMPRK43 is not adaptive and needs to be removed + #TODO: Add MPDeC prod_1! = (P, u, p, t) -> begin fill!(P, zero(eltype(P))) for i in 1:(length(u) - 1) @@ -1434,6 +1449,7 @@ end # Here we check that in-place and out-of-place implementations # deliver the same results @testset "Different matrix types (nonconservative)" begin + #TODO: Add MPDeC prod_1! = (P, u, p, t) -> begin fill!(P, zero(eltype(P))) for i in 1:(length(u) - 1) @@ -1559,6 +1575,8 @@ end end @testset "Different matrix types (nonconservative, adaptive)" begin + #TODO: SSPMPRK43 is not adaptive and needs to be removed + #TODO: Add MPDeC - Requires PDSProblem and adaptivity prod_1! = (P, u, p, t) -> begin fill!(P, zero(eltype(P))) for i in 1:(length(u) - 1) @@ -1690,6 +1708,7 @@ end # defines the types of the Ps inside the algorithm caches. # We test sparse, tridiagonal, and dense matrices. @testset "Prototype type check" begin + #TODO: Add MPDeC #prod and dest functions prod_inner! = (P, u, p, t) -> begin fill!(P, zero(eltype(P))) @@ -1770,7 +1789,8 @@ end # Here we check the convergence order of pth-order schemes for which # also an interpolation of order p is available @testset "Convergence tests (conservative)" begin - algs = (MPE(), MPRK22(0.5), MPRK22(1.0), MPRK22(2.0), SSPMPRK22(0.5, 1.0)) + algs = (MPE(), MPRK22(0.5), MPRK22(1.0), MPRK22(2.0), SSPMPRK22(0.5, 1.0), + MPDeC(2), MPDeC(2, nodes = :lagrange)) dts = 0.5 .^ (4:15) problems = (prob_pds_linmod, prob_pds_linmod_array, prob_pds_linmod_mvector, prob_pds_linmod_inplace) @@ -1806,6 +1826,7 @@ end # Here we check the convergence order of pth-order schemes for which # also an interpolation of order p is available @testset "Convergence tests (nonconservative)" begin + #TODO: Add MPDeC - PDSProblem algs = (MPE(), MPRK22(0.5), MPRK22(1.0), MPRK22(2.0), SSPMPRK22(0.5, 1.0)) dts = 0.5 .^ (4:15) problems = (prob_pds_linmod_nonconservative, @@ -1837,6 +1858,7 @@ end end end + #TODO: Check order of MPDeC(K) for K ≥ 4 # Here we check the convergence order of pth-order schemes for which # no interpolation of order p is available @testset "Convergence tests (conservative)" begin @@ -1844,7 +1866,8 @@ end problems = (prob_pds_linmod, prob_pds_linmod_array, prob_pds_linmod_mvector, prob_pds_linmod_inplace) algs = (MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), - MPRK43II(0.5), MPRK43II(2.0 / 3.0), SSPMPRK43()) + MPRK43II(0.5), MPRK43II(2.0 / 3.0), SSPMPRK43(), MPDeC(3), + MPDeC(3, nodes = :lagrange)) for alg in algs, prob in problems orders = experimental_orders_of_convergence(prob, alg, dts) @test check_order(orders, PositiveIntegrators.alg_order(alg), atol = 0.2) @@ -1854,6 +1877,7 @@ end # Here we check the convergence order of pth-order schemes for which # no interpolation of order p is available @testset "Convergence tests (nonconservative)" begin + #TODO: Add MPDeC(3) - PDSProblem dts = 0.5 .^ (4:12) problems = (prob_pds_linmod_nonconservative, prob_pds_linmod_nonconservative_inplace) From 80f699af2f0539cd8f00817051750569e0710d9f Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Mon, 10 Mar 2025 21:03:04 +0100 Subject: [PATCH 04/38] fix error test --- src/mpdec.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mpdec.jl b/src/mpdec.jl index a7b08cf0..65d0beb9 100644 --- a/src/mpdec.jl +++ b/src/mpdec.jl @@ -234,7 +234,7 @@ function alg_cache(alg::MPDeC, u, rate_prototype, ::Type{uEltypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck, ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} if !(f isa PDSFunction || f isa ConservativePDSFunction) - throw(ArgumentError("MPDEC can only be applied to production-destruction systems")) + throw(ArgumentError("MPDeC can only be applied to production-destruction systems")) end theta = get_constant_parameters(alg) From 723e02dc4b29cf023231aae687abe479cded1070 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Mon, 10 Mar 2025 21:10:01 +0100 Subject: [PATCH 05/38] added TODOs in tests --- test/runtests.jl | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/test/runtests.jl b/test/runtests.jl index e3a8fcb5..aebe292d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1890,6 +1890,7 @@ end end @testset "Interpolation tests (conservative)" begin + #TODO: Add MPDeC algs = (MPE(), MPRK22(0.5), MPRK22(1.0), MPRK22(2.0), MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), MPRK43II(0.5), MPRK43II(2.0 / 3.0), SSPMPRK22(0.5, 1.0), SSPMPRK43()) @@ -1908,6 +1909,7 @@ end end @testset "Check convergence order (nonautonomous conservative PDS)" begin + #TODO: Add MPDeC prod! = (P, u, p, t) -> begin fill!(P, zero(eltype(P))) P[1, 2] = sin(t)^2 * u[2] @@ -1936,6 +1938,7 @@ end end @testset "Check convergence order (nonautonomous nonconservative PDS)" begin + #TODO: Add MPDeC prod! = (P, u, p, t) -> begin fill!(P, zero(eltype(P))) P[1, 2] = sin(t)^2 * u[2] @@ -1994,6 +1997,7 @@ end # Check that the schemes accept zero initial values @testset "Zero initial values" begin + #TODO: Add MPDeC - Requires PDSProblem # Do a single step and check that no NaNs occur u0 = [1.0, 0.0] dt = 1.0 @@ -2042,6 +2046,7 @@ end # Check that approximations, and thus the Patankar weights, # remain positive to avoid division by zero. @testset "Positvity check" begin + #TODO: Add MPDeC - Requires PDSProblem # For this problem u[1] decreases montonically to 0 very fast. # We perform 10^5 steps and check that u[end] does not contain any NaNs u0 = [0.9, 0.1] @@ -2092,6 +2097,7 @@ end # Here we check that the implemented schemes can solve the predefined PDS # (at least for specific parameters) @testset "PDS problem library (adaptive schemes)" begin + #TODO: Add MPDeC algs = (MPRK22(0.5), MPRK22(1.0), MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), MPRK43II(2.0 / 3.0), MPRK43II(0.5), SSPMPRK22(0.5, 1.0)) probs = (prob_pds_linmod, prob_pds_linmod_inplace, prob_pds_nonlinmod, @@ -2135,6 +2141,7 @@ end #Here we check the different possibilities to define small_constant. @testset "Different possibilities to set small_constant" begin + #TODO: Add MPDeC - Requires PDSProblem # For this problem u[1] decreases montonically to 0 very fast. u0 = [0.9, 0.1] tspan = (0.0, 100.0) @@ -2188,6 +2195,7 @@ end #Here we check if the RK methods on which the MPRK schemes are based integrate # u'(t) = q * t^(q-1) exactly for q from 1 to the order of the method. @testset "Exact solutions (RK)" begin + #TODO: Add MPDeC - Requires PDSProblem algs = (MPE(), MPRK22(0.5), MPRK22(1.0), MPRK22(2.0), MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), MPRK43II(0.5), MPRK43II(2.0 / 3.0), SSPMPRK22(0.5, 1.0), SSPMPRK43()) From 6e42eabe34574ec1232d58ccb5b9c46ca625fc3e Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Mon, 10 Mar 2025 21:44:27 +0100 Subject: [PATCH 06/38] added mpdec api reference --- docs/src/api_reference.md | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/src/api_reference.md b/docs/src/api_reference.md index ac1a6f9f..c9f26925 100644 --- a/docs/src/api_reference.md +++ b/docs/src/api_reference.md @@ -35,6 +35,7 @@ SSPMPRK22 MPRK43I MPRK43II SSPMPRK43 +MPDeC ``` ## Auxiliary functions From 1163864ad02f5d2192013626346293f27a79ae86 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Wed, 12 Mar 2025 12:22:36 +0100 Subject: [PATCH 07/38] Added MPDeC for non-autonomous PDS --- src/mpdec.jl | 92 +++++++++++++++++++++++++++++++++++++++--------- test/runtests.jl | 24 ++++++++----- 2 files changed, 90 insertions(+), 26 deletions(-) diff --git a/src/mpdec.jl b/src/mpdec.jl index 65d0beb9..51c1caa7 100644 --- a/src/mpdec.jl +++ b/src/mpdec.jl @@ -74,6 +74,7 @@ isfsal(::MPDeC) = false function get_constant_parameters(alg::MPDeC) if alg.nodes == :lagrange + nodes = collect(0.0:(1 / alg.M):1.0) if alg.M == 1 theta = [0.0 0.5; 0.0 0.5] elseif alg.M == 2 @@ -153,23 +154,35 @@ function get_constant_parameters(alg::MPDeC) end else # alg.nodes == :gausslobatto if alg.M == 1 + nodes = [0.0, 1.0] theta = [0.0 0.5; 0.0 0.5] elseif alg.M == 2 + nodes = [0.0, 0.5, 1.0] theta = [0.0 0.20833333333333337 0.16666666666666652; 0.0 0.33333333333333337 0.6666666666666667; 0.0 -0.04166666666666667 0.16666666666666663] elseif alg.M == 3 + nodes = [0.0, 0.27639320225002106, 0.7236067977499789, 1.0] theta = [0.0 0.11030056647916493 0.07303276685416865 0.08333333333333393; 0.0 0.1896994335208352 0.45057403089581083 0.41666666666666785; 0.0 -0.033907364229143935 0.22696723314583123 0.4166666666666661; 0.0 0.010300566479164913 -0.026967233145831604 0.08333333333333326] elseif alg.M == 4 + nodes = [0.0, 0.1726731646460114, 0.5, 0.8273268353539887, 1.0] theta = [0.0 0.0677284321861569 0.04062499999999991 0.05370013924241501 0.05000000000000071; 0.0 0.11974476934341162 0.30318418332304287 0.2615863979968083 0.27222222222222214; 0.0 -0.021735721866558116 0.17777777777777748 0.3772912774221129 0.35555555555555296; 0.0 0.010635824225415487 -0.0309619611008205 0.1524774528788102 0.27222222222222037; 0.0 -0.0037001392424145354 0.009375000000000022 -0.017728432186157494 0.04999999999999982] elseif alg.M == 5 + nodes = [ + 0.0, + 0.11747233803526769, + 0.3573842417596774, + 0.6426157582403226, + 0.8825276619647323, + 1.0 + ] theta = [0.0 0.04567980513375505 0.025908385387879762 0.03746264288972734 0.03168990732937349 0.033333333333333215; 0.0 0.08186781700897068 0.2138408086328255 0.177429781771262 0.19370925858950017 0.18923747814892522; 0.0 -0.01487460578908985 0.13396073565086075 0.30143326325089315 0.2698015123994857 0.2774291885177327; @@ -177,6 +190,15 @@ function get_constant_parameters(alg::MPDeC) 0.0 -0.004471780440573713 0.011807696377659743 -0.02460333048390262 0.10736966113994595 0.18923747814891745; 0.0 0.0016434260039545345 -0.004129309556393734 0.007424947945453564 -0.012346471800418257 0.03333333333333588] elseif alg.M == 6 + nodes = [ + 0.0, + 0.08488805186071652, + 0.2655756032646429, + 0.5, + 0.7344243967353571, + 0.9151119481392835, + 1.0 + ] theta = [0.0 0.03284626432829264 0.018002223201815104 0.027529761904763195 0.02170954269630787 0.02464779425215191 0.023809523809520172; 0.0 0.059322894027551365 0.15770113064168867 0.1277882555598353 0.14409488824697547 0.13619592709181916 0.1384130236807124; 0.0 -0.01076859445118926 0.10235481204686203 0.23748565272164512 0.20629511050420746 0.2193616205757678 0.21587269060495373; @@ -185,6 +207,16 @@ function get_constant_parameters(alg::MPDeC) 0.0 0.002217096588914541 -0.005681864566224326 0.010624768120945982 -0.019288106960911655 0.07909012965322404 0.13841302368079056; 0.0 -0.000838270442615105 0.0020999811132187685 -0.0037202380952379155 0.005807300607708399 -0.00903674051876635 0.02380952380952195] elseif alg.M == 7 + nodes = [ + 0.0, + 0.06412992574519666, + 0.20414990928342885, + 0.3953503910487606, + 0.6046496089512394, + 0.7958500907165711, + 0.9358700742548034, + 1.0 + ] theta = [0.0 0.024737514438875514 0.013258719822130019 0.021034356725210923 0.01577972028613317 0.019036721343624663 0.017385669593011244 0.01785714285716722; 0.0 0.044892662602755 0.12064973282409763 0.09628017831208524 0.11096097521429726 0.10224716355590147 0.10657833455800869 0.10535211357202456; 0.0 -0.008140767742389747 0.07999763578665159 0.1890416498443277 0.16114433643933257 0.17539401516496778 0.1687159595472414 0.17056134624175545; @@ -194,6 +226,17 @@ function get_constant_parameters(alg::MPDeC) 0.0 -0.0012262209861064882 0.00310495001592085 -0.005608861642537821 0.00907193525967287 -0.015297619252294226 0.060459450968835426 0.10535211357171193; 0.0 0.0004714732640502682 -0.001179578486445107 0.002077422571009513 -0.003177213868071016 0.0045984230350075705 -0.006880371581776679 0.017857142857105046] elseif alg.M == 8 + nodes = [ + 0.0, + 0.05012100229426997, + 0.16140686024463113, + 0.3184412680869109, + 0.5, + 0.6815587319130891, + 0.8385931397553689, + 0.94987899770573, + 1.0 + ] theta = [0.0 0.019293838201043245 0.01018408040822739 0.01656936984357027 0.011990017361096061 0.01514127656190567 0.013175811859724718 0.01417408466352299 0.013888888888914153; 0.0 0.03512552097762184 0.09508650844936217 0.07509351709620726 0.08786983372308654 0.07945805621039881 0.08459518591985216 0.0820139081133675 0.08274768078126726; 0.0 -0.006364102418704803 0.0639319956284379 0.15288206102488378 0.12868222230659399 0.14236568211174472 0.13451403216049584 0.13834341694882824 0.13726935625072656; @@ -204,6 +247,18 @@ function get_constant_parameters(alg::MPDeC) 0.0 0.0007337726665392911 -0.0018475051393835457 0.0032896245700402282 -0.005122152942657721 0.007654163684208015 -0.012338827668713748 0.04762215980304063 0.08274768078103989; 0.0 -0.00028519577496630263 0.0007130770290413452 -0.0012523876730271555 0.0018988715277794554 -0.002680480954676767 0.0037048084806841075 -0.005404949312023177 0.013888888889184159] elseif alg.M == 9 + nodes = [ + 0.0, + 0.04023304591677057, + 0.13061306744724743, + 0.26103752509477773, + 0.4173605211668065, + 0.5826394788331936, + 0.7389624749052223, + 0.8693869325527526, + 0.9597669540832294, + 1.0 + ] theta = [0.0 0.015465148886122208 0.008074564828863144 0.013376018525479871 0.009424872116881033 0.012319010010855891 0.010311035538052238 0.01156730916181914 0.010928591594165482 0.011111111110039928; 0.0 0.02821797600073163 0.07677451467523239 0.060184542475466785 0.07119971523047752 0.06348355183104104 0.06872200531256567 0.06548282490166457 0.06711913714752882 0.0666529954141879; 0.0 -0.0051090759506765785 0.05211803626519543 0.12565307727248143 0.10482643149413295 0.11734393918791586 0.10937232662209873 0.11414511256730009 0.11177491077136281 0.1124446710354885; @@ -218,12 +273,13 @@ function get_constant_parameters(alg::MPDeC) error("Theta Gauss-Lobatto M > 9") end end - return theta + return nodes, theta end -struct MPDeCConstantCache{KType, T, T2} <: OrdinaryDiffEqConstantCache +struct MPDeCConstantCache{KType, NType, T, T2} <: OrdinaryDiffEqConstantCache K::KType M::KType + nodes::NType theta::T2 small_constant::T end @@ -237,15 +293,16 @@ function alg_cache(alg::MPDeC, u, rate_prototype, ::Type{uEltypeNoUnits}, throw(ArgumentError("MPDeC can only be applied to production-destruction systems")) end - theta = get_constant_parameters(alg) - MPDeCConstantCache(alg.K, alg.M, theta, alg.small_constant_function(uEltypeNoUnits)) + nodes, theta = get_constant_parameters(alg) + MPDeCConstantCache(alg.K, alg.M, nodes, theta, + alg.small_constant_function(uEltypeNoUnits)) end function initialize!(integrator, cache::MPDeCConstantCache) end # out-of-place -function build_mpdec_matrix(m, prod, C, p, t, dt, theta, small_constant) +function build_mpdec_matrix(m, prod, C, p, t, dt, nodes, theta, small_constant) N, M = size(C) M = M - 1 @@ -253,7 +310,7 @@ function build_mpdec_matrix(m, prod, C, p, t, dt, theta, small_constant) # P stores the production matrix in in-place computations. # In out-of-place computations it is not needed P = nothing - build_mpdec_matrix!(Mmat, m, prod, P, C, p, t, dt, theta, small_constant) + build_mpdec_matrix!(Mmat, m, prod, P, C, p, t, dt, nodes, theta, small_constant) if C isa StaticArray return SMatrix(Mmat) @@ -263,7 +320,7 @@ function build_mpdec_matrix(m, prod, C, p, t, dt, theta, small_constant) end # in-place for dense arrays -function build_mpdec_matrix!(Mmat, m, prod, P, C, p, t, dt, theta, small_constant) +function build_mpdec_matrix!(Mmat, m, prod, P, C, p, t, dt, nodes, theta, small_constant) N, M = size(C) M = M - 1 @@ -280,9 +337,9 @@ function build_mpdec_matrix!(Mmat, m, prod, P, C, p, t, dt, theta, small_constan dt_th = dt * th #TODO: This should be checked earlier and only once if isinplace(prod, 4) - prod(P, C[:, r], p, t) + prod(P, C[:, r], p, t + nodes[r] * dt) else - P = prod(C[:, r], p, t) + P = prod(C[:, r], p, t + nodes[r] * dt) end for i in 1:N for j in 1:N @@ -300,7 +357,7 @@ end @muladd function perform_step!(integrator, cache::MPDeCConstantCache, repeat_step = false) @unpack alg, t, dt, uprev, f, p = integrator - @unpack K, M, theta, small_constant = cache + @unpack K, M, nodes, theta, small_constant = cache C2 = zeros(length(uprev), M + 1) for i in 1:(M + 1) @@ -314,9 +371,9 @@ end C = copy(C2) end - for _ in 1:K + for k in 1:K for m in 2:(M + 1) - Mmat = build_mpdec_matrix(m, f.p, C, p, t, dt, theta, small_constant) + Mmat = build_mpdec_matrix(m, f.p, C, p, t, dt, nodes, theta, small_constant) # solve linear system linprob = LinearProblem(Mmat, uprev) @@ -329,7 +386,8 @@ end u2 = C[:, M + 1] #TODO: Remove this - #@assert u2 ≈ MPDeC_check(K, M, uprev, theta, f, p, t, dt) + #u1 = MPDeC_check(K, M, uprev, theta, f, p, t, dt) + #@assert u2 ≈ u1 integrator.u = u2 end @@ -366,8 +424,8 @@ function alg_cache(alg::MPDeC, u, rate_prototype, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck, ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} - theta = get_constant_parameters(alg) - tab = MPDeCConstantCache(alg.K, alg.M, theta, + nodes, theta = get_constant_parameters(alg) + tab = MPDeCConstantCache(alg.K, alg.M, nodes, theta, alg.small_constant_function(uEltypeNoUnits)) tmp = zero(u) @@ -509,7 +567,7 @@ end repeat_step = false) @unpack t, dt, uprev, u, f, p = integrator @unpack P, P2, σ, C, C2, linsolve_rhs, linsolve = cache - @unpack K, M, theta, small_constant = cache.tab + @unpack K, M, nodes, theta, small_constant = cache.tab # Set right hand side of linear system linsolve_rhs .= uprev @@ -522,7 +580,7 @@ end for _ in 1:K for m in 2:(M + 1) - build_mpdec_matrix!(P2, m, f.p, P, C, p, t, dt, theta, small_constant) + build_mpdec_matrix!(P2, m, f.p, P, C, p, t, dt, nodes, theta, small_constant) # Same as linres = P2 \ linsolve_rhs linsolve.A = P2 diff --git a/test/runtests.jl b/test/runtests.jl index aebe292d..dace8cbe 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -38,7 +38,7 @@ solution is computed using `ref_alg`. """ function experimental_orders_of_convergence(prob, alg, dts; test_time = nothing, only_first_index = false, - ref_alg = TRBDF2(autodiff = AutoFiniteDiff())) + ref_alg = Vern7()) @assert length(dts) > 1 errors = zeros(eltype(dts), length(dts)) @@ -1826,7 +1826,7 @@ end # Here we check the convergence order of pth-order schemes for which # also an interpolation of order p is available @testset "Convergence tests (nonconservative)" begin - #TODO: Add MPDeC - PDSProblem + #TODO: Add MPDeC(2) - PDSProblem algs = (MPE(), MPRK22(0.5), MPRK22(1.0), MPRK22(2.0), SSPMPRK22(0.5, 1.0)) dts = 0.5 .^ (4:15) problems = (prob_pds_linmod_nonconservative, @@ -1877,7 +1877,7 @@ end # Here we check the convergence order of pth-order schemes for which # no interpolation of order p is available @testset "Convergence tests (nonconservative)" begin - #TODO: Add MPDeC(3) - PDSProblem + #TODO: Add MPDeC(K), K≥ 3 - PDSProblem dts = 0.5 .^ (4:12) problems = (prob_pds_linmod_nonconservative, prob_pds_linmod_nonconservative_inplace) @@ -1890,10 +1890,15 @@ end end @testset "Interpolation tests (conservative)" begin - #TODO: Add MPDeC algs = (MPE(), MPRK22(0.5), MPRK22(1.0), MPRK22(2.0), MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), MPRK43II(0.5), MPRK43II(2.0 / 3.0), - SSPMPRK22(0.5, 1.0), SSPMPRK43()) + SSPMPRK22(0.5, 1.0), SSPMPRK43(), + MPDeC(2), MPDeC(2, nodes = :gausslobatto), + MPDeC(3), MPDeC(3, nodes = :gausslobatto), + MPDeC(4), MPDeC(4, nodes = :gausslobatto), + MPDeC(5), MPDeC(5, nodes = :gausslobatto), + MPDeC(6), MPDeC(6, nodes = :gausslobatto), + MPDeC(7), MPDeC(7, nodes = :gausslobatto)) dt = 0.5^6 problems = (prob_pds_linmod, prob_pds_linmod_array, prob_pds_linmod_mvector, prob_pds_linmod_inplace) @@ -1909,11 +1914,11 @@ end end @testset "Check convergence order (nonautonomous conservative PDS)" begin - #TODO: Add MPDeC + #TODO: Check MPDeC(K) order for K ≥ 4 prod! = (P, u, p, t) -> begin fill!(P, zero(eltype(P))) P[1, 2] = sin(t)^2 * u[2] - P[2, 1] = cos(2 * t)^2 * u[1] + P[2, 1] = cos(t)^2 * u[1] return nothing end prod = (u, p, t) -> begin @@ -1921,14 +1926,15 @@ end prod!(P, u, p, t) return P end - u0 = [1.0; 0.0] + u0 = [1.1; 0.9] # values close to zero may decrease the order tspan = (0.0, 1.0) prob_oop = ConservativePDSProblem(prod, u0, tspan) #out-of-place prob_ip = ConservativePDSProblem(prod!, u0, tspan) #in-place dts = 0.5 .^ (4:15) algs = (MPE(), MPRK22(0.5), MPRK22(1.0), MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), - MPRK43II(2.0 / 3.0), MPRK43II(0.5), SSPMPRK22(0.5, 1.0), SSPMPRK43()) + MPRK43II(2.0 / 3.0), MPRK43II(0.5), SSPMPRK22(0.5, 1.0), SSPMPRK43(), + MPDeC(2), MPDeC(2, nodes = :lagrange), MPDeC(3), MPDeC(3, nodes = :lagrange)) @testset "$alg" for alg in algs orders = experimental_orders_of_convergence(prob_oop, alg, dts) @test check_order(orders, PositiveIntegrators.alg_order(alg), atol = 0.2) From b9bdc019fe90d3ba54de5275c4b8750f8a3bc506 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Wed, 12 Mar 2025 14:49:02 +0100 Subject: [PATCH 08/38] added adaptive MPDeC --- src/mpdec.jl | 56 +++++++++++++++++++++++++++++------------------- test/runtests.jl | 21 ++++++++++++++---- 2 files changed, 51 insertions(+), 26 deletions(-) diff --git a/src/mpdec.jl b/src/mpdec.jl index 51c1caa7..5a887550 100644 --- a/src/mpdec.jl +++ b/src/mpdec.jl @@ -32,7 +32,7 @@ to avoid divisions by zero. You can pass a value explicitly, otherwise `small_co "Arbitrary high-order, conservative and positivity preserving Patankar-type deferred correction schemes." Applied Numerical Mathematics 153 (2020): 15-34. """ -struct MPDeC{T, N, F, T2} <: OrdinaryDiffEqAlgorithm +struct MPDeC{T, N, F, T2} <: OrdinaryDiffEqAdaptiveAlgorithm K::T M::T nodes::N @@ -359,19 +359,14 @@ end @unpack alg, t, dt, uprev, f, p = integrator @unpack K, M, nodes, theta, small_constant = cache + C = zeros(length(uprev), M + 1) C2 = zeros(length(uprev), M + 1) for i in 1:(M + 1) C2[:, i] = uprev end - if uprev isa StaticArray - #TODO make C static - C = copy(C2) - else - C = copy(C2) - end - - for k in 1:K + for _ in 1:K + C .= C2 for m in 2:(M + 1) Mmat = build_mpdec_matrix(m, f.p, C, p, t, dt, nodes, theta, small_constant) @@ -381,15 +376,22 @@ end C2[:, m] = sol.u integrator.stats.nsolve += 1 end - C = copy(C2) + #C = copy(C2) end - u2 = C[:, M + 1] + u = C2[:, M + 1] + u1 = C[:, M + 1] # one order less accurate - #TODO: Remove this - #u1 = MPDeC_check(K, M, uprev, theta, f, p, t, dt) - #@assert u2 ≈ u1 + #TODO: Remove this check + # Check only valid for autonomous conserverative PDS + #u_check = MPDeC_check(K, M, uprev, theta, f, p, t, dt) + #@assert u ≈ u_check - integrator.u = u2 + tmp = u - u1 + atmp = calculate_residuals(tmp, uprev, u, integrator.opts.abstol, + integrator.opts.reltol, integrator.opts.internalnorm, t) + integrator.EEst = integrator.opts.internalnorm(atmp, t) + + integrator.u = u end struct MPDeCCache{uType, PType, CType, tabType, F} <: MPRKMutableCache @@ -407,6 +409,7 @@ struct MPDeCCache{uType, PType, CType, tabType, F} <: MPRKMutableCache end struct MPDeCConservativeCache{uType, PType, CType, tabType, F} <: MPRKMutableCache + tmp::uType P::PType P2::PType σ::uType @@ -445,7 +448,7 @@ function alg_cache(alg::MPDeC, u, rate_prototype, ::Type{uEltypeNoUnits}, linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true, assumptions = LinearSolve.OperatorAssumptions(true)) - MPDeCConservativeCache(P, P2, σ, C, C2, + MPDeCConservativeCache(tmp, P, P2, σ, C, C2, tab, #MPDeCConstantCache linsolve_rhs, linsolve) @@ -566,19 +569,19 @@ end @muladd function perform_step!(integrator, cache::MPDeCConservativeCache, repeat_step = false) @unpack t, dt, uprev, u, f, p = integrator - @unpack P, P2, σ, C, C2, linsolve_rhs, linsolve = cache + @unpack tmp, P, P2, σ, C, C2, linsolve_rhs, linsolve = cache @unpack K, M, nodes, theta, small_constant = cache.tab # Set right hand side of linear system linsolve_rhs .= uprev - # Initializ C matrices + # Initialize C matrices for i in 1:(M + 1) - C2[:, i] = uprev - C[:, i] = uprev + C2[:, i] .= uprev end for _ in 1:K + C .= C2 for m in 2:(M + 1) build_mpdec_matrix!(P2, m, f.p, P, C, p, t, dt, nodes, theta, small_constant) @@ -588,10 +591,19 @@ end C2[:, m] .= linres integrator.stats.nsolve += 1 end - C .= C2 end - u .= C[:, M + 1] + u .= C2[:, M + 1] + σ .= C[:, M + 1] # one order less accurate + + # Now σ stores the error estimate + @.. broadcast=false σ=u - σ + + # Now tmp stores error residuals + calculate_residuals!(tmp, σ, uprev, u, integrator.opts.abstol, + integrator.opts.reltol, integrator.opts.internalnorm, t, + False()) + integrator.EEst = integrator.opts.internalnorm(tmp, t) end ######################################################################################################## diff --git a/test/runtests.jl b/test/runtests.jl index dace8cbe..575e3dd9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1195,7 +1195,6 @@ end # Here we check that MPRK22(α) = SSPMPRK22(0,α) @testset "MPRK22(α) = SSPMPRK22(0, α)" begin - #TODO: Add MPDeC(2) and MPDeC(2, nodes=:lagrange) REQUIRES ADAPTIVITY for α in (0.5, 2.0 / 3.0, 1.0, 2.0) # conservative PDS sol1 = solve(prob_pds_linmod, MPRK22(α)) @@ -1213,6 +1212,19 @@ end end end + # Here we check that MPRK22(1.0) = MPDeC(2) + @testset "MPRK22(1.0) = MPDeC(2)" begin + #TODO: Add MPDeC(2) for nonconservative PDS - REQUIRES PDSProblem + # conservative PDS + sol1 = solve(prob_pds_linmod, MPRK22(1.0)) + sol2 = solve(prob_pds_linmod, MPDeC(2)) + sol3 = solve(prob_pds_linmod, MPDeC(2, nodes = :lagrange)) + sol4 = solve(prob_pds_linmod_inplace, MPRK22(1.0)) + sol5 = solve(prob_pds_linmod_inplace, MPDeC(2)) + sol6 = solve(prob_pds_linmod_inplace, MPDeC(2, nodes = :lagrange)) + @test sol1.u ≈ sol2.u ≈ sol3.u ≈ sol4.u ≈ sol5.u ≈ sol6.u + end + # Here we check that different linear solvers can be used @testset "Different linear solvers" begin #TODO: Add MPDeC - REQUIRES PDSProblem @@ -1361,7 +1373,7 @@ end end @testset "Different matrix types (conservative, adaptive)" begin - #TODO: SSPMPRK43 is not adaptive and needs to be removed + #TODO: MPE, SSPMPRK43 are not adaptive and need to be removed (solve without giving dt!) #TODO: Add MPDeC prod_1! = (P, u, p, t) -> begin fill!(P, zero(eltype(P))) @@ -1575,7 +1587,7 @@ end end @testset "Different matrix types (nonconservative, adaptive)" begin - #TODO: SSPMPRK43 is not adaptive and needs to be removed + #TODO: MPE, SSPMPRK43 are not adaptive and need to be removed (solve without giving dt) #TODO: Add MPDeC - Requires PDSProblem and adaptivity prod_1! = (P, u, p, t) -> begin fill!(P, zero(eltype(P))) @@ -1934,7 +1946,8 @@ end dts = 0.5 .^ (4:15) algs = (MPE(), MPRK22(0.5), MPRK22(1.0), MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), MPRK43II(2.0 / 3.0), MPRK43II(0.5), SSPMPRK22(0.5, 1.0), SSPMPRK43(), - MPDeC(2), MPDeC(2, nodes = :lagrange), MPDeC(3), MPDeC(3, nodes = :lagrange)) + MPDeC(2), MPDeC(2, nodes = :lagrange), MPDeC(3), + MPDeC(3, nodes = :lagrange)) @testset "$alg" for alg in algs orders = experimental_orders_of_convergence(prob_oop, alg, dts) @test check_order(orders, PositiveIntegrators.alg_order(alg), atol = 0.2) From 280cffaeb4d1a5905a0945b1ea436e2a1143328c Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Wed, 12 Mar 2025 15:52:39 +0100 Subject: [PATCH 09/38] mpdec interpolation --- src/interpolation.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/interpolation.jl b/src/interpolation.jl index e86cdf84..d257009d 100644 --- a/src/interpolation.jl +++ b/src/interpolation.jl @@ -46,7 +46,8 @@ const MPRKCaches = Union{MPEConstantCache, MPECache, MPEConservativeCache, MPRK22ConstantCache, MPRK22Cache, MPRK22ConservativeCache, MPRK43ConstantCache, MPRK43Cache, MPRK43ConservativeCache, SSPMPRK22ConstantCache, SSPMPRK22Cache, SSPMPRK22ConservativeCache, - SSPMPRK43ConstantCache, SSPMPRK43Cache, SSPMPRK43ConservativeCache} + SSPMPRK43ConstantCache, SSPMPRK43Cache, SSPMPRK43ConservativeCache, + MPDeCConstantCache, MPDeCCache, MPDeCConservativeCache} function interp_summary(::Type{cacheType}, dense::Bool) where {cacheType <: MPRKCaches} From d417f726ff679e5bfc16ddd06aec94bdee2ff071 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Wed, 12 Mar 2025 16:55:36 +0100 Subject: [PATCH 10/38] MPDeC for PDSProblem (out-of-place) --- src/mpdec.jl | 47 ++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 38 insertions(+), 9 deletions(-) diff --git a/src/mpdec.jl b/src/mpdec.jl index 5a887550..9f6735a3 100644 --- a/src/mpdec.jl +++ b/src/mpdec.jl @@ -302,25 +302,35 @@ function initialize!(integrator, cache::MPDeCConstantCache) end # out-of-place -function build_mpdec_matrix(m, prod, C, p, t, dt, nodes, theta, small_constant) +function build_mpdec_matrix(uprev, m, f, C, p, t, dt, nodes, theta, small_constant) N, M = size(C) M = M - 1 Mmat = zeros(eltype(C), N, N) + rhs = similar(uprev) + rhs .= uprev + # P stores the production matrix in in-place computations. # In out-of-place computations it is not needed P = nothing - build_mpdec_matrix!(Mmat, m, prod, P, C, p, t, dt, nodes, theta, small_constant) + if f isa PDSFunction + build_mpdec_matrix!(Mmat, rhs, m, f.p, P, C, p, t, dt, nodes, theta, + small_constant, f.d) + else + build_mpdec_matrix!(Mmat, rhs, m, f.p, P, C, p, t, dt, nodes, theta, small_constant) + end if C isa StaticArray - return SMatrix(Mmat) + return SMatrix(Mmat), SVector(rhs) else - return Mmat + return Mmat, rhs end end # in-place for dense arrays -function build_mpdec_matrix!(Mmat, m, prod, P, C, p, t, dt, nodes, theta, small_constant) +function build_mpdec_matrix!(Mmat, rhs, m, prod, P, C, p, t, dt, nodes, theta, + small_constant, + dest = nothing) N, M = size(C) M = M - 1 @@ -338,10 +348,27 @@ function build_mpdec_matrix!(Mmat, m, prod, P, C, p, t, dt, nodes, theta, small_ #TODO: This should be checked earlier and only once if isinplace(prod, 4) prod(P, C[:, r], p, t + nodes[r] * dt) + #if !isnothing(dest) + # dest(D, C[:, r], p, t + nodes[r] * dt) + #end else P = prod(C[:, r], p, t + nodes[r] * dt) + if !isnothing(dest) + d = dest(C[:, r], p, t + nodes[r] * dt) + end end + for i in 1:N + # Add nonconservative destruction terms to diagonal (PDSFunctions only!) + if !isnothing(dest) + if th >= 0 + Mmat[i, i] += dt_th * d[i] / σ[i] + rhs[i] += dt_th * P[i, i] + else + Mmat[i, i] -= dt_th * P[i, i] / σ[i] + rhs[i] -= dt_th * d[i] + end + end for j in 1:N if th >= 0 Mmat[i, j] -= dt_th * P[i, j] / σ[j] @@ -368,10 +395,11 @@ end for _ in 1:K C .= C2 for m in 2:(M + 1) - Mmat = build_mpdec_matrix(m, f.p, C, p, t, dt, nodes, theta, small_constant) + Mmat, rhs = build_mpdec_matrix(uprev, m, f, C, p, t, dt, nodes, theta, + small_constant) # solve linear system - linprob = LinearProblem(Mmat, uprev) + linprob = LinearProblem(Mmat, rhs) sol = solve(linprob, alg.linsolve) C2[:, m] = sol.u integrator.stats.nsolve += 1 @@ -572,7 +600,7 @@ end @unpack tmp, P, P2, σ, C, C2, linsolve_rhs, linsolve = cache @unpack K, M, nodes, theta, small_constant = cache.tab - # Set right hand side of linear system + # Initialize right hand side of linear system linsolve_rhs .= uprev # Initialize C matrices @@ -583,7 +611,8 @@ end for _ in 1:K C .= C2 for m in 2:(M + 1) - build_mpdec_matrix!(P2, m, f.p, P, C, p, t, dt, nodes, theta, small_constant) + build_mpdec_matrix!(P2, linsolve_rhs, m, f.p, P, C, p, t, dt, nodes, theta, + small_constant) # Same as linres = P2 \ linsolve_rhs linsolve.A = P2 From 6bb4bfc5a6b863c2d5ee9bfabdcd41c55c23e057 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Thu, 13 Mar 2025 16:23:25 +0100 Subject: [PATCH 11/38] added mpdec for non-conservative PDS (in-place) --- src/mpdec.jl | 153 ++++++++++++++++------------------------------- test/runtests.jl | 73 +++++++++++++--------- 2 files changed, 99 insertions(+), 127 deletions(-) diff --git a/src/mpdec.jl b/src/mpdec.jl index 9f6735a3..61e5e1a1 100644 --- a/src/mpdec.jl +++ b/src/mpdec.jl @@ -40,8 +40,17 @@ struct MPDeC{T, N, F, T2} <: OrdinaryDiffEqAdaptiveAlgorithm small_constant_function::T2 end +function small_constant_function_MPDeC(type) + if type == Float64 + small_constant = 1e-200 + else + small_constant = floatmin(type) + end + return small_constant +end + function MPDeC(K; nodes = :gausslobatto, linsolve = LUFactorization(), - small_constant = nothing) + small_constant = small_constant_function_MPDeC) if !(isinteger(K)) throw(ArgumentError("MPDeC requires the parameter K to be an integer.")) end @@ -49,9 +58,7 @@ function MPDeC(K; nodes = :gausslobatto, linsolve = LUFactorization(), K = Int(K) end - if isnothing(small_constant) - small_constant_function = floatmin - elseif small_constant isa Number + if small_constant isa Number small_constant_function = Returns(small_constant) else # assume small_constant isa Function small_constant_function = small_constant @@ -302,7 +309,7 @@ function initialize!(integrator, cache::MPDeCConstantCache) end # out-of-place -function build_mpdec_matrix(uprev, m, f, C, p, t, dt, nodes, theta, small_constant) +function build_mpdec_matrix_and_rhs(uprev, m, f, C, p, t, dt, nodes, theta, small_constant) N, M = size(C) M = M - 1 @@ -314,10 +321,11 @@ function build_mpdec_matrix(uprev, m, f, C, p, t, dt, nodes, theta, small_consta # In out-of-place computations it is not needed P = nothing if f isa PDSFunction - build_mpdec_matrix!(Mmat, rhs, m, f.p, P, C, p, t, dt, nodes, theta, - small_constant, f.d) + build_mpdec_matrix_and_rhs!(Mmat, rhs, m, f.p, P, C, p, t, dt, nodes, theta, + small_constant, f.d) else - build_mpdec_matrix!(Mmat, rhs, m, f.p, P, C, p, t, dt, nodes, theta, small_constant) + build_mpdec_matrix_and_rhs!(Mmat, rhs, m, f.p, P, C, p, t, dt, nodes, theta, + small_constant) end if C isa StaticArray @@ -328,9 +336,9 @@ function build_mpdec_matrix(uprev, m, f, C, p, t, dt, nodes, theta, small_consta end # in-place for dense arrays -function build_mpdec_matrix!(Mmat, rhs, m, prod, P, C, p, t, dt, nodes, theta, - small_constant, - dest = nothing) +function build_mpdec_matrix_and_rhs!(Mmat, rhs, m, prod, P, C, p, t, dt, nodes, theta, + small_constant, + dest = nothing, d = nothing) N, M = size(C) M = M - 1 @@ -340,6 +348,7 @@ function build_mpdec_matrix!(Mmat, rhs, m, prod, P, C, p, t, dt, nodes, theta, Mmat[i, i] = oneMmat end + #TODO: use predefined memory in inplace computations σ = add_small_constant(C[:, m], small_constant) for r in 1:(M + 1) @@ -348,9 +357,9 @@ function build_mpdec_matrix!(Mmat, rhs, m, prod, P, C, p, t, dt, nodes, theta, #TODO: This should be checked earlier and only once if isinplace(prod, 4) prod(P, C[:, r], p, t + nodes[r] * dt) - #if !isnothing(dest) - # dest(D, C[:, r], p, t + nodes[r] * dt) - #end + if !isnothing(dest) + dest(d, C[:, r], p, t + nodes[r] * dt) + end else P = prod(C[:, r], p, t + nodes[r] * dt) if !isnothing(dest) @@ -395,8 +404,8 @@ end for _ in 1:K C .= C2 for m in 2:(M + 1) - Mmat, rhs = build_mpdec_matrix(uprev, m, f, C, p, t, dt, nodes, theta, - small_constant) + Mmat, rhs = build_mpdec_matrix_and_rhs(uprev, m, f, C, p, t, dt, nodes, theta, + small_constant) # solve linear system linprob = LinearProblem(Mmat, rhs) @@ -404,7 +413,6 @@ end C2[:, m] = sol.u integrator.stats.nsolve += 1 end - #C = copy(C2) end u = C2[:, M + 1] u1 = C[:, M + 1] # one order less accurate @@ -426,8 +434,7 @@ struct MPDeCCache{uType, PType, CType, tabType, F} <: MPRKMutableCache tmp::uType P::PType P2::PType - D::uType - D2::uType + d::uType σ::uType C::CType C2::CType @@ -462,6 +469,7 @@ function alg_cache(alg::MPDeC, u, rate_prototype, ::Type{uEltypeNoUnits}, tmp = zero(u) P = p_prototype(u, f) # stores evaluation of the production matrix P2 = p_prototype(u, f) # stores the linear system matrix + d = zero(u) σ = zero(u) C = zeros(eltype(u), length(u), alg.M + 1) C2 = zeros(eltype(u), length(u), alg.M + 1) @@ -485,11 +493,7 @@ function alg_cache(alg::MPDeC, u, rate_prototype, ::Type{uEltypeNoUnits}, linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true, assumptions = LinearSolve.OperatorAssumptions(true)) - MPDeCCache(tmp, P, P2, - similar(u), # D - similar(u), # D2 - σ, - C, C2, + MPDeCCache(tmp, P, P2, d, σ, C, C2, tab, #MPDeCConstantCache linsolve_rhs, linsolve) @@ -501,89 +505,36 @@ end function initialize!(integrator, cache::Union{MPDeCCache, MPDeCConservativeCache}) end -#= COPY OF MPRK22 @muladd function perform_step!(integrator, cache::MPDeCCache, repeat_step = false) @unpack t, dt, uprev, u, f, p = integrator - @unpack tmp, P, P2, D, D2, σ, linsolve = cache - @unpack a21, b1, b2, small_constant = cache.tab - - # We use P2 to store the last evaluation of the PDS - # as well as to store the system matrix of the linear system - - f.p(P, uprev, p, t) # evaluate production terms - f.d(D, uprev, p, t) # evaluate nonconservative destruction terms - integrator.stats.nf += 1 - if issparse(P) - # We need to keep the structural nonzeros of the production terms. - # However, this is not guaranteed by broadcasting, see - # https://github.com/JuliaSparse/SparseArrays.jl/issues/190 - nz_P = nonzeros(P) - nz_P2 = nonzeros(P2) - @.. broadcast=false nz_P2=a21 * nz_P - else - @.. broadcast=false P2=a21 * P - end - @.. broadcast=false D2=a21 * D - - # avoid division by zero due to zero Patankar weights - @.. broadcast=false σ=uprev + small_constant + @unpack tmp, P, P2, d, σ, C, C2, linsolve_rhs, linsolve = cache + @unpack K, M, nodes, theta, small_constant = cache.tab - # tmp holds the right hand side of the linear system - tmp .= uprev - @inbounds for i in eachindex(tmp) - tmp[i] += dt * P2[i, i] + # Initialize C matrices + for i in 1:(M + 1) + C2[:, i] .= uprev end - build_mprk_matrix!(P2, P2, σ, dt, D2) - - # Same as linres = P2 \ tmp - linsolve.A = P2 - linres = solve!(linsolve) - - u .= linres - integrator.stats.nsolve += 1 - - if isone(a21) - σ .= u - else - @.. broadcast=false σ=σ^(1 - 1 / a21) * u^(1 / a21) - end - @.. broadcast=false σ=σ + small_constant - - f.p(P2, u, p, t + a21 * dt) # evaluate production terms - f.d(D2, u, p, t + a21 * dt) # evaluate nonconservative destruction terms - integrator.stats.nf += 1 - - if issparse(P) - # We need to keep the structural nonzeros of the production terms. - # However, this is not guaranteed by broadcasting, see - # https://github.com/JuliaSparse/SparseArrays.jl/issues/190 - nz_P = nonzeros(P) - nz_P2 = nonzeros(P2) - @.. broadcast=false nz_P2=b1 * nz_P + b2 * nz_P2 - else - @.. broadcast=false P2=b1 * P + b2 * P2 - end - @.. broadcast=false D2=b1 * D + b2 * D2 + for _ in 1:K + C .= C2 + for m in 2:(M + 1) + linsolve_rhs .= uprev + build_mpdec_matrix_and_rhs!(P2, linsolve_rhs, m, f.p, P, C, p, t, dt, nodes, + theta, + small_constant, f.d, d) - # tmp holds the right hand side of the linear system - tmp .= uprev - @inbounds for i in eachindex(tmp) - tmp[i] += dt * P2[i, i] + # Same as linres = P2 \ linsolve_rhs + linsolve.A = P2 + linres = solve!(linsolve) + C2[:, m] .= linres + integrator.stats.nsolve += 1 + end end - build_mprk_matrix!(P2, P2, σ, dt, D2) - - # Same as linres = P2 \ tmp - linsolve.A = P2 - linres = solve!(linsolve) - - u .= linres - integrator.stats.nsolve += 1 + u .= C2[:, M + 1] + σ .= C[:, M + 1] # one order less accurate # Now σ stores the error estimate - # If a21 = 1, then σ is the MPE approximation, i.e. suited for stiff problems. - # Otherwise, this is not clear. @.. broadcast=false σ=u - σ # Now tmp stores error residuals @@ -592,7 +543,6 @@ end False()) integrator.EEst = integrator.opts.internalnorm(tmp, t) end -=# @muladd function perform_step!(integrator, cache::MPDeCConservativeCache, repeat_step = false) @@ -611,8 +561,10 @@ end for _ in 1:K C .= C2 for m in 2:(M + 1) - build_mpdec_matrix!(P2, linsolve_rhs, m, f.p, P, C, p, t, dt, nodes, theta, - small_constant) + linsolve_rhs .= uprev + build_mpdec_matrix_and_rhs!(P2, linsolve_rhs, m, f.p, P, C, p, t, dt, nodes, + theta, + small_constant) # Same as linres = P2 \ linsolve_rhs linsolve.A = P2 @@ -635,6 +587,7 @@ end integrator.EEst = integrator.opts.internalnorm(tmp, t) end +# TODO: Remove the code below after testing is complete ######################################################################################################## ######################################################################################################## #### The functions below are provided with the original paper and used for validation. #### diff --git a/test/runtests.jl b/test/runtests.jl index 575e3dd9..eece2cd1 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1214,7 +1214,6 @@ end # Here we check that MPRK22(1.0) = MPDeC(2) @testset "MPRK22(1.0) = MPDeC(2)" begin - #TODO: Add MPDeC(2) for nonconservative PDS - REQUIRES PDSProblem # conservative PDS sol1 = solve(prob_pds_linmod, MPRK22(1.0)) sol2 = solve(prob_pds_linmod, MPDeC(2)) @@ -1223,11 +1222,20 @@ end sol5 = solve(prob_pds_linmod_inplace, MPDeC(2)) sol6 = solve(prob_pds_linmod_inplace, MPDeC(2, nodes = :lagrange)) @test sol1.u ≈ sol2.u ≈ sol3.u ≈ sol4.u ≈ sol5.u ≈ sol6.u + + # nonconservative PDS + sol1 = solve(prob_pds_linmod_nonconservative, MPRK22(1.0)) + sol2 = solve(prob_pds_linmod_nonconservative, MPDeC(2)) + sol3 = solve(prob_pds_linmod_nonconservative, MPDeC(2, nodes = :lagrange)) + sol4 = solve(prob_pds_linmod_nonconservative_inplace, MPRK22(1.0)) + sol5 = solve(prob_pds_linmod_nonconservative_inplace, MPDeC(2)) + sol6 = solve(prob_pds_linmod_nonconservative_inplace, + MPDeC(2, nodes = :lagrange)) + @test sol1.u ≈ sol2.u ≈ sol3.u ≈ sol4.u ≈ sol5.u ≈ sol6.u end # Here we check that different linear solvers can be used @testset "Different linear solvers" begin - #TODO: Add MPDeC - REQUIRES PDSProblem # problem data u0 = [0.9, 0.1] tspan = (0.0, 2.0) @@ -1257,15 +1265,19 @@ end prob_ip_2 = ConservativePDSProblem(linmodP!, u0, tspan, p; analytic = f_analytic) - algs = (MPE, (; kwargs...) -> MPRK22(1.0; kwargs...), - (; kwargs...) -> MPRK22(0.5; kwargs...), - (; kwargs...) -> MPRK22(2.0; kwargs...), - (; kwargs...) -> MPRK43I(1.0, 0.5; kwargs...), - (; kwargs...) -> MPRK43I(0.5, 0.75; kwargs...), - (; kwargs...) -> MPRK43II(0.5; kwargs...), - (; kwargs...) -> MPRK43II(2.0 / 3.0; kwargs...), - (; kwargs...) -> SSPMPRK22(0.5, 1.0; kwargs...), - (; kwargs...) -> SSPMPRK43(; kwargs...)) + algs = [MPE, (; kwargs...) -> MPRK22(1.0; kwargs...), + (; kwargs...) -> MPRK22(0.5; kwargs...), + (; kwargs...) -> MPRK22(2.0; kwargs...), + (; kwargs...) -> MPRK43I(1.0, 0.5; kwargs...), + (; kwargs...) -> MPRK43I(0.5, 0.75; kwargs...), + (; kwargs...) -> MPRK43II(0.5; kwargs...), + (; kwargs...) -> MPRK43II(2.0 / 3.0; kwargs...), + (; kwargs...) -> SSPMPRK22(0.5, 1.0; kwargs...), + (; kwargs...) -> SSPMPRK43(; kwargs...)] + for k in 2:9 + push!(algs, (; kwargs...) -> MPDeC(k; kwargs...), + (; kwargs...) -> MPDeC(k; nodes = :lagrange, kwargs...)) + end for alg in algs # Check different linear solvers @@ -1838,13 +1850,12 @@ end # Here we check the convergence order of pth-order schemes for which # also an interpolation of order p is available @testset "Convergence tests (nonconservative)" begin - #TODO: Add MPDeC(2) - PDSProblem - algs = (MPE(), MPRK22(0.5), MPRK22(1.0), MPRK22(2.0), SSPMPRK22(0.5, 1.0)) + algs = (MPE(), MPRK22(0.5), MPRK22(1.0), MPRK22(2.0), SSPMPRK22(0.5, 1.0), + MPDeC(2), MPDeC(2; nodes = :lagrange)) dts = 0.5 .^ (4:15) problems = (prob_pds_linmod_nonconservative, prob_pds_linmod_nonconservative_inplace) @testset "$alg" for alg in algs - alg = MPRK22(1.0) for prob in problems orders = experimental_orders_of_convergence(prob, alg, dts) @test check_order(orders, PositiveIntegrators.alg_order(alg)) @@ -1858,13 +1869,13 @@ end dts; test_time) @test check_order(orders, PositiveIntegrators.alg_order(alg), - atol = 0.2) + atol = 0.3) orders = experimental_orders_of_convergence(prob, alg, dts; test_time, only_first_index = true) @test check_order(orders, PositiveIntegrators.alg_order(alg), - atol = 0.2) + N = 2, atol = 0.3) end end end @@ -1888,13 +1899,15 @@ end # Here we check the convergence order of pth-order schemes for which # no interpolation of order p is available + #TODO: Check order of MPDeC(K), K≥ 5 @testset "Convergence tests (nonconservative)" begin - #TODO: Add MPDeC(K), K≥ 3 - PDSProblem dts = 0.5 .^ (4:12) problems = (prob_pds_linmod_nonconservative, prob_pds_linmod_nonconservative_inplace) algs = (MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), - MPRK43II(0.5), MPRK43II(2.0 / 3.0), SSPMPRK43()) + MPRK43II(0.5), MPRK43II(2.0 / 3.0), SSPMPRK43(), + MPDeC(3), MPDeC(3; nodes = :lagrange), + MPDeC(4), MPDeC(4; nodes = :lagrange)) for alg in algs, prob in problems orders = experimental_orders_of_convergence(prob, alg, dts) @test check_order(orders, PositiveIntegrators.alg_order(alg), atol = 0.2) @@ -1925,8 +1938,8 @@ end end end + #TODO: Check MPDeC(K) order for K ≥ 4 @testset "Check convergence order (nonautonomous conservative PDS)" begin - #TODO: Check MPDeC(K) order for K ≥ 4 prod! = (P, u, p, t) -> begin fill!(P, zero(eltype(P))) P[1, 2] = sin(t)^2 * u[2] @@ -1956,8 +1969,8 @@ end end end + #TODO: Check order of MPDeC(K), K≥ 6 @testset "Check convergence order (nonautonomous nonconservative PDS)" begin - #TODO: Add MPDeC prod! = (P, u, p, t) -> begin fill!(P, zero(eltype(P))) P[1, 2] = sin(t)^2 * u[2] @@ -1981,14 +1994,18 @@ end dest!(d, u, p, t) return d end - u0 = [1.0; 0.0] + u0 = [1.1; 0.9] # stay away from zero tspan = (0.0, 1.0) prob_oop = PDSProblem(prod, dest, u0, tspan) #out-of-place prob_ip = PDSProblem(prod!, dest!, u0, tspan) #in-place - dts = 0.5 .^ (4:15) + dts = 0.5 .^ (4:10) algs = (MPE(), MPRK22(0.5), MPRK22(1.0), MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), - MPRK43II(2.0 / 3.0), MPRK43II(0.5), SSPMPRK22(0.5, 1.0), SSPMPRK43()) + MPRK43II(2.0 / 3.0), MPRK43II(0.5), SSPMPRK22(0.5, 1.0), SSPMPRK43(), + MPDeC(2), MPDeC(2; nodes = :lagrange), MPDeC(3), + MPDeC(3; nodes = :lagrange), + MPDeC(4), MPDeC(4; nodes = :lagrange), MPDeC(5), + MPDeC(5; nodes = :lagrange)) @testset "$alg" for alg in algs orders = experimental_orders_of_convergence(prob_oop, alg, dts) @test check_order(orders, PositiveIntegrators.alg_order(alg), atol = 0.2) @@ -2016,7 +2033,6 @@ end # Check that the schemes accept zero initial values @testset "Zero initial values" begin - #TODO: Add MPDeC - Requires PDSProblem # Do a single step and check that no NaNs occur u0 = [1.0, 0.0] dt = 1.0 @@ -2046,9 +2062,12 @@ end prob_oop = ConservativePDSProblem(prod, u0, tspan, p) prob_oop_2 = PDSProblem(prod, dest, u0, tspan, p) - algs = (MPE(), MPRK22(0.5), MPRK22(1.0), MPRK22(2.0), - MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), MPRK43II(0.5), - MPRK43II(2.0 / 3.0), SSPMPRK22(0.5, 1.0), SSPMPRK43()) + algs = [MPE(), MPRK22(0.5), MPRK22(1.0), MPRK22(2.0), + MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), MPRK43II(0.5), + MPRK43II(2.0 / 3.0), SSPMPRK22(0.5, 1.0), SSPMPRK43()] + for k in 2:9 + push!(algs, MPDeC(k), MPDeC(k; nodes = :lagrange)) + end for alg in algs sol = solve(prob_ip, alg; dt = dt, adaptive = false) From 9861f79488a4c3d4e94cb0a5fcb9d8dedd523349 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Fri, 14 Mar 2025 13:54:56 +0100 Subject: [PATCH 12/38] bugfix MPDeC for PDSProblem & additional tests --- src/mpdec.jl | 37 ++++++++++++++++++++++++---------- test/runtests.jl | 52 ++++++++++++++++++++++++++++++++---------------- 2 files changed, 61 insertions(+), 28 deletions(-) diff --git a/src/mpdec.jl b/src/mpdec.jl index 61e5e1a1..e29dcece 100644 --- a/src/mpdec.jl +++ b/src/mpdec.jl @@ -1,8 +1,6 @@ """ MPDeC(K; [nodes = :gausslobatto, linsolve = ..., small_constant = ...]) -EVERYTHING BELOW IS COPIED FROM MPRK22 - A family of arbitrary order modified Patankar-Runge-Kutta algorithms for production-destruction systems. Each member of this family is an adaptive, one-step method which is Kth order accurate, unconditionally positivity-preserving, and linearly @@ -31,6 +29,7 @@ to avoid divisions by zero. You can pass a value explicitly, otherwise `small_co - Davide Torlo, and Philipp Öffner. "Arbitrary high-order, conservative and positivity preserving Patankar-type deferred correction schemes." Applied Numerical Mathematics 153 (2020): 15-34. +- TODO: Literature non-autonomous DeC """ struct MPDeC{T, N, F, T2} <: OrdinaryDiffEqAdaptiveAlgorithm K::T @@ -42,7 +41,8 @@ end function small_constant_function_MPDeC(type) if type == Float64 - small_constant = 1e-200 + #TODO: Adjust default value based on K + small_constant = floatmin(type) else small_constant = floatmin(type) end @@ -124,6 +124,7 @@ function get_constant_parameters(alg::MPDeC) 0.0 -0.013405848450491272 -0.009863945578231281 -0.011894132653058165 -0.009674981103547253 -0.014615221088391195 0.04183673469395899 0.2070023148147584; 0.0 0.001623913454270594 0.0012093726379440242 0.0014349489795916492 0.0012093726379429626 0.0016239134542663791 -1.0658141036401503e-14 0.04346064814814099] elseif alg.M == 8 + #= theta = [0.0 0.03685850005511468 0.035688932980599594 0.03594029017857134 0.0358289241622568 0.035914937444895045 0.035803571428594694 0.03605492862659787 0.034885361552085214; 0.0 0.15387641920194012 0.2012610229276906 0.19782924107142996 0.1990828924162642 0.19819740685622378 0.1992857142857929 0.1969121334875581 0.2076895943564523; 0.0 -0.15861283344356245 -0.046840828924162636 0.009592633928562577 0.0021516754850381403 0.0065018050044045594 0.0016071428572104196 0.011744309413188603 -0.03273368606460281; @@ -133,6 +134,16 @@ function get_constant_parameters(alg::MPDeC) 0.0 -0.044477995480599525 -0.03434082892416224 -0.03923549107142321 -0.03488536155200972 -0.04232631999563363 0.01410714285702852 0.1258791473759011 -0.03273368606702931; 0.0 0.010777460868606693 0.008403880070546516 0.009492187499996696 0.008606701940021111 0.009860353284821599 0.006428571428614305 0.053813175154459714 0.2076895943564523; 0.0 -0.0011695670745149895 -0.0009182098765432678 -0.0010295758928574872 -0.0009435626102302086 -0.0010549286265453262 -0.0008035714285732354 -0.0019731385031036552 0.03488536155197153] + =# + theta = [0 1070017/29030400 32377/907200 12881/358400 4063/113400 41705/1161216 401/11200 149527/4147200 989/28350; + 0 2233547/14515200 22823/113400 35451/179200 2822/14175 115075/580608 279/1400 408317/2073600 2944/14175; + 0 -2302297/14515200 -21247/453600 1719/179200 61/28350 3775/580608 9/5600 24353/2073600 -464/14175; + 0 2797679/14515200 15011/113400 39967/179200 4094/14175 159175/580608 403/1400 542969/2073600 5248/14175; + 0 -31457/181440 -2903/22680 -351/2240 -227/2835 -125/36288 -9/280 343/25920 -454/2835; + 0 1573169/14515200 9341/113400 17217/179200 1154/14175 85465/580608 333/1400 368039/2073600 5248/14175; + 0 -645607/14515200 -15577/453600 -7031/179200 -989/28350 -24575/580608 79/5600 261023/2073600 -464/14175; + 0 156437/14515200 953/113400 243/25600 122/14175 5725/580608 9/1400 111587/2073600 2944/14175; + 0 -33953/29030400 -119/129600 -369/358400 -107/113400 -175/165888 -9/11200 -8183/4147200 989/28350] elseif alg.M == 9 theta = [0.0 0.03188616071428564 0.031009210268469034 0.03117187499999874 0.031111111111105316 0.031149339236719698 0.031111111111073342 0.03117187499992724 0.031009210268393872 0.031886160714066136; 0.0 0.14467159330295903 0.18532725847540765 0.1828236607142859 0.18359396433471176 0.1831509191895293 0.18357142857178133 0.18292556155574857 0.18461297276007826 0.17568080357159488; @@ -157,7 +168,7 @@ function get_constant_parameters(alg::MPDeC) 0.0 0.007575105385869255 0.006322003099781336 0.006733969155847452 0.006509967398880434 0.0066988293985943415 0.006454545454857907 0.006958222269304315 0.005062262842329801 0.04054565746628214 0.17753594143141527; 0.0 -0.0006785849984634729 -0.0005679145956923939 -0.000603642451298736 -0.0005846828069051568 -0.0006001284755665637 -0.000581168831175205 -0.0006168966868074222 -0.0005062262839601317 -0.0011848112826555735 0.026834148362013366] else - error("Theta Lagrange M > 9") + error("MPDeC requires 2 ≤ K ≤ 9.") end else # alg.nodes == :gausslobatto if alg.M == 1 @@ -277,7 +288,7 @@ function get_constant_parameters(alg::MPDeC) 0.0 -0.000466141727447515 0.0011701705230945658 -0.002069009887745693 0.0031694435943515065 -0.004546719804944366 0.0064684529496542575 -0.010121519251002686 0.03843501942257066 0.0666529954255306; 0.0 0.0001825195178684848 -0.0004561980512006905 0.0008000755736511378 -0.0012078988997608064 0.00168623899423892 -0.002264907414485151 0.00303654628237382 -0.004354037773737218 0.011111111112313665] else - error("Theta Gauss-Lobatto M > 9") + error("MPDeC requires 2 ≤ K ≤ 9.") end end return nodes, theta @@ -321,9 +332,11 @@ function build_mpdec_matrix_and_rhs(uprev, m, f, C, p, t, dt, nodes, theta, smal # In out-of-place computations it is not needed P = nothing if f isa PDSFunction + # Additional destruction terms build_mpdec_matrix_and_rhs!(Mmat, rhs, m, f.p, P, C, p, t, dt, nodes, theta, small_constant, f.d) else + # No additional destruction terms build_mpdec_matrix_and_rhs!(Mmat, rhs, m, f.p, P, C, p, t, dt, nodes, theta, small_constant) end @@ -379,12 +392,14 @@ function build_mpdec_matrix_and_rhs!(Mmat, rhs, m, prod, P, C, p, t, dt, nodes, end end for j in 1:N - if th >= 0 - Mmat[i, j] -= dt_th * P[i, j] / σ[j] - Mmat[i, i] += dt_th * P[j, i] / σ[i] - else - Mmat[i, j] += dt_th * P[j, i] / σ[j] - Mmat[i, i] -= dt_th * P[i, j] / σ[i] + if j != i + if th >= 0 + Mmat[i, j] -= dt_th * P[i, j] / σ[j] + Mmat[i, i] += dt_th * P[j, i] / σ[i] # P[j, i] = D[i, j] + else + Mmat[i, j] += dt_th * P[j, i] / σ[j] # P[j, i] = D[i, j] + Mmat[i, i] -= dt_th * P[i, j] / σ[i] + end end end end diff --git a/test/runtests.jl b/test/runtests.jl index eece2cd1..286eff75 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2084,7 +2084,6 @@ end # Check that approximations, and thus the Patankar weights, # remain positive to avoid division by zero. @testset "Positvity check" begin - #TODO: Add MPDeC - Requires PDSProblem # For this problem u[1] decreases montonically to 0 very fast. # We perform 10^5 steps and check that u[end] does not contain any NaNs u0 = [0.9, 0.1] @@ -2114,9 +2113,12 @@ end prob_oop = ConservativePDSProblem(prod, u0, tspan, p) prob_oop_2 = PDSProblem(prod, dest, u0, tspan, p) - algs = (MPE(), MPRK22(0.5), MPRK22(1.0), MPRK22(2.0), - MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), MPRK43II(0.5), - MPRK43II(2.0 / 3.0), SSPMPRK22(0.5, 1.0), SSPMPRK43()) + algs = [MPE(), MPRK22(0.5), MPRK22(1.0), MPRK22(2.0), + MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), MPRK43II(0.5), + MPRK43II(2.0 / 3.0), SSPMPRK22(0.5, 1.0), SSPMPRK43()] + for k in 2:9 + push!(algs, MPDeC(k), MPDeC(k; nodes = :lagrange)) + end dt = 1e-3 for alg in algs @@ -2135,9 +2137,11 @@ end # Here we check that the implemented schemes can solve the predefined PDS # (at least for specific parameters) @testset "PDS problem library (adaptive schemes)" begin - #TODO: Add MPDeC - algs = (MPRK22(0.5), MPRK22(1.0), MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), - MPRK43II(2.0 / 3.0), MPRK43II(0.5), SSPMPRK22(0.5, 1.0)) + algs = [MPRK22(0.5), MPRK22(1.0), MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), + MPRK43II(2.0 / 3.0), MPRK43II(0.5), SSPMPRK22(0.5, 1.0)] + for k in 2:9 + push!(algs, MPDeC(k), MPDeC(k; nodes = :lagrange)) + end probs = (prob_pds_linmod, prob_pds_linmod_inplace, prob_pds_nonlinmod, prob_pds_robertson, prob_pds_bertolazzi, prob_pds_brusselator, prob_pds_npzd, @@ -2151,6 +2155,9 @@ end elseif prob == prob_pds_stratreac && alg == MPRK43I(0.5, 0.75) # Not successful on Julia 1.9 break + elseif prob == prob_pds_stratreac && alg == MPDeC(9; nodes = :lagrange) + # unstable + break end # later versions of OrdinaryDiffEq.jl use dtmin = 0 by default, # see https://github.com/SciML/OrdinaryDiffEq.jl/pull/2098 @@ -2179,7 +2186,6 @@ end #Here we check the different possibilities to define small_constant. @testset "Different possibilities to set small_constant" begin - #TODO: Add MPDeC - Requires PDSProblem # For this problem u[1] decreases montonically to 0 very fast. u0 = [0.9, 0.1] tspan = (0.0, 100.0) @@ -2213,11 +2219,15 @@ end prob_pds_npzd, prob_pds_sir, prob_pds_stratreac) - algs = (MPE, (; kwargs...) -> MPRK22(1.0; kwargs...), - (; kwargs...) -> MPRK43I(1.0, 0.5; kwargs...), - (; kwargs...) -> MPRK43II(0.5; kwargs...), - (; kwargs...) -> SSPMPRK22(0.5, 1.0; kwargs...), - (; kwargs...) -> SSPMPRK43(; kwargs...)) + algs = [MPE, (; kwargs...) -> MPRK22(1.0; kwargs...), + (; kwargs...) -> MPRK43I(1.0, 0.5; kwargs...), + (; kwargs...) -> MPRK43II(0.5; kwargs...), + (; kwargs...) -> SSPMPRK22(0.5, 1.0; kwargs...), + (; kwargs...) -> SSPMPRK43(; kwargs...)] + for k in 2:9 + push!(algs, (; kwargs...) -> MPDeC(k; kwargs...), + (; kwargs...) -> MPDeC(k; nodes = :lagrange, kwargs...)) + end for alg in algs for prob in probs sol1 = solve(prob_pds_linmod, alg(), dt = 0.1) @@ -2230,13 +2240,21 @@ end end end - #Here we check if the RK methods on which the MPRK schemes are based integrate + # Here we check if the RK methods on which the MPRK schemes are based integrate # u'(t) = q * t^(q-1) exactly for q from 1 to the order of the method. + # This is also true for MPDeC as long as the theta matrix is nonnegative, i.e. K = 2. + # Nevertheless, the results of most MPDeC schemes are good enough to pass this test @testset "Exact solutions (RK)" begin #TODO: Add MPDeC - Requires PDSProblem - algs = (MPE(), MPRK22(0.5), MPRK22(1.0), MPRK22(2.0), - MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), MPRK43II(0.5), - MPRK43II(2.0 / 3.0), SSPMPRK22(0.5, 1.0), SSPMPRK43()) + algs = [MPE(), MPRK22(0.5), MPRK22(1.0), MPRK22(2.0), + MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), MPRK43II(0.5), + MPRK43II(2.0 / 3.0), SSPMPRK22(0.5, 1.0), SSPMPRK43()] + for k in 2:10 + push!(algs, MPDeC(k)) + if k != 9 + push!(algs, MPDeC(k, nodes = :lagrange)) + end + end @testset "$alg, $q" for alg in algs, q in 1:PositiveIntegrators.alg_order(alg) f(t) = q * t^(q - 1) function prod!(P, u, p, t) From a6b9254c69b5b87888335d32f790a6f4e98a19a0 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Fri, 14 Mar 2025 14:57:46 +0100 Subject: [PATCH 13/38] fix testset 'zero intial values' --- src/mpdec.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/mpdec.jl b/src/mpdec.jl index e29dcece..b167c4c4 100644 --- a/src/mpdec.jl +++ b/src/mpdec.jl @@ -41,8 +41,9 @@ end function small_constant_function_MPDeC(type) if type == Float64 - #TODO: Adjust default value based on K - small_constant = floatmin(type) + # small_constant is chosen such that + # the testet "Zero initial values" passes. + small_constant = 1e-300 else small_constant = floatmin(type) end From 0373369d3b26284bc791c01c80b6338a48456787 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Fri, 14 Mar 2025 18:03:07 +0100 Subject: [PATCH 14/38] code restructuring --- src/mpdec.jl | 199 ++++++++++++++++++++++++++++++++++----------------- 1 file changed, 133 insertions(+), 66 deletions(-) diff --git a/src/mpdec.jl b/src/mpdec.jl index b167c4c4..cbddd618 100644 --- a/src/mpdec.jl +++ b/src/mpdec.jl @@ -1,18 +1,19 @@ +#TODO: Comment on choice of small_consant """ MPDeC(K; [nodes = :gausslobatto, linsolve = ..., small_constant = ...]) A family of arbitrary order modified Patankar-Runge-Kutta algorithms for production-destruction systems. Each member of this family is an adaptive, one-step method which is Kth order accurate, unconditionally positivity-preserving, and linearly -implicit. +implicit. The integer K must be chosen to satisfy 2 ≤ K ≤ 10. +Available node choices are Lagrange or Gauss-Lobatto nodes, with the latter being the default. +These methods support adaptive time stepping using the numerical solution obtained with one correction step less, as lower order approximation to estimate the error. -This method supports adaptive time stepping, using TODO: ???, as lower order approximations to estimate the error. - -These schemes were introduced by Torlo and Öffner (2020) for conservative production-destruction systems. +The MPDeC schemes were introduced by Torlo and Öffner (2020) for autonomous conservative production-destruction systems. For nonconservative production–destruction systems we use a straight forward extension analogous to [`MPE`](@ref). - -TODO: We need a sentences about the nodes choices. Literature for Gauss-Lobatto? +A general discussion of DeC schemes applied to non-autonomous differential equations +and using general integration nodes is given by Ong and Spiteri (2020). The MPDeC methods require the special structure of a [`PDSProblem`](@ref) or a [`ConservativePDSProblem`](@ref). @@ -29,7 +30,9 @@ to avoid divisions by zero. You can pass a value explicitly, otherwise `small_co - Davide Torlo, and Philipp Öffner. "Arbitrary high-order, conservative and positivity preserving Patankar-type deferred correction schemes." Applied Numerical Mathematics 153 (2020): 15-34. -- TODO: Literature non-autonomous DeC +- Benjamin W. Ong & Raymond J. Spiteri. + "Deferred Correction Methods for Ordinary Differential Equations." + Journal of Scientific Computing 83 (2020): Article 60 """ struct MPDeC{T, N, F, T2} <: OrdinaryDiffEqAdaptiveAlgorithm K::T @@ -321,88 +324,150 @@ function initialize!(integrator, cache::MPDeCConstantCache) end # out-of-place -function build_mpdec_matrix_and_rhs(uprev, m, f, C, p, t, dt, nodes, theta, small_constant) - N, M = size(C) - M = M - 1 - - Mmat = zeros(eltype(C), N, N) - rhs = similar(uprev) - rhs .= uprev - - # P stores the production matrix in in-place computations. - # In out-of-place computations it is not needed - P = nothing +function build_mpdec_matrix_and_rhs_oop(uprev, m, f, C, p, t, dt, nodes, theta, + small_constant) + N = length(uprev) if f isa PDSFunction # Additional destruction terms - build_mpdec_matrix_and_rhs!(Mmat, rhs, m, f.p, P, C, p, t, dt, nodes, theta, - small_constant, f.d) + Mmat, rhs = _build_mpdec_matrix_and_rhs_oop(uprev, m, f.p, C, p, t, dt, nodes, + theta, + small_constant, f.d) else # No additional destruction terms - build_mpdec_matrix_and_rhs!(Mmat, rhs, m, f.p, P, C, p, t, dt, nodes, theta, - small_constant) + Mmat, rhs = _build_mpdec_matrix_and_rhs_oop(uprev, m, f.p, C, p, t, dt, nodes, + theta, + small_constant) end - if C isa StaticArray - return SMatrix(Mmat), SVector(rhs) + if uprev isa StaticArray + return SMatrix{N, N}(Mmat), SVector{N}(rhs) else return Mmat, rhs end end -# in-place for dense arrays -function build_mpdec_matrix_and_rhs!(Mmat, rhs, m, prod, P, C, p, t, dt, nodes, theta, - small_constant, - dest = nothing, d = nothing) +# out-of-place for dense arrays +@muladd function _build_mpdec_matrix_and_rhs_oop(uprev, m, prod, C, p, t, dt, nodes, theta, + small_constant, dest = nothing) N, M = size(C) M = M - 1 + # Create linear system matrix and rhs + Mmat = zeros(eltype(C), N, N) + rhs = similar(uprev) + + # Initialize linear system matrix and rhs fill!(Mmat, zero(eltype(Mmat))) oneMmat = one(eltype(Mmat)) - for i in 1:N + @inbounds for i in 1:N Mmat[i, i] = oneMmat end + rhs .= uprev - #TODO: use predefined memory in inplace computations σ = add_small_constant(C[:, m], small_constant) - for r in 1:(M + 1) + @fastmath @inbounds @simd for r in 1:(M + 1) th = theta[r, m] dt_th = dt * th - #TODO: This should be checked earlier and only once - if isinplace(prod, 4) - prod(P, C[:, r], p, t + nodes[r] * dt) - if !isnothing(dest) - dest(d, C[:, r], p, t + nodes[r] * dt) - end + P = prod(C[:, r], p, t + nodes[r] * dt) + if !isnothing(dest) + d = dest(C[:, r], p, t + nodes[r] * dt) else - P = prod(C[:, r], p, t + nodes[r] * dt) - if !isnothing(dest) - d = dest(C[:, r], p, t + nodes[r] * dt) + d = nothing + end + _build_mpdec_matrix_and_rhs!(Mmat, rhs, P, dt_th, σ, d) + end + + return Mmat, rhs +end + +# in-place for dense arrays +@muladd function build_mpdec_matrix_and_rhs_ip!(Mmat, rhs, m, prod, P, C, p, t, dt, σ, + nodes, theta, small_constant, + dest = nothing, d = nothing) + N, M = size(C) + M = M - 1 + + fill!(Mmat, zero(eltype(Mmat))) + oneMmat = one(eltype(Mmat)) + @inbounds for i in 1:N + Mmat[i, i] = oneMmat + end + + σ .= C[:, m] .+ small_constant + + @fastmath @inbounds @simd for r in 1:(M + 1) + th = theta[r, m] + dt_th = dt * th + + prod(P, C[:, r], p, t + nodes[r] * dt) + if !isnothing(dest) + dest(d, C[:, r], p, t + nodes[r] * dt) + end + + _build_mpdec_matrix_and_rhs!(Mmat, rhs, P, dt_th, σ, d) + end +end + +#= +function _build_mpdec_matrix_and_rhs_old!(Mmat, rhs, P, dt_th, σ, d = nothing) + N = length(rhs) + @fastmath @inbounds @simd for i in 1:N + # Add nonconservative destruction terms to diagonal (PDSFunctions only!) + if !isnothing(d) + if dt_th >= 0 + Mmat[i, i] += dt_th * d[i] / σ[i] + rhs[i] += dt_th * P[i, i] + else + Mmat[i, i] -= dt_th * P[i, i] / σ[i] + rhs[i] -= dt_th * d[i] + end + end + @fastmath @inbounds @simd for j in 1:N + if j != i + if dt_th >= 0 + Mmat[i, j] -= dt_th * P[i, j] / σ[j] + Mmat[i, i] += dt_th * P[j, i] / σ[i] # P[j, i] = D[i, j] + else + Mmat[i, j] += dt_th * P[j, i] / σ[j] # P[j, i] = D[i, j] + Mmat[i, i] -= dt_th * P[i, j] / σ[i] + end + end + end + end +end +=# +function _build_mpdec_matrix_and_rhs!(Mmat, rhs, P, dt_th, σ, d = nothing) + @fastmath @inbounds @simd for I in CartesianIndices(P) + if !iszero(P[I]) + dt_th_P = dt_th * P[I] + if I[1] != I[2] + if dt_th ≥ 0 + Mmat[I] -= dt_th_P / σ[I[2]] + Mmat[I[2], I[2]] += dt_th_P / σ[I[2]] + else + Mmat[I[2], I[1]] += dt_th_P / σ[I[1]] + Mmat[I[1], I[1]] -= dt_th_P / σ[I[1]] + end + else # diagonal elements + if dt_th >= 0 + rhs[I[1]] += dt_th_P + else + Mmat[I] -= dt_th_P / σ[I[1]] + end end end + end - for i in 1:N - # Add nonconservative destruction terms to diagonal (PDSFunctions only!) - if !isnothing(dest) - if th >= 0 + if !isnothing(d) + @fastmath @inbounds @simd for i in eachindex(d) + if !iszero(d[i]) + if dt_th >= 0 Mmat[i, i] += dt_th * d[i] / σ[i] - rhs[i] += dt_th * P[i, i] else - Mmat[i, i] -= dt_th * P[i, i] / σ[i] rhs[i] -= dt_th * d[i] end end - for j in 1:N - if j != i - if th >= 0 - Mmat[i, j] -= dt_th * P[i, j] / σ[j] - Mmat[i, i] += dt_th * P[j, i] / σ[i] # P[j, i] = D[i, j] - else - Mmat[i, j] += dt_th * P[j, i] / σ[j] # P[j, i] = D[i, j] - Mmat[i, i] -= dt_th * P[i, j] / σ[i] - end - end - end end end end @@ -420,9 +485,9 @@ end for _ in 1:K C .= C2 for m in 2:(M + 1) - Mmat, rhs = build_mpdec_matrix_and_rhs(uprev, m, f, C, p, t, dt, nodes, theta, - small_constant) - + Mmat, rhs = build_mpdec_matrix_and_rhs_oop(uprev, m, f, C, p, t, dt, nodes, + theta, + small_constant) # solve linear system linprob = LinearProblem(Mmat, rhs) sol = solve(linprob, alg.linsolve) @@ -535,9 +600,10 @@ end C .= C2 for m in 2:(M + 1) linsolve_rhs .= uprev - build_mpdec_matrix_and_rhs!(P2, linsolve_rhs, m, f.p, P, C, p, t, dt, nodes, - theta, - small_constant, f.d, d) + build_mpdec_matrix_and_rhs_ip!(P2, linsolve_rhs, m, f.p, P, C, p, t, dt, σ, + nodes, + theta, + small_constant, f.d, d) # Same as linres = P2 \ linsolve_rhs linsolve.A = P2 @@ -578,9 +644,10 @@ end C .= C2 for m in 2:(M + 1) linsolve_rhs .= uprev - build_mpdec_matrix_and_rhs!(P2, linsolve_rhs, m, f.p, P, C, p, t, dt, nodes, - theta, - small_constant) + build_mpdec_matrix_and_rhs_ip!(P2, linsolve_rhs, m, f.p, P, C, p, t, dt, σ, + nodes, + theta, + small_constant) # Same as linres = P2 \ linsolve_rhs linsolve.A = P2 From e681d504dcb24501c367c10f1d3ad533a03cfe2d Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Mon, 17 Mar 2025 11:01:30 +0100 Subject: [PATCH 15/38] performance enhancement --- src/PositiveIntegrators.jl | 2 +- src/mpdec.jl | 23 +++++++++++++++++------ 2 files changed, 18 insertions(+), 7 deletions(-) diff --git a/src/PositiveIntegrators.jl b/src/PositiveIntegrators.jl index 3ac47a5d..737fac19 100644 --- a/src/PositiveIntegrators.jl +++ b/src/PositiveIntegrators.jl @@ -6,7 +6,7 @@ using Statistics: median using SparseArrays: SparseArrays, AbstractSparseMatrix, issparse, nonzeros, nzrange, rowvals, spdiagm -using StaticArrays: SVector, SMatrix, StaticArray, @SVector, @SMatrix +using StaticArrays: SVector, SMatrix, StaticArray, @SVector, @SMatrix, MMatrix using FastBroadcast: @.. using MuladdMacro: @muladd diff --git a/src/mpdec.jl b/src/mpdec.jl index cbddd618..a038c733 100644 --- a/src/mpdec.jl +++ b/src/mpdec.jl @@ -352,12 +352,15 @@ end N, M = size(C) M = M - 1 - # Create linear system matrix and rhs - Mmat = zeros(eltype(C), N, N) + # Create linear system matrix and rhs + if uprev isa StaticArray + Mmat = MMatrix{N, N}(zeros(eltype(uprev), N, N)) + else + Mmat = zeros(eltype(uprev), N, N) + end rhs = similar(uprev) - # Initialize linear system matrix and rhs - fill!(Mmat, zero(eltype(Mmat))) + # Initialize oneMmat = one(eltype(Mmat)) @inbounds for i in 1:N Mmat[i, i] = oneMmat @@ -476,8 +479,16 @@ end @unpack alg, t, dt, uprev, f, p = integrator @unpack K, M, nodes, theta, small_constant = cache - C = zeros(length(uprev), M + 1) - C2 = zeros(length(uprev), M + 1) + N = length(uprev) + + if uprev isa StaticArray + C = MMatrix{N, M + 1}(zeros(N, M + 1)) + C2 = MMatrix{N, M + 1}(zeros(N, M + 1)) + else + C = zeros(N, M + 1) + C2 = zeros(N, M + 1) + end + for i in 1:(M + 1) C2[:, i] = uprev end From f859596326ff25f36609119cbd313caf64449963 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Mon, 17 Mar 2025 14:17:28 +0100 Subject: [PATCH 16/38] check theta < 0 only once --- src/mpdec.jl | 67 ++++++++++++++++++++++++++++++++++++++++++++++------ 1 file changed, 60 insertions(+), 7 deletions(-) diff --git a/src/mpdec.jl b/src/mpdec.jl index a038c733..8fe524b9 100644 --- a/src/mpdec.jl +++ b/src/mpdec.jl @@ -440,23 +440,76 @@ function _build_mpdec_matrix_and_rhs_old!(Mmat, rhs, P, dt_th, σ, d = nothing) end end =# -function _build_mpdec_matrix_and_rhs!(Mmat, rhs, P, dt_th, σ, d = nothing) +function _build_mpdec_matrix_and_rhs!(M, rhs, P, dt_th, σ, d = nothing) + Base.require_one_based_indexing(M, P, σ) + @assert size(M, 1) == size(M, 2) == size(P, 1) == size(P, 2) == length(σ) + + if dt_th ≥ 0 + @fastmath @inbounds @simd for I in CartesianIndices(P) + if !iszero(P[I]) + dt_th_P = dt_th * P[I] + if I[1] != I[2] + M[I] -= dt_th_P / σ[I[2]] + M[I[2], I[2]] += dt_th_P / σ[I[2]] + else # diagonal elements + rhs[I[1]] += dt_th_P + end + end + end + + if !isnothing(d) + @fastmath @inbounds @simd for i in eachindex(d) + if !iszero(d[i]) + M[i, i] += dt_th * d[i] / σ[i] + end + end + end + else # dt_th ≤ 0 + @fastmath @inbounds @simd for I in CartesianIndices(P) + if !iszero(P[I]) + dt_th_P = dt_th * P[I] + if I[1] != I[2] + M[I[2], I[1]] += dt_th_P / σ[I[1]] + M[I[1], I[1]] -= dt_th_P / σ[I[1]] + else # diagonal elements + M[I] -= dt_th_P / σ[I[1]] + end + end + end + + if !isnothing(d) + @fastmath @inbounds @simd for i in eachindex(d) + if !iszero(d[i]) + rhs[i] -= dt_th * d[i] + end + end + end + end +end + +# optimized version for Tridiagonal matrices +function _build_mpdec_matrix_and_rhs!(M::Tridiagonal, rhs, P::Tridiagonal, dt_th, σ, + d = nothing) + Base.require_one_based_indexing(M.dl, M.d, M.du, P.dl, P.d, P.du, σ) + @assert length(M.dl) + 1 == length(M.d) == length(M.du) + 1 == + length(P.dl) + 1 == length(P.d) == length(P.du) + 1 == length(σ) + @fastmath @inbounds @simd for I in CartesianIndices(P) if !iszero(P[I]) dt_th_P = dt_th * P[I] if I[1] != I[2] if dt_th ≥ 0 - Mmat[I] -= dt_th_P / σ[I[2]] - Mmat[I[2], I[2]] += dt_th_P / σ[I[2]] + M[I] -= dt_th_P / σ[I[2]] + M[I[2], I[2]] += dt_th_P / σ[I[2]] else - Mmat[I[2], I[1]] += dt_th_P / σ[I[1]] - Mmat[I[1], I[1]] -= dt_th_P / σ[I[1]] + M[I[2], I[1]] += dt_th_P / σ[I[1]] + M[I[1], I[1]] -= dt_th_P / σ[I[1]] end else # diagonal elements if dt_th >= 0 rhs[I[1]] += dt_th_P else - Mmat[I] -= dt_th_P / σ[I[1]] + M[I] -= dt_th_P / σ[I[1]] end end end @@ -466,7 +519,7 @@ function _build_mpdec_matrix_and_rhs!(Mmat, rhs, P, dt_th, σ, d = nothing) @fastmath @inbounds @simd for i in eachindex(d) if !iszero(d[i]) if dt_th >= 0 - Mmat[i, i] += dt_th * d[i] / σ[i] + M[i, i] += dt_th * d[i] / σ[i] else rhs[i] -= dt_th * d[i] end From 68c7eca1f7bc5da802841156c2d94e1f5f8b6a9d Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Mon, 17 Mar 2025 18:33:05 +0100 Subject: [PATCH 17/38] Support Tridiagonal --- src/mpdec.jl | 86 ++++++++++++++++++++++++++++++++++++---------------- 1 file changed, 60 insertions(+), 26 deletions(-) diff --git a/src/mpdec.jl b/src/mpdec.jl index 8fe524b9..f0401db5 100644 --- a/src/mpdec.jl +++ b/src/mpdec.jl @@ -391,10 +391,18 @@ end N, M = size(C) M = M - 1 - fill!(Mmat, zero(eltype(Mmat))) oneMmat = one(eltype(Mmat)) - @inbounds for i in 1:N - Mmat[i, i] = oneMmat + zeroMmat = zero(eltype(Mmat)) + + if Mmat isa Tridiagonal + Mmat.d .= oneMmat + Mmat.du .= zeroMmat + Mmat.dl .= zeroMmat + else + fill!(Mmat, zeroMmat) + @inbounds for i in 1:N + Mmat[i, i] = oneMmat + end end σ .= C[:, m] .+ small_constant @@ -493,36 +501,62 @@ function _build_mpdec_matrix_and_rhs!(M::Tridiagonal, rhs, P::Tridiagonal, dt_th Base.require_one_based_indexing(M.dl, M.d, M.du, P.dl, P.d, P.du, σ) @assert length(M.dl) + 1 == length(M.d) == length(M.du) + 1 == length(P.dl) + 1 == length(P.d) == length(P.du) + 1 == length(σ) + + if dt_th ≥ 0 + @fastmath @inbounds @simd for i in eachindex(P.d, rhs) + rhs[i] += dt_th * P.d[i] + end - @fastmath @inbounds @simd for I in CartesianIndices(P) - if !iszero(P[I]) - dt_th_P = dt_th * P[I] - if I[1] != I[2] - if dt_th ≥ 0 - M[I] -= dt_th_P / σ[I[2]] - M[I[2], I[2]] += dt_th_P / σ[I[2]] - else + for i in eachindex(M.dl, P.dl) + dt_th_P = dt_th * P.dl[i] + M.dl[i] -= dt_th_P / σ[i] + M.d[i] += dt_th_P / σ[i] + end + + for i in eachindex(M.du, P.du) + dt_th_P = dt_th * P.du[i] + M.du[i] -= dt_th_P / σ[i + 1] + M.d[i + 1] += dt_th_P / σ[i + 1] + end + + if !isnothing(d) + @fastmath @inbounds @simd for i in eachindex(M.d, σ, d) + M.d[i] += dt_th * d[i] / σ[i] + end + end + else # dt_th ≤ 0 + #= + @fastmath @inbounds @simd for I in CartesianIndices(P) + if !iszero(P[I]) + dt_th_P = dt_th * P[I] + if I[1] != I[2] M[I[2], I[1]] += dt_th_P / σ[I[1]] M[I[1], I[1]] -= dt_th_P / σ[I[1]] end - else # diagonal elements - if dt_th >= 0 - rhs[I[1]] += dt_th_P - else - M[I] -= dt_th_P / σ[I[1]] - end end end - end + =# - if !isnothing(d) - @fastmath @inbounds @simd for i in eachindex(d) - if !iszero(d[i]) - if dt_th >= 0 - M[i, i] += dt_th * d[i] / σ[i] - else - rhs[i] -= dt_th * d[i] - end + @fastmath @inbounds @simd for i in eachindex(M.d, P.d, σ) + M.d[i] -= dt_th * P.d[i] / σ[i] + end + + for i in eachindex(M.dl, P.dl) + dt_th_P = dt_th * P.dl[i] + M.du[i] += dt_th_P / σ[i+1] + M.d[i+1] -= dt_th_P / σ[i+1] + end + + for i in eachindex(M.du, P.du) + dt_th_P = dt_th * P.du[i] + M.dl[i] += dt_th_P / σ[i] + M.d[i] -= dt_th_P / σ[i] + end + + + if !isnothing(d) + @fastmath @inbounds @simd for i in eachindex(rhs, d) + rhs[i] -= dt_th * d[i] end end end From d857f1976eee6bf12d2ebef6aa6d8a18233a145b Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Tue, 18 Mar 2025 16:05:32 +0100 Subject: [PATCH 18/38] Support for sparse matrices --- src/mpdec.jl | 132 +++++++++++++++++++++++++++++++++++++++-------- test/runtests.jl | 86 +++++++++++++++++++----------- 2 files changed, 165 insertions(+), 53 deletions(-) diff --git a/src/mpdec.jl b/src/mpdec.jl index f0401db5..bd2fd91e 100644 --- a/src/mpdec.jl +++ b/src/mpdec.jl @@ -385,7 +385,7 @@ end end # in-place for dense arrays -@muladd function build_mpdec_matrix_and_rhs_ip!(Mmat, rhs, m, prod, P, C, p, t, dt, σ, +@muladd function build_mpdec_matrix_and_rhs_ip!(Mmat, rhs, m, prod, P, C, p, t, dt, σ, tmp, nodes, theta, small_constant, dest = nothing, d = nothing) N, M = size(C) @@ -398,6 +398,21 @@ end Mmat.d .= oneMmat Mmat.du .= zeroMmat Mmat.dl .= zeroMmat + elseif issparse(Mmat) + # Fill sparse matrix with zeros without changing the sparsity pattern, see + # https://github.com/JuliaSparse/SparseArrays.jl/issues/190#issuecomment-1186690008. + fill!(Mmat.nzval, false) + + M_rows = rowvals(Mmat) + M_vals = nonzeros(Mmat) + for j in 1:N + for idx_M in nzrange(Mmat, j) + i = M_rows[idx_M] + if i == j + M_vals[idx_M] = oneMmat + end + end + end else fill!(Mmat, zeroMmat) @inbounds for i in 1:N @@ -416,7 +431,11 @@ end dest(d, C[:, r], p, t + nodes[r] * dt) end - _build_mpdec_matrix_and_rhs!(Mmat, rhs, P, dt_th, σ, d) + if issparse(Mmat) + _build_mpdec_matrix_and_rhs!(Mmat, rhs, P, dt_th, σ, d, tmp) + else + _build_mpdec_matrix_and_rhs!(Mmat, rhs, P, dt_th, σ, d) + end end end @@ -501,7 +520,7 @@ function _build_mpdec_matrix_and_rhs!(M::Tridiagonal, rhs, P::Tridiagonal, dt_th Base.require_one_based_indexing(M.dl, M.d, M.du, P.dl, P.d, P.du, σ) @assert length(M.dl) + 1 == length(M.d) == length(M.du) + 1 == length(P.dl) + 1 == length(P.d) == length(P.du) + 1 == length(σ) - + if dt_th ≥ 0 @fastmath @inbounds @simd for i in eachindex(P.d, rhs) rhs[i] += dt_th * P.d[i] @@ -512,7 +531,7 @@ function _build_mpdec_matrix_and_rhs!(M::Tridiagonal, rhs, P::Tridiagonal, dt_th M.dl[i] -= dt_th_P / σ[i] M.d[i] += dt_th_P / σ[i] end - + for i in eachindex(M.du, P.du) dt_th_P = dt_th * P.du[i] M.du[i] -= dt_th_P / σ[i + 1] @@ -525,35 +544,22 @@ function _build_mpdec_matrix_and_rhs!(M::Tridiagonal, rhs, P::Tridiagonal, dt_th end end else # dt_th ≤ 0 - #= - @fastmath @inbounds @simd for I in CartesianIndices(P) - if !iszero(P[I]) - dt_th_P = dt_th * P[I] - if I[1] != I[2] - M[I[2], I[1]] += dt_th_P / σ[I[1]] - M[I[1], I[1]] -= dt_th_P / σ[I[1]] - end - end - end - =# - @fastmath @inbounds @simd for i in eachindex(M.d, P.d, σ) M.d[i] -= dt_th * P.d[i] / σ[i] end for i in eachindex(M.dl, P.dl) dt_th_P = dt_th * P.dl[i] - M.du[i] += dt_th_P / σ[i+1] - M.d[i+1] -= dt_th_P / σ[i+1] + M.du[i] += dt_th_P / σ[i + 1] + M.d[i + 1] -= dt_th_P / σ[i + 1] end - + for i in eachindex(M.du, P.du) dt_th_P = dt_th * P.du[i] M.dl[i] += dt_th_P / σ[i] M.d[i] -= dt_th_P / σ[i] end - if !isnothing(d) @fastmath @inbounds @simd for i in eachindex(rhs, d) rhs[i] -= dt_th * d[i] @@ -562,6 +568,88 @@ function _build_mpdec_matrix_and_rhs!(M::Tridiagonal, rhs, P::Tridiagonal, dt_th end end +# optimized version for sparse matrices +function _build_mpdec_matrix_and_rhs!(M::AbstractSparseMatrix, rhs, P::AbstractSparseMatrix, + dt_th, σ, + d = nothing, tmp = nothing) + Base.require_one_based_indexing(M, P, σ) + @assert size(M, 1) == size(M, 2) == size(P, 1) == size(P, 2) == length(σ) + if !isnothing(d) + Base.require_one_based_indexing(d) + @assert length(σ) == length(d) + end + + # By construction M and P share the same sparsity pattern. + M_rows = rowvals(M) + M_vals = nonzeros(M) + P_rows = rowvals(P) + P_vals = nonzeros(P) + n = size(M, 2) + + # tmp[j] = M[j,j] + zeroP = zero(eltype(P)) + for i in eachindex(tmp) + tmp[i] = zeroP + end + + if dt_th ≥ 0 + for j in 1:n # j is column index + for idx_P in nzrange(P, j) + dt_th_P = dt_th * P_vals[idx_P] + i = P_rows[idx_P] # i is row index + if i != j + M_vals[idx_P] -= dt_th_P / σ[j] + tmp[j] += dt_th_P / σ[j] + else + rhs[i] += dt_th_P + end + end + end + + if !isnothing(d) + for i in eachindex(d) + tmp[i] += dt_th * d[i] / σ[i] + end + end + + for j in 1:n + for idx_P in nzrange(P, j) + i = P_rows[idx_P] + if i == j + M_vals[idx_P] += tmp[j] + end + end + end + else # dt ≤ 0 + for j in 1:n # j is column index + for idx_P in nzrange(P, j) + dt_th_P = dt_th * P_vals[idx_P] + i = P_rows[idx_P] # i is row index + if i != j + #TODO: In general (j,i) is not in the sparsity pattern! + M[j, i] += dt_th_P / σ[i] + tmp[i] -= dt_th_P / σ[i] + else + M_vals[idx_P] -= dt_th_P / σ[i] + end + end + end + + for j in 1:n + for idx_P in nzrange(P, j) + i = P_rows[idx_P] + if i == j + M_vals[idx_P] += tmp[j] + end + end + end + + if !isnothing(d) + @.. broadcast=false rhs-=dt_th * d + end + end +end + @muladd function perform_step!(integrator, cache::MPDeCConstantCache, repeat_step = false) @unpack alg, t, dt, uprev, f, p = integrator @unpack K, M, nodes, theta, small_constant = cache @@ -698,7 +786,7 @@ end C .= C2 for m in 2:(M + 1) linsolve_rhs .= uprev - build_mpdec_matrix_and_rhs_ip!(P2, linsolve_rhs, m, f.p, P, C, p, t, dt, σ, + build_mpdec_matrix_and_rhs_ip!(P2, linsolve_rhs, m, f.p, P, C, p, t, dt, σ, tmp, nodes, theta, small_constant, f.d, d) @@ -742,7 +830,7 @@ end C .= C2 for m in 2:(M + 1) linsolve_rhs .= uprev - build_mpdec_matrix_and_rhs_ip!(P2, linsolve_rhs, m, f.p, P, C, p, t, dt, σ, + build_mpdec_matrix_and_rhs_ip!(P2, linsolve_rhs, m, f.p, P, C, p, t, dt, σ, tmp, nodes, theta, small_constant) diff --git a/test/runtests.jl b/test/runtests.jl index 286eff75..8f084e9c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1300,7 +1300,6 @@ end # Here we check that in-place and out-of-place implementations # deliver the same results @testset "Different matrix types (conservative)" begin - #TODO: Add MPDeC - REQUIRES spare and tridiagonal matrices prod_1! = (P, u, p, t) -> begin fill!(P, zero(eltype(P))) for i in 1:(length(u) - 1) @@ -1336,11 +1335,21 @@ end tspan = (0.0, 1.0) dt = 0.25 - @testset "$alg" for alg in (MPE(), - MPRK22(0.5), MPRK22(1.0), - MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), - MPRK43II(2.0 / 3.0), MPRK43II(0.5), - SSPMPRK22(0.5, 1.0), SSPMPRK43()) + algs = [ + MPE(), + MPRK22(0.5), + MPRK22(1.0), + MPRK43I(1.0, 0.5), + MPRK43I(0.5, 0.75), + MPRK43II(2.0 / 3.0), + MPRK43II(0.5), + SSPMPRK22(0.5, 1.0), + SSPMPRK43() + ] + for k in 2:10 + push!(algs, MPDeC(k), MPDeC(k; nodes = :lagrange)) + end + @testset "$alg" for alg in algs for prod! in (prod_1!, prod_2!, prod_3!) prod = (u, p, t) -> begin P = similar(u, (length(u), length(u))) @@ -1386,7 +1395,6 @@ end @testset "Different matrix types (conservative, adaptive)" begin #TODO: MPE, SSPMPRK43 are not adaptive and need to be removed (solve without giving dt!) - #TODO: Add MPDeC prod_1! = (P, u, p, t) -> begin fill!(P, zero(eltype(P))) for i in 1:(length(u) - 1) @@ -1424,11 +1432,15 @@ end rtol = sqrt(eps(Float32)) - @testset "$alg" for alg in (MPE(), - MPRK22(0.5), MPRK22(1.0), - MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), - MPRK43II(2.0 / 3.0), MPRK43II(0.5), - SSPMPRK22(0.5, 1.0), SSPMPRK43()) + algs = [MPE(), + MPRK22(0.5), MPRK22(1.0), + MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), + MPRK43II(2.0 / 3.0), MPRK43II(0.5), + SSPMPRK22(0.5, 1.0), SSPMPRK43()] + for k in 2:10 + push!(algs, MPDeC(k), MPDeC(k; nodes = :lagrange)) + end + @testset "$alg" for alg in algs for prod! in (prod_1!, prod_2!, prod_3!) prod = (u, p, t) -> begin P = similar(u, (length(u), length(u))) @@ -1473,7 +1485,6 @@ end # Here we check that in-place and out-of-place implementations # deliver the same results @testset "Different matrix types (nonconservative)" begin - #TODO: Add MPDeC prod_1! = (P, u, p, t) -> begin fill!(P, zero(eltype(P))) for i in 1:(length(u) - 1) @@ -1540,11 +1551,16 @@ end tspan = (0.0, 1.0) dt = 0.25 - @testset "$alg" for alg in (MPE(), - MPRK22(0.5), MPRK22(1.0), - MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), - MPRK43II(2.0 / 3.0), MPRK43II(0.5), - SSPMPRK22(0.5, 1.0), SSPMPRK43()) + algs = [MPE(), + MPRK22(0.5), MPRK22(1.0), + MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), + MPRK43II(2.0 / 3.0), MPRK43II(0.5), + SSPMPRK22(0.5, 1.0), SSPMPRK43()] + for k in 2:10 + push!(algs, MPDeC(k), MPDeC(k; nodes = :lagrange)) + end + + @testset "$alg" for alg in algs for (prod!, dest!) in zip((prod_1!, prod_2!, prod_3!), (dest_1!, dest_2!, dest_3!)) prod = (u, p, t) -> begin @@ -1600,7 +1616,6 @@ end @testset "Different matrix types (nonconservative, adaptive)" begin #TODO: MPE, SSPMPRK43 are not adaptive and need to be removed (solve without giving dt) - #TODO: Add MPDeC - Requires PDSProblem and adaptivity prod_1! = (P, u, p, t) -> begin fill!(P, zero(eltype(P))) for i in 1:(length(u) - 1) @@ -1667,12 +1682,17 @@ end tspan = (0.0, 1.0) dt = 0.25 + algs = [MPE(), + MPRK22(0.5), MPRK22(1.0), + MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), + MPRK43II(2.0 / 3.0), MPRK43II(0.5), + SSPMPRK22(0.5, 1.0), SSPMPRK43()] + for k in 2:10 + push!(algs, MPDeC(k), MPDeC(k; nodes = :lagrange)) + end + rtol = sqrt(eps(Float32)) - @testset "$alg" for alg in (MPE(), - MPRK22(0.5), MPRK22(1.0), - MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), - MPRK43II(2.0 / 3.0), MPRK43II(0.5), - SSPMPRK22(0.5, 1.0), SSPMPRK43()) + @testset "$alg" for alg in algs for (prod!, dest!) in zip((prod_1!, prod_2!, prod_3!), (dest_1!, dest_2!, dest_3!)) prod! = prod_3! @@ -1732,7 +1752,6 @@ end # defines the types of the Ps inside the algorithm caches. # We test sparse, tridiagonal, and dense matrices. @testset "Prototype type check" begin - #TODO: Add MPDeC #prod and dest functions prod_inner! = (P, u, p, t) -> begin fill!(P, zero(eltype(P))) @@ -1785,11 +1804,17 @@ end p_prototype = P_dense) prob_sparse2 = PDSProblem(prod_sparse!, dest!, u0, tspan; p_prototype = P_sparse) + + algs = [MPE(), + MPRK22(0.5), MPRK22(1.0), + MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), + MPRK43II(2.0 / 3.0), MPRK43II(0.5), + SSPMPRK22(0.5, 1.0), SSPMPRK43()] + for k in 2:10 + push!(algs, MPDeC(k), MPDeC(k; nodes = :lagrange)) + end #solve and test - for alg in (MPE(), MPRK22(0.5), MPRK22(1.0), MPRK43I(1.0, 0.5), - MPRK43I(0.5, 0.75), - MPRK43II(2.0 / 3.0), MPRK43II(0.5), SSPMPRK22(0.5, 1.0), - SSPMPRK43()) + for alg in algs for prob in (prob_default, prob_tridiagonal, prob_dense, prob_sparse, prob_default2, prob_tridiagonal2, prob_dense2, prob_sparse2) @@ -2065,7 +2090,7 @@ end algs = [MPE(), MPRK22(0.5), MPRK22(1.0), MPRK22(2.0), MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), MPRK43II(0.5), MPRK43II(2.0 / 3.0), SSPMPRK22(0.5, 1.0), SSPMPRK43()] - for k in 2:9 + for k in 2:10 push!(algs, MPDeC(k), MPDeC(k; nodes = :lagrange)) end @@ -2245,7 +2270,6 @@ end # This is also true for MPDeC as long as the theta matrix is nonnegative, i.e. K = 2. # Nevertheless, the results of most MPDeC schemes are good enough to pass this test @testset "Exact solutions (RK)" begin - #TODO: Add MPDeC - Requires PDSProblem algs = [MPE(), MPRK22(0.5), MPRK22(1.0), MPRK22(2.0), MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), MPRK43II(0.5), MPRK43II(2.0 / 3.0), SSPMPRK22(0.5, 1.0), SSPMPRK43()] From 1905e22bde1054b6232629a54f971dce2382975e Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Tue, 18 Mar 2025 23:04:04 +0100 Subject: [PATCH 19/38] fixed sparse matrix support --- src/mpdec.jl | 64 ++++++++++++++++++++++++++++++++---------------- test/runtests.jl | 8 +++--- 2 files changed, 47 insertions(+), 25 deletions(-) diff --git a/src/mpdec.jl b/src/mpdec.jl index bd2fd91e..77af6a68 100644 --- a/src/mpdec.jl +++ b/src/mpdec.jl @@ -394,6 +394,7 @@ end oneMmat = one(eltype(Mmat)) zeroMmat = zero(eltype(Mmat)) + #Initialize Mmat as identity matrix if Mmat isa Tridiagonal Mmat.d .= oneMmat Mmat.du .= zeroMmat @@ -579,6 +580,9 @@ function _build_mpdec_matrix_and_rhs!(M::AbstractSparseMatrix, rhs, P::AbstractS @assert length(σ) == length(d) end + #TODO: It is crucial that evaluation of the production function + # does not alter the sparsity pattern of p_prototype. This should be checked earlier. + # By construction M and P share the same sparsity pattern. M_rows = rowvals(M) M_vals = nonzeros(M) @@ -587,59 +591,72 @@ function _build_mpdec_matrix_and_rhs!(M::AbstractSparseMatrix, rhs, P::AbstractS n = size(M, 2) # tmp[j] = M[j,j] - zeroP = zero(eltype(P)) - for i in eachindex(tmp) - tmp[i] = zeroP - end + fill!(tmp, zero(eltype(tmp))) if dt_th ≥ 0 - for j in 1:n # j is column index - for idx_P in nzrange(P, j) + for j in 1:n # run through columns of P + for idx_P in nzrange(P, j) # run through rows of P + i = P_rows[idx_P] dt_th_P = dt_th * P_vals[idx_P] - i = P_rows[idx_P] # i is row index if i != j - M_vals[idx_P] -= dt_th_P / σ[j] - tmp[j] += dt_th_P / σ[j] + for idx_M in nzrange(M, j) + if M_rows[idx_M] == i + M_vals[idx_M] -= dt_th_P / σ[j] # M_ij <- P_ij + #break + end + end + tmp[j] += dt_th_P / σ[j] # M_jj <- P_ij = D_ji else - rhs[i] += dt_th_P + rhs[i] += dt_th_P # rhs_i <- P_ii end end end if !isnothing(d) for i in eachindex(d) - tmp[i] += dt_th * d[i] / σ[i] + tmp[i] += dt_th * d[i] / σ[i] # M_ii <- D_i end end for j in 1:n - for idx_P in nzrange(P, j) - i = P_rows[idx_P] + for idx_M in nzrange(M, j) + i = M_rows[idx_M] if i == j - M_vals[idx_P] += tmp[j] + M_vals[idx_M] += tmp[j] + #break end end end else # dt ≤ 0 for j in 1:n # j is column index for idx_P in nzrange(P, j) - dt_th_P = dt_th * P_vals[idx_P] i = P_rows[idx_P] # i is row index + dt_th_P = dt_th * P_vals[idx_P] if i != j - #TODO: In general (j,i) is not in the sparsity pattern! - M[j, i] += dt_th_P / σ[i] + for idx_M in nzrange(M, i) + if M_rows[idx_M] == j + M_vals[idx_M] += dt_th_P / σ[i] # M_ji <- P_ij + end + #break + end tmp[i] -= dt_th_P / σ[i] else - M_vals[idx_P] -= dt_th_P / σ[i] + for idx_M in nzrange(M, j) + if i == M_rows[idx_M] + M_vals[idx_M] -= dt_th_P / σ[i] # M_ij <- P_ij + #break + end + end end end end for j in 1:n - for idx_P in nzrange(P, j) - i = P_rows[idx_P] + for idx_M in nzrange(M, j) + i = M_rows[idx_M] if i == j - M_vals[idx_P] += tmp[j] + M_vals[idx_M] += tmp[j] + #break end end end @@ -736,6 +753,11 @@ function alg_cache(alg::MPDeC, u, rate_prototype, ::Type{uEltypeNoUnits}, tmp = zero(u) P = p_prototype(u, f) # stores evaluation of the production matrix P2 = p_prototype(u, f) # stores the linear system matrix + if issparse(P2) && alg.K > 2 + # Negative weights of MPDeC(K) , K >=3 require + # a symmetric sparsity pattern + P2 = P2 + P2' + end d = zero(u) σ = zero(u) C = zeros(eltype(u), length(u), alg.M + 1) diff --git a/test/runtests.jl b/test/runtests.jl index 8f084e9c..76ad2bc0 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1274,7 +1274,7 @@ end (; kwargs...) -> MPRK43II(2.0 / 3.0; kwargs...), (; kwargs...) -> SSPMPRK22(0.5, 1.0; kwargs...), (; kwargs...) -> SSPMPRK43(; kwargs...)] - for k in 2:9 + for k in 2:10 push!(algs, (; kwargs...) -> MPDeC(k; kwargs...), (; kwargs...) -> MPDeC(k; nodes = :lagrange, kwargs...)) end @@ -2141,7 +2141,7 @@ end algs = [MPE(), MPRK22(0.5), MPRK22(1.0), MPRK22(2.0), MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), MPRK43II(0.5), MPRK43II(2.0 / 3.0), SSPMPRK22(0.5, 1.0), SSPMPRK43()] - for k in 2:9 + for k in 2:10 push!(algs, MPDeC(k), MPDeC(k; nodes = :lagrange)) end @@ -2164,7 +2164,7 @@ end @testset "PDS problem library (adaptive schemes)" begin algs = [MPRK22(0.5), MPRK22(1.0), MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), MPRK43II(2.0 / 3.0), MPRK43II(0.5), SSPMPRK22(0.5, 1.0)] - for k in 2:9 + for k in 2:10 push!(algs, MPDeC(k), MPDeC(k; nodes = :lagrange)) end probs = (prob_pds_linmod, prob_pds_linmod_inplace, prob_pds_nonlinmod, @@ -2249,7 +2249,7 @@ end (; kwargs...) -> MPRK43II(0.5; kwargs...), (; kwargs...) -> SSPMPRK22(0.5, 1.0; kwargs...), (; kwargs...) -> SSPMPRK43(; kwargs...)] - for k in 2:9 + for k in 2:10 push!(algs, (; kwargs...) -> MPDeC(k; kwargs...), (; kwargs...) -> MPDeC(k; nodes = :lagrange, kwargs...)) end From 8b6cbdd98df7a433c27cb05f661311ca5f5c60d1 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Wed, 19 Mar 2025 12:22:12 +0100 Subject: [PATCH 20/38] minor changes and tests --- src/mpdec.jl | 256 +++++------------------------------------------ test/runtests.jl | 13 +++ 2 files changed, 39 insertions(+), 230 deletions(-) diff --git a/src/mpdec.jl b/src/mpdec.jl index 77af6a68..d0d5baf2 100644 --- a/src/mpdec.jl +++ b/src/mpdec.jl @@ -1,4 +1,3 @@ -#TODO: Comment on choice of small_consant """ MPDeC(K; [nodes = :gausslobatto, linsolve = ..., small_constant = ...]) @@ -8,8 +7,9 @@ Kth order accurate, unconditionally positivity-preserving, and linearly implicit. The integer K must be chosen to satisfy 2 ≤ K ≤ 10. Available node choices are Lagrange or Gauss-Lobatto nodes, with the latter being the default. These methods support adaptive time stepping using the numerical solution obtained with one correction step less, as lower order approximation to estimate the error. +The MPDeC schemes were introduced by Torlo and Öffner (2020) for autonomous conservative production-destruction systems and +further investigated in Torlo, Öffner and Ranocha (2022). -The MPDeC schemes were introduced by Torlo and Öffner (2020) for autonomous conservative production-destruction systems. For nonconservative production–destruction systems we use a straight forward extension analogous to [`MPE`](@ref). A general discussion of DeC schemes applied to non-autonomous differential equations @@ -23,13 +23,18 @@ algorithm from [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl) as keyword argument `linsolve`. You can also choose the parameter `small_constant` which is added to all Patankar-weight denominators to avoid divisions by zero. You can pass a value explicitly, otherwise `small_constant` is set to -`floatmin` of the floating point type used. +1e-300 in double precision computations or `floatmin` of the floating point type used. ## References - Davide Torlo, and Philipp Öffner. "Arbitrary high-order, conservative and positivity preserving Patankar-type deferred correction schemes." Applied Numerical Mathematics 153 (2020): 15-34. +- Davide Torlo, Philipp Öffner, and Hendrik Ranocha. + "Issues with positivity-preserving Patankar-type schemes." + Applied Numerical Mathematics 182 (2022): 117-147. + [DOI: 10.1016/j.apnum.2022.07.014](https://doi.org/10.1016/j.apnum.2022.07.014) + - Benjamin W. Ong & Raymond J. Spiteri. "Deferred Correction Methods for Ordinary Differential Equations." Journal of Scientific Computing 83 (2020): Article 60 @@ -85,9 +90,10 @@ isfsal(::MPDeC) = false function get_constant_parameters(alg::MPDeC) if alg.nodes == :lagrange - nodes = collect(0.0:(1 / alg.M):1.0) + nodes = collect(0:(1 / alg.M):1) + # M is one less than the methods order if alg.M == 1 - theta = [0.0 0.5; 0.0 0.5] + theta = [0 0.5; 0 0.5] elseif alg.M == 2 theta = [0.0 0.20833333333333337 0.16666666666666652; 0.0 0.33333333333333337 0.6666666666666667; @@ -128,7 +134,6 @@ function get_constant_parameters(alg::MPDeC) 0.0 -0.013405848450491272 -0.009863945578231281 -0.011894132653058165 -0.009674981103547253 -0.014615221088391195 0.04183673469395899 0.2070023148147584; 0.0 0.001623913454270594 0.0012093726379440242 0.0014349489795916492 0.0012093726379429626 0.0016239134542663791 -1.0658141036401503e-14 0.04346064814814099] elseif alg.M == 8 - #= theta = [0.0 0.03685850005511468 0.035688932980599594 0.03594029017857134 0.0358289241622568 0.035914937444895045 0.035803571428594694 0.03605492862659787 0.034885361552085214; 0.0 0.15387641920194012 0.2012610229276906 0.19782924107142996 0.1990828924162642 0.19819740685622378 0.1992857142857929 0.1969121334875581 0.2076895943564523; 0.0 -0.15861283344356245 -0.046840828924162636 0.009592633928562577 0.0021516754850381403 0.0065018050044045594 0.0016071428572104196 0.011744309413188603 -0.03273368606460281; @@ -138,16 +143,7 @@ function get_constant_parameters(alg::MPDeC) 0.0 -0.044477995480599525 -0.03434082892416224 -0.03923549107142321 -0.03488536155200972 -0.04232631999563363 0.01410714285702852 0.1258791473759011 -0.03273368606702931; 0.0 0.010777460868606693 0.008403880070546516 0.009492187499996696 0.008606701940021111 0.009860353284821599 0.006428571428614305 0.053813175154459714 0.2076895943564523; 0.0 -0.0011695670745149895 -0.0009182098765432678 -0.0010295758928574872 -0.0009435626102302086 -0.0010549286265453262 -0.0008035714285732354 -0.0019731385031036552 0.03488536155197153] - =# - theta = [0 1070017/29030400 32377/907200 12881/358400 4063/113400 41705/1161216 401/11200 149527/4147200 989/28350; - 0 2233547/14515200 22823/113400 35451/179200 2822/14175 115075/580608 279/1400 408317/2073600 2944/14175; - 0 -2302297/14515200 -21247/453600 1719/179200 61/28350 3775/580608 9/5600 24353/2073600 -464/14175; - 0 2797679/14515200 15011/113400 39967/179200 4094/14175 159175/580608 403/1400 542969/2073600 5248/14175; - 0 -31457/181440 -2903/22680 -351/2240 -227/2835 -125/36288 -9/280 343/25920 -454/2835; - 0 1573169/14515200 9341/113400 17217/179200 1154/14175 85465/580608 333/1400 368039/2073600 5248/14175; - 0 -645607/14515200 -15577/453600 -7031/179200 -989/28350 -24575/580608 79/5600 261023/2073600 -464/14175; - 0 156437/14515200 953/113400 243/25600 122/14175 5725/580608 9/1400 111587/2073600 2944/14175; - 0 -33953/29030400 -119/129600 -369/358400 -107/113400 -175/165888 -9/11200 -8183/4147200 989/28350] + elseif alg.M == 9 theta = [0.0 0.03188616071428564 0.031009210268469034 0.03117187499999874 0.031111111111105316 0.031149339236719698 0.031111111111073342 0.03117187499992724 0.031009210268393872 0.031886160714066136; 0.0 0.14467159330295903 0.18532725847540765 0.1828236607142859 0.18359396433471176 0.1831509191895293 0.18357142857178133 0.18292556155574857 0.18461297276007826 0.17568080357159488; @@ -159,20 +155,8 @@ function get_constant_parameters(alg::MPDeC) 0.0 0.04115018126592196 0.03318440133255063 0.036339285714291236 0.034175974916735186 0.03666654418941562 0.03142857142916 0.07940414951872299 0.18461297276189725 0.0120535714286234; 0.0 -0.008932169189692364 -0.007244757985499331 -0.007890625000002531 -0.007470115618263051 -0.00791316076327453 -0.007142857142902415 -0.009646454903545987 0.03100921026901915 0.17568080356954852; 0.0 0.0008769504458161903 0.0007142857142857367 0.000775049603174427 0.0007368214775622661 0.0007750496031686538 0.0007142857142810044 0.0008769504457859512 5.684341886080802e-14 0.03188616071417982] - elseif alg.M == 10 - theta = [0.0 0.028018959644393687 0.027340374645930202 0.027451045048698663 0.0274153171930962 0.02743427683751065 0.027418831168816382 0.027437790813053198 0.02740206295788994 0.027512733360708808 0.02683414836155862; - 0.0 0.13699028395729754 0.17247367858478835 0.17057771915585285 0.17108139597029037 0.17083711202639051 0.17102597402617903 0.17080197227005556 0.17121393832712783 0.16996083603953593 0.17753594143869122; - 0.0 -0.18583978613015117 -0.08617167708834472 -0.04460141030846643 -0.048462401795699606 -0.04691594453948866 -0.048009740260582134 -0.046778097819711206 -0.04896713164639266 -0.04246829342450553 -0.08104357062261158; - 0.0 0.30192072009780363 0.22804745871412313 0.3094549512987534 0.35692031425357507 0.34993098144275514 0.35402597402539016 0.34980383697848083 0.3569305756250287 0.33648092529256246 0.45494628807136905; - 0.0 -0.3806483247655124 -0.3026606541606571 -0.3400126826297907 -0.2703953823952361 -0.21667333679079093 -0.2287597402537358 -0.2184080650258693 -0.23442039440851659 -0.19077242288039997 -0.43515512242447585; - 0.0 0.35715424082090674 0.29001218534554596 0.31687012987021035 0.29602437069246434 0.3568823152170353 0.41774025973563766 0.3968945005565274 0.4237524451623358 0.35661038970238224 0.7137646304354348; - 0.0 -0.24438269976550997 -0.2007347282347136 -0.21674705762993796 -0.20639538239549893 -0.2184817858628776 -0.16475974026434415 -0.0951424400191172 -0.13249446850386448 -0.05450679792556912 -0.43515512242447585; - 0.0 0.11846536295494622 0.09801571268237676 0.10514245129872465 0.10092031425360393 0.10501530683775329 0.09802597402490676 0.14549133697096295 0.22689882957820373 0.15302556815186108 0.45494628810411086; - 0.0 -0.038575277201579265 -0.03207643899310476 -0.03426547280844172 -0.03303383036714891 -0.03412762608719788 -0.032581168831256946 -0.03644216031966607 0.005128106447045866 0.10479621548802243 -0.08104357064439682; - 0.0 0.007575105385869255 0.006322003099781336 0.006733969155847452 0.006509967398880434 0.0066988293985943415 0.006454545454857907 0.006958222269304315 0.005062262842329801 0.04054565746628214 0.17753594143141527; - 0.0 -0.0006785849984634729 -0.0005679145956923939 -0.000603642451298736 -0.0005846828069051568 -0.0006001284755665637 -0.000581168831175205 -0.0006168966868074222 -0.0005062262839601317 -0.0011848112826555735 0.026834148362013366] else - error("MPDeC requires 2 ≤ K ≤ 9.") + error("MPDeC requires 2 ≤ K ≤ 10.") end else # alg.nodes == :gausslobatto if alg.M == 1 @@ -211,88 +195,8 @@ function get_constant_parameters(alg::MPDeC) 0.0 0.007627676118250971 -0.024004074733154912 0.14346845286688126 0.2923037943068376 0.27742918851774334; 0.0 -0.004471780440573713 0.011807696377659743 -0.02460333048390262 0.10736966113994595 0.18923747814891745; 0.0 0.0016434260039545345 -0.004129309556393734 0.007424947945453564 -0.012346471800418257 0.03333333333333588] - elseif alg.M == 6 - nodes = [ - 0.0, - 0.08488805186071652, - 0.2655756032646429, - 0.5, - 0.7344243967353571, - 0.9151119481392835, - 1.0 - ] - theta = [0.0 0.03284626432829264 0.018002223201815104 0.027529761904763195 0.02170954269630787 0.02464779425215191 0.023809523809520172; - 0.0 0.059322894027551365 0.15770113064168867 0.1277882555598353 0.14409488824697547 0.13619592709181916 0.1384130236807124; - 0.0 -0.01076859445118926 0.10235481204686203 0.23748565272164512 0.20629511050420746 0.2193616205757678 0.21587269060495373; - 0.0 0.0055975917805697745 -0.018478259273458753 0.12190476190476146 0.2622877830829893 0.23821193202898883 0.24380952380948884; - 0.0 -0.0034889299708074674 0.009577580100741362 -0.02161296211671493 0.1135178785580706 0.22664128505616077 0.2158726906049253; - 0.0 0.002217096588914541 -0.005681864566224326 0.010624768120945982 -0.019288106960911655 0.07909012965322404 0.13841302368079056; - 0.0 -0.000838270442615105 0.0020999811132187685 -0.0037202380952379155 0.005807300607708399 -0.00903674051876635 0.02380952380952195] - elseif alg.M == 7 - nodes = [ - 0.0, - 0.06412992574519666, - 0.20414990928342885, - 0.3953503910487606, - 0.6046496089512394, - 0.7958500907165711, - 0.9358700742548034, - 1.0 - ] - theta = [0.0 0.024737514438875514 0.013258719822130019 0.021034356725210923 0.01577972028613317 0.019036721343624663 0.017385669593011244 0.01785714285716722; - 0.0 0.044892662602755 0.12064973282409763 0.09628017831208524 0.11096097521429726 0.10224716355590147 0.10657833455800869 0.10535211357202456; - 0.0 -0.008140767742389747 0.07999763578665159 0.1890416498443277 0.16114433643933257 0.17539401516496778 0.1687159595472414 0.17056134624175545; - 0.0 0.004254082548545937 -0.014480830962625313 0.10124595564769123 0.22436671754081416 0.19859744812260516 0.2089336024046844 0.20622939732966472; - 0.0 -0.002704205075015816 0.0076319492068338685 -0.018137320211455643 0.10498344168165086 0.22071022829197284 0.20197531478076014 0.20622939732958345; - 0.0 0.0018453866944819865 -0.004832668923134664 0.009417009802429988 -0.01848030360257269 0.09056371045502942 0.1787021139839453 0.17056134624169772; - 0.0 -0.0012262209861064882 0.00310495001592085 -0.005608861642537821 0.00907193525967287 -0.015297619252294226 0.060459450968835426 0.10535211357171193; - 0.0 0.0004714732640502682 -0.001179578486445107 0.002077422571009513 -0.003177213868071016 0.0045984230350075705 -0.006880371581776679 0.017857142857105046] - elseif alg.M == 8 - nodes = [ - 0.0, - 0.05012100229426997, - 0.16140686024463113, - 0.3184412680869109, - 0.5, - 0.6815587319130891, - 0.8385931397553689, - 0.94987899770573, - 1.0 - ] - theta = [0.0 0.019293838201043245 0.01018408040822739 0.01656936984357027 0.011990017361096061 0.01514127656190567 0.013175811859724718 0.01417408466352299 0.013888888888914153; - 0.0 0.03512552097762184 0.09508650844936217 0.07509351709620726 0.08786983372308654 0.07945805621039881 0.08459518591985216 0.0820139081133675 0.08274768078126726; - 0.0 -0.006364102418704803 0.0639319956284379 0.15288206102488378 0.12868222230659399 0.14236568211174472 0.13451403216049584 0.13834341694882824 0.13726935625072656; - 0.0 0.0033337771969983894 -0.011585731333842608 0.0840854786889147 0.18975808031592045 0.16521365191611892 0.17719508127352412 0.1717199563247931 0.17321425548699665; - 0.0 -0.002136847017608247 0.006149936900551798 -0.015130673172334241 0.09287981859410088 0.20089031036043536 0.17960970028764223 0.18789648420624872 0.1857596371874024; - 0.0 0.0014942991627282308 -0.003980825787182701 0.008000603570397224 -0.0165438248294123 0.08912877679762232 0.18479998682039422 0.16988047828840536 0.17321425548630032; - 0.0 -0.001074060699381661 0.0027553240894192185 -0.0050963258617406915 0.008587133943507297 -0.015612704774753183 0.0733373606215082 0.14363345866922828 0.13726935624890757; - 0.0 0.0007337726665392911 -0.0018475051393835457 0.0032896245700402282 -0.005122152942657721 0.007654163684208015 -0.012338827668713748 0.04762215980304063 0.08274768078103989; - 0.0 -0.00028519577496630263 0.0007130770290413452 -0.0012523876730271555 0.0018988715277794554 -0.002680480954676767 0.0037048084806841075 -0.005404949312023177 0.013888888889184159] - elseif alg.M == 9 - nodes = [ - 0.0, - 0.04023304591677057, - 0.13061306744724743, - 0.26103752509477773, - 0.4173605211668065, - 0.5826394788331936, - 0.7389624749052223, - 0.8693869325527526, - 0.9597669540832294, - 1.0 - ] - theta = [0.0 0.015465148886122208 0.008074564828863144 0.013376018525479871 0.009424872116881033 0.012319010010855891 0.010311035538052238 0.01156730916181914 0.010928591594165482 0.011111111110039928; - 0.0 0.02821797600073163 0.07677451467523239 0.060184542475466785 0.07119971523047752 0.06348355183104104 0.06872200531256567 0.06548282490166457 0.06711913714752882 0.0666529954141879; - 0.0 -0.0051090759506765785 0.05211803626519543 0.12565307727248143 0.10482643149413295 0.11734393918791586 0.10937232662209873 0.11414511256730009 0.11177491077136281 0.1124446710354885; - 0.0 0.0026797026134601862 -0.009449705457448665 0.0703555323098648 0.1607094677113423 0.1383496630545551 0.15043195838711654 0.14368298780891564 0.14692275221750606 0.1460213418404237; - 0.0 -0.0017246233123523063 0.005034230848115648 -0.012689136533024836 0.08096672386450621 0.17827052516622643 0.15700518309495237 0.1670603392740304 0.16255069057842775 0.1637698805952823; - 0.0 0.0012191900098376621 -0.003290458683960145 0.006764697496483674 -0.014500644574205523 0.08280315672743654 0.17645901712467094 0.15873564974754117 0.16549450390448328 0.1637698805971013; - 0.0 -0.0009014103789669575 0.002338354035534229 -0.004410616548205404 0.007671678785357017 -0.014688125871543889 0.07566580953016455 0.15547104729739658 0.1433416392246727 0.1460213418395142; - 0.0 0.0006697602581937584 -0.00170044153617834 0.0030723444103197828 -0.0048992681562261 0.007618239537521276 -0.013208406240892145 0.06032663476707967 0.11755374698441301 0.1124446710346092; - 0.0 -0.000466141727447515 0.0011701705230945658 -0.002069009887745693 0.0031694435943515065 -0.004546719804944366 0.0064684529496542575 -0.010121519251002686 0.03843501942257066 0.0666529954255306; - 0.0 0.0001825195178684848 -0.0004561980512006905 0.0008000755736511378 -0.0012078988997608064 0.00168623899423892 -0.002264907414485151 0.00303654628237382 -0.004354037773737218 0.011111111112313665] else - error("MPDeC requires 2 ≤ K ≤ 9.") + error("MPDeC requires 2 ≤ K ≤ 10.") end end return nodes, theta @@ -440,34 +344,6 @@ end end end -#= -function _build_mpdec_matrix_and_rhs_old!(Mmat, rhs, P, dt_th, σ, d = nothing) - N = length(rhs) - @fastmath @inbounds @simd for i in 1:N - # Add nonconservative destruction terms to diagonal (PDSFunctions only!) - if !isnothing(d) - if dt_th >= 0 - Mmat[i, i] += dt_th * d[i] / σ[i] - rhs[i] += dt_th * P[i, i] - else - Mmat[i, i] -= dt_th * P[i, i] / σ[i] - rhs[i] -= dt_th * d[i] - end - end - @fastmath @inbounds @simd for j in 1:N - if j != i - if dt_th >= 0 - Mmat[i, j] -= dt_th * P[i, j] / σ[j] - Mmat[i, i] += dt_th * P[j, i] / σ[i] # P[j, i] = D[i, j] - else - Mmat[i, j] += dt_th * P[j, i] / σ[j] # P[j, i] = D[i, j] - Mmat[i, i] -= dt_th * P[i, j] / σ[i] - end - end - end - end -end -=# function _build_mpdec_matrix_and_rhs!(M, rhs, P, dt_th, σ, d = nothing) Base.require_one_based_indexing(M, P, σ) @assert size(M, 1) == size(M, 2) == size(P, 1) == size(P, 2) == length(σ) @@ -580,9 +456,6 @@ function _build_mpdec_matrix_and_rhs!(M::AbstractSparseMatrix, rhs, P::AbstractS @assert length(σ) == length(d) end - #TODO: It is crucial that evaluation of the production function - # does not alter the sparsity pattern of p_prototype. This should be checked earlier. - # By construction M and P share the same sparsity pattern. M_rows = rowvals(M) M_vals = nonzeros(M) @@ -701,11 +574,6 @@ end u = C2[:, M + 1] u1 = C[:, M + 1] # one order less accurate - #TODO: Remove this check - # Check only valid for autonomous conserverative PDS - #u_check = MPDeC_check(K, M, uprev, theta, f, p, t, dt) - #@assert u ≈ u_check - tmp = u - u1 atmp = calculate_residuals(tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t) @@ -753,10 +621,18 @@ function alg_cache(alg::MPDeC, u, rate_prototype, ::Type{uEltypeNoUnits}, tmp = zero(u) P = p_prototype(u, f) # stores evaluation of the production matrix P2 = p_prototype(u, f) # stores the linear system matrix - if issparse(P2) && alg.K > 2 - # Negative weights of MPDeC(K) , K >=3 require - # a symmetric sparsity pattern - P2 = P2 + P2' + if issparse(P2) + # We need to ensure that evaluating the production function does + # not alter the sparsity pattern given by the production matrix prototype + f.p(P2, uprev, p, t) + @assert P.rowval == P2.rowval&&P.colptr == P2.colptr "Evaluation of the production terms must not alter the sparsity pattern given by the prototype." + + if alg.K > 2 + # Negative weights of MPDeC(K) , K > 2 require + # a symmetric sparsity pattern of the linear system matrix P2 + P2 = P2 + P2' + end + else end d = zero(u) σ = zero(u) @@ -877,83 +753,3 @@ end False()) integrator.EEst = integrator.opts.internalnorm(tmp, t) end - -# TODO: Remove the code below after testing is complete -######################################################################################################## -######################################################################################################## -#### The functions below are provided with the original paper and used for validation. #### -#### They will be removed in the future. #### -######################################################################################################## -######################################################################################################## -function patanker_type_dec(prod_p, dest_p, delta_t, m_substep::Int, M_sub, theta, u_p, dim) - #This function builds and solves the system u^{(k+1)}=Mu^{(k)} for a subtimestep m - - mass = Matrix{Float64}(I, dim, dim) - #println(mass) - for i in 1:dim - for r in 1:(M_sub + 1) - if theta[r, m_substep] > 0 - for j in 1:dim - mass[i, j] = mass[i, j] - - delta_t * theta[r, m_substep] * - (prod_p[i, j, r] / u_p[j, m_substep]) - mass[i, i] = mass[i, i] + - delta_t * theta[r, m_substep] * - (dest_p[i, j, r] / u_p[i, m_substep]) - end - elseif theta[r, m_substep] < 0 - for j in 1:dim - mass[i, i] = mass[i, i] - - delta_t * theta[r, m_substep] * - (prod_p[i, j, r] / u_p[i, m_substep]) - mass[i, j] = mass[i, j] + - delta_t * theta[r, m_substep] * - (dest_p[i, j, r] / u_p[j, m_substep]) - end - end - end - end - return mass \ u_p[:, 1] -end - -function MPDeC_check(K, M, uprev, theta, f, p, t, dt) - M_sub = M - K_corr = K - dim = length(uprev) - - #Variables of the DeC procedure: one for each subtimesteps and one for previous correction up, one for current correction ua - u_p = zeros(dim, M_sub + 1) - u_a = zeros(dim, M_sub + 1) - - #Matrices of production and destructions applied to up at different subtimesteps - prod_p = zeros(dim, dim, M_sub + 1) - dest_p = zeros(dim, dim, M_sub + 1) - - #coefficients for time integration - Theta = theta - - # Subtimesteps loop to update variables at iteration 0 - for m in 1:(M_sub + 1) - u_a[:, m] = uprev - u_p[:, m] = uprev - end - - #Loop for iterations K - for k in 2:(K_corr + 1) - # Updating previous iteration variabls - u_p = copy(u_a) - #Loop on subtimesteps to compute the production destruction terms - for r in 1:(M_sub + 1) - #prod_p[:,:,r], dest_p[:,:,r]=prod_dest(u_p[:,r]) - prod_p[:, :, r] = f.p(u_p[:, r], p, t) - dest_p[:, :, r] = prod_p[:, :, r]' - end - # Loop on subtimesteps to compute the new variables - for m in 2:(M_sub + 1) - u_a[:, m] = patanker_type_dec(prod_p, dest_p, dt, m, M_sub, Theta, u_p, - dim) - end - end - u1 = u_a[:, M_sub + 1] - return u1 -end diff --git a/test/runtests.jl b/test/runtests.jl index 76ad2bc0..f3c8cca3 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1101,6 +1101,19 @@ end @test_throws "MPDeC requires the parameter K to be an integer." solve(prob_pds_linmod, MPDeC(2.1), dt = 0.1) + @test_throws "MPDeC requires 2 ≤ K ≤ 10." solve(prob_pds_linmod, MPDeC(123)) + @test_throws "MPDeC requires 2 ≤ K ≤ 10." solve(prob_pds_linmod, MPDeC(123, nodes =:lagrange)) + P = spdiagm(1 => [1.0]) + function prod!(P, u, p, t) + P[2, 1] = one(eltype(P)) + end + prob = ConservativePDSProblem(prod!, ones(2, 1), (0.0, 1.0); p_prototype = P) + @test_throws "must not alter the sparsity pattern" solve(prob, MPDeC(2)) + end + + # Check parameter type conversions + @testset "Parameter input" begin + @test_nowarn solve(prob_pds_linmod, MPDeC(2.0)) end # Here we check that algorithms which accept input parameters return constants From 4759126af812f51fa33a252c52f1970f9176e41a Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Wed, 19 Mar 2025 12:26:16 +0100 Subject: [PATCH 21/38] format --- test/runtests.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index f3c8cca3..a7128445 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1102,7 +1102,8 @@ end MPDeC(2.1), dt = 0.1) @test_throws "MPDeC requires 2 ≤ K ≤ 10." solve(prob_pds_linmod, MPDeC(123)) - @test_throws "MPDeC requires 2 ≤ K ≤ 10." solve(prob_pds_linmod, MPDeC(123, nodes =:lagrange)) + @test_throws "MPDeC requires 2 ≤ K ≤ 10." solve(prob_pds_linmod, + MPDeC(123, nodes = :lagrange)) P = spdiagm(1 => [1.0]) function prod!(P, u, p, t) P[2, 1] = one(eltype(P)) From 35ff52a29ef9d493c1619d050a489d3c322f5ea5 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Wed, 19 Mar 2025 15:32:04 +0100 Subject: [PATCH 22/38] fix - change in sparsity pattern --- src/mpdec.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/mpdec.jl b/src/mpdec.jl index d0d5baf2..3c8ef7ec 100644 --- a/src/mpdec.jl +++ b/src/mpdec.jl @@ -630,9 +630,8 @@ function alg_cache(alg::MPDeC, u, rate_prototype, ::Type{uEltypeNoUnits}, if alg.K > 2 # Negative weights of MPDeC(K) , K > 2 require # a symmetric sparsity pattern of the linear system matrix P2 - P2 = P2 + P2' + P2 = P + P' end - else end d = zero(u) σ = zero(u) From e63d9a58cde695a30aee12c19f9011c12f399dd4 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Thu, 20 Mar 2025 13:13:14 +0100 Subject: [PATCH 23/38] Use POSITIVE PDS in tests --- test/runtests.jl | 42 +++++++++++++++++++++--------------------- 1 file changed, 21 insertions(+), 21 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index a7128445..13d33e40 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -343,22 +343,22 @@ end prod_1! = (P, u, p, t) -> begin fill!(P, zero(eltype(P))) for i in 1:(length(u) - 1) - P[i, i + 1] = i * u[i] + P[i, i + 1] = i * u[i + 1] end return nothing end prod_2! = (P, u, p, t) -> begin fill!(P, zero(eltype(P))) for i in 1:(length(u) - 1) - P[i + 1, i] = i * u[i + 1] + P[i + 1, i] = i * u[i] end return nothing end prod_3! = (P, u, p, t) -> begin fill!(P, zero(eltype(P))) for i in 1:(length(u) - 1) - P[i, i + 1] = i * u[i] - P[i + 1, i] = i * u[i + 1] + P[i, i + 1] = i * u[i + 1] + P[i + 1, i] = i * u[i] end return nothing end @@ -1317,7 +1317,7 @@ end prod_1! = (P, u, p, t) -> begin fill!(P, zero(eltype(P))) for i in 1:(length(u) - 1) - P[i, i + 1] = i * u[i] + P[i, i + 1] = i * u[i + 1] end return nothing end @@ -1325,7 +1325,7 @@ end prod_2! = (P, u, p, t) -> begin fill!(P, zero(eltype(P))) for i in 1:(length(u) - 1) - P[i + 1, i] = i * u[i + 1] + P[i + 1, i] = i * u[i] end return nothing end @@ -1333,8 +1333,8 @@ end prod_3! = (P, u, p, t) -> begin fill!(P, zero(eltype(P))) for i in 1:(length(u) - 1) - P[i, i + 1] = i * u[i] - P[i + 1, i] = i * u[i + 1] + P[i, i + 1] = i * u[i + 1] + P[i + 1, i] = i * u[i] end return nothing end @@ -1412,7 +1412,7 @@ end prod_1! = (P, u, p, t) -> begin fill!(P, zero(eltype(P))) for i in 1:(length(u) - 1) - P[i, i + 1] = i * u[i] + P[i, i + 1] = i * u[i + 1] end return nothing end @@ -1420,7 +1420,7 @@ end prod_2! = (P, u, p, t) -> begin fill!(P, zero(eltype(P))) for i in 1:(length(u) - 1) - P[i + 1, i] = i * u[i + 1] + P[i + 1, i] = i * u[i] end return nothing end @@ -1428,8 +1428,8 @@ end prod_3! = (P, u, p, t) -> begin fill!(P, zero(eltype(P))) for i in 1:(length(u) - 1) - P[i, i + 1] = i * u[i] - P[i + 1, i] = i * u[i + 1] + P[i, i + 1] = i * u[i + 1] + P[i + 1, i] = i * u[i] end return nothing end @@ -1502,7 +1502,7 @@ end prod_1! = (P, u, p, t) -> begin fill!(P, zero(eltype(P))) for i in 1:(length(u) - 1) - P[i, i + 1] = i * u[i] + P[i, i + 1] = i * u[i + 1] end for i in 1:length(u) P[i, i] = i * u[i] @@ -1520,7 +1520,7 @@ end prod_2! = (P, u, p, t) -> begin fill!(P, zero(eltype(P))) for i in 1:(length(u) - 1) - P[i + 1, i] = i * u[i + 1] + P[i + 1, i] = i * u[i] end for i in 1:length(u) P[i, i] = (i - 1) * u[i] @@ -1538,8 +1538,8 @@ end prod_3! = (P, u, p, t) -> begin fill!(P, zero(eltype(P))) for i in 1:(length(u) - 1) - P[i, i + 1] = i * u[i] - P[i + 1, i] = i * u[i + 1] + P[i, i + 1] = i * u[i + 1] + P[i + 1, i] = i * u[i] end for i in 1:length(u) P[i, i] = (i + 1) * u[i] @@ -1633,7 +1633,7 @@ end prod_1! = (P, u, p, t) -> begin fill!(P, zero(eltype(P))) for i in 1:(length(u) - 1) - P[i, i + 1] = i * u[i] + P[i, i + 1] = i * u[i + 1] end for i in 1:length(u) P[i, i] = i * u[i] @@ -1651,7 +1651,7 @@ end prod_2! = (P, u, p, t) -> begin fill!(P, zero(eltype(P))) for i in 1:(length(u) - 1) - P[i + 1, i] = i * u[i + 1] + P[i + 1, i] = i * u[i] end for i in 1:length(u) P[i, i] = (i - 1) * u[i] @@ -1669,8 +1669,8 @@ end prod_3! = (P, u, p, t) -> begin fill!(P, zero(eltype(P))) for i in 1:(length(u) - 1) - P[i, i + 1] = i * u[i] - P[i + 1, i] = i * u[i + 1] + P[i, i + 1] = i * u[i + 1] + P[i + 1, i] = i * u[i] end for i in 1:length(u) P[i, i] = (i + 1) * u[i] @@ -1770,7 +1770,7 @@ end prod_inner! = (P, u, p, t) -> begin fill!(P, zero(eltype(P))) for i in 1:(length(u) - 1) - P[i, i + 1] = i * u[i] + P[i, i + 1] = i * u[i + 1] end return nothing end From 8d9e12871fb0fa258af66594d621e7ea3ea2c819 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Fri, 21 Mar 2025 14:29:29 +0100 Subject: [PATCH 24/38] Apply suggestions from code review Co-authored-by: Hendrik Ranocha --- src/mpdec.jl | 46 ++++++++++++++++++++++------------------------ 1 file changed, 22 insertions(+), 24 deletions(-) diff --git a/src/mpdec.jl b/src/mpdec.jl index 3c8ef7ec..fba0c498 100644 --- a/src/mpdec.jl +++ b/src/mpdec.jl @@ -6,7 +6,7 @@ production-destruction systems. Each member of this family is an adaptive, one-s Kth order accurate, unconditionally positivity-preserving, and linearly implicit. The integer K must be chosen to satisfy 2 ≤ K ≤ 10. Available node choices are Lagrange or Gauss-Lobatto nodes, with the latter being the default. -These methods support adaptive time stepping using the numerical solution obtained with one correction step less, as lower order approximation to estimate the error. +These methods support adaptive time stepping, using the numerical solution obtained with one correction step less as a lower-order approximation to estimate the error. The MPDeC schemes were introduced by Torlo and Öffner (2020) for autonomous conservative production-destruction systems and further investigated in Torlo, Öffner and Ranocha (2022). @@ -23,25 +23,27 @@ algorithm from [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl) as keyword argument `linsolve`. You can also choose the parameter `small_constant` which is added to all Patankar-weight denominators to avoid divisions by zero. You can pass a value explicitly, otherwise `small_constant` is set to -1e-300 in double precision computations or `floatmin` of the floating point type used. +`1e-300` in double precision computations or `floatmin` of the floating point type used. ## References -- Davide Torlo, and Philipp Öffner. +- Davide Torlo and Philipp Öffner. "Arbitrary high-order, conservative and positivity preserving Patankar-type deferred correction schemes." Applied Numerical Mathematics 153 (2020): 15-34. + [DOI: 10.1016/j.apnum.2020.01.025](https://doi.org/10.1016/j.apnum.2020.01.025) - Davide Torlo, Philipp Öffner, and Hendrik Ranocha. "Issues with positivity-preserving Patankar-type schemes." Applied Numerical Mathematics 182 (2022): 117-147. [DOI: 10.1016/j.apnum.2022.07.014](https://doi.org/10.1016/j.apnum.2022.07.014) -- Benjamin W. Ong & Raymond J. Spiteri. +- Benjamin W. Ong and Raymond J. Spiteri. "Deferred Correction Methods for Ordinary Differential Equations." - Journal of Scientific Computing 83 (2020): Article 60 + Journal of Scientific Computing 83 (2020): Article 60. + [DOI: 10.1007/s10915-020-01235-8](https://doi.org/10.1007/s10915-020-01235-8) """ -struct MPDeC{T, N, F, T2} <: OrdinaryDiffEqAdaptiveAlgorithm - K::T - M::T +struct MPDeC{N, F, T2} <: OrdinaryDiffEqAdaptiveAlgorithm + K::Int + M::Int nodes::N linsolve::F small_constant_function::T2 @@ -58,14 +60,10 @@ function small_constant_function_MPDeC(type) return small_constant end -function MPDeC(K; nodes = :gausslobatto, linsolve = LUFactorization(), +function MPDeC(K::Integer; + nodes = :gausslobatto, + linsolve = LUFactorization(), small_constant = small_constant_function_MPDeC) - if !(isinteger(K)) - throw(ArgumentError("MPDeC requires the parameter K to be an integer.")) - end - if !(typeof(K) <: Integer) - K = Int(K) - end if small_constant isa Number small_constant_function = Returns(small_constant) @@ -73,16 +71,16 @@ function MPDeC(K; nodes = :gausslobatto, linsolve = LUFactorization(), small_constant_function = small_constant end - if nodes == :lagrange + if nodes === :lagrange M = K - 1 else # :gausslobatto - M = ceil(Integer, K / 2) + M = cld(K, 2) end - MPDeC{typeof(K), typeof(nodes), typeof(linsolve), typeof(small_constant_function)}(K, M, - nodes, - linsolve, - small_constant_function) + MPDeC{typeof(nodes), typeof(linsolve), typeof(small_constant_function)}(K, M, + nodes, + linsolve, + small_constant_function) end alg_order(alg::MPDeC) = alg.K @@ -202,9 +200,9 @@ function get_constant_parameters(alg::MPDeC) return nodes, theta end -struct MPDeCConstantCache{KType, NType, T, T2} <: OrdinaryDiffEqConstantCache - K::KType - M::KType +struct MPDeCConstantCache{NType, T, T2} <: OrdinaryDiffEqConstantCache + K::Int + M::Int nodes::NType theta::T2 small_constant::T From a1cf60e697fc3bb14c4962567d50de5e8d0c8a74 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Fri, 21 Mar 2025 14:39:57 +0100 Subject: [PATCH 25/38] additiional changes after review --- src/mpdec.jl | 14 ++++++-------- test/runtests.jl | 47 ++++++++++++++++------------------------------- 2 files changed, 22 insertions(+), 39 deletions(-) diff --git a/src/mpdec.jl b/src/mpdec.jl index fba0c498..a617a5e5 100644 --- a/src/mpdec.jl +++ b/src/mpdec.jl @@ -64,7 +64,6 @@ function MPDeC(K::Integer; nodes = :gausslobatto, linsolve = LUFactorization(), small_constant = small_constant_function_MPDeC) - if small_constant isa Number small_constant_function = Returns(small_constant) else # assume small_constant isa Function @@ -454,14 +453,13 @@ function _build_mpdec_matrix_and_rhs!(M::AbstractSparseMatrix, rhs, P::AbstractS @assert length(σ) == length(d) end - # By construction M and P share the same sparsity pattern. M_rows = rowvals(M) M_vals = nonzeros(M) P_rows = rowvals(P) P_vals = nonzeros(P) n = size(M, 2) - # tmp[j] = M[j,j] + # We use tmp as a buffer for the diagonal elements of M. fill!(tmp, zero(eltype(tmp))) if dt_th ≥ 0 @@ -473,7 +471,7 @@ function _build_mpdec_matrix_and_rhs!(M::AbstractSparseMatrix, rhs, P::AbstractS for idx_M in nzrange(M, j) if M_rows[idx_M] == i M_vals[idx_M] -= dt_th_P / σ[j] # M_ij <- P_ij - #break + break end end tmp[j] += dt_th_P / σ[j] # M_jj <- P_ij = D_ji @@ -494,7 +492,7 @@ function _build_mpdec_matrix_and_rhs!(M::AbstractSparseMatrix, rhs, P::AbstractS i = M_rows[idx_M] if i == j M_vals[idx_M] += tmp[j] - #break + break end end end @@ -508,14 +506,14 @@ function _build_mpdec_matrix_and_rhs!(M::AbstractSparseMatrix, rhs, P::AbstractS if M_rows[idx_M] == j M_vals[idx_M] += dt_th_P / σ[i] # M_ji <- P_ij end - #break + break end tmp[i] -= dt_th_P / σ[i] else for idx_M in nzrange(M, j) if i == M_rows[idx_M] M_vals[idx_M] -= dt_th_P / σ[i] # M_ij <- P_ij - #break + break end end end @@ -527,7 +525,7 @@ function _build_mpdec_matrix_and_rhs!(M::AbstractSparseMatrix, rhs, P::AbstractS i = M_rows[idx_M] if i == j M_vals[idx_M] += tmp[j] - #break + break end end end diff --git a/test/runtests.jl b/test/runtests.jl index 13d33e40..cfba8af6 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1112,11 +1112,6 @@ end @test_throws "must not alter the sparsity pattern" solve(prob, MPDeC(2)) end - # Check parameter type conversions - @testset "Parameter input" begin - @test_nowarn solve(prob_pds_linmod, MPDeC(2.0)) - end - # Here we check that algorithms which accept input parameters return constants # of the same type as the inputs @testset "Constant types" begin @@ -1408,7 +1403,6 @@ end end @testset "Different matrix types (conservative, adaptive)" begin - #TODO: MPE, SSPMPRK43 are not adaptive and need to be removed (solve without giving dt!) prod_1! = (P, u, p, t) -> begin fill!(P, zero(eltype(P))) for i in 1:(length(u) - 1) @@ -1446,11 +1440,10 @@ end rtol = sqrt(eps(Float32)) - algs = [MPE(), - MPRK22(0.5), MPRK22(1.0), + algs = [MPRK22(0.5), MPRK22(1.0), MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), MPRK43II(2.0 / 3.0), MPRK43II(0.5), - SSPMPRK22(0.5, 1.0), SSPMPRK43()] + SSPMPRK22(0.5, 1.0)] for k in 2:10 push!(algs, MPDeC(k), MPDeC(k; nodes = :lagrange)) end @@ -1474,12 +1467,12 @@ end prob_sparse_op = ConservativePDSProblem(prod, u0, tspan; p_prototype = P_sparse) - sol_tridiagonal_ip = solve(prob_tridiagonal_ip, alg; dt) - sol_tridiagonal_op = solve(prob_tridiagonal_op, alg; dt) - sol_dense_ip = solve(prob_dense_ip, alg; dt) - sol_dense_op = solve(prob_dense_op, alg; dt) - sol_sparse_ip = solve(prob_sparse_ip, alg; dt) - sol_sparse_op = solve(prob_sparse_op, alg; dt) + sol_tridiagonal_ip = solve(prob_tridiagonal_ip, alg) + sol_tridiagonal_op = solve(prob_tridiagonal_op, alg) + sol_dense_ip = solve(prob_dense_ip, alg) + sol_dense_op = solve(prob_dense_op, alg) + sol_sparse_ip = solve(prob_sparse_ip, alg) + sol_sparse_op = solve(prob_sparse_op, alg) @test isapprox(sol_tridiagonal_ip.t, sol_tridiagonal_op.t; rtol) @test isapprox(sol_dense_ip.t, sol_dense_op.t; rtol) @@ -1629,7 +1622,6 @@ end end @testset "Different matrix types (nonconservative, adaptive)" begin - #TODO: MPE, SSPMPRK43 are not adaptive and need to be removed (solve without giving dt) prod_1! = (P, u, p, t) -> begin fill!(P, zero(eltype(P))) for i in 1:(length(u) - 1) @@ -1696,11 +1688,10 @@ end tspan = (0.0, 1.0) dt = 0.25 - algs = [MPE(), - MPRK22(0.5), MPRK22(1.0), + algs = [MPRK22(0.5), MPRK22(1.0), MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), MPRK43II(2.0 / 3.0), MPRK43II(0.5), - SSPMPRK22(0.5, 1.0), SSPMPRK43()] + SSPMPRK22(0.5, 1.0)] for k in 2:10 push!(algs, MPDeC(k), MPDeC(k; nodes = :lagrange)) end @@ -1734,18 +1725,12 @@ end prob_sparse_op = PDSProblem(prod, dest, u0, tspan; p_prototype = P_sparse) - sol_tridiagonal_ip = solve(prob_tridiagonal_ip, alg; - dt) - sol_tridiagonal_op = solve(prob_tridiagonal_op, alg; - dt) - sol_dense_ip = solve(prob_dense_ip, alg; - dt) - sol_dense_op = solve(prob_dense_op, alg; - dt) - sol_sparse_ip = solve(prob_sparse_ip, alg; - dt) - sol_sparse_op = solve(prob_sparse_op, alg; - dt) + sol_tridiagonal_ip = solve(prob_tridiagonal_ip, alg) + sol_tridiagonal_op = solve(prob_tridiagonal_op, alg) + sol_dense_ip = solve(prob_dense_ip, alg) + sol_dense_op = solve(prob_dense_op, alg) + sol_sparse_ip = solve(prob_sparse_ip, alg) + sol_sparse_op = solve(prob_sparse_op, alg) @test isapprox(sol_tridiagonal_ip.t, sol_tridiagonal_op.t; rtol) @test isapprox(sol_dense_ip.t, sol_dense_op.t; rtol) From 03ee723fe2f55a13857d909098dbb43c874fdcf8 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Fri, 21 Mar 2025 18:14:57 +0100 Subject: [PATCH 26/38] Removed breaks in loops for sparse matrices --- src/mpdec.jl | 10 +++++----- test/runtests.jl | 3 --- 2 files changed, 5 insertions(+), 8 deletions(-) diff --git a/src/mpdec.jl b/src/mpdec.jl index a617a5e5..464c032f 100644 --- a/src/mpdec.jl +++ b/src/mpdec.jl @@ -471,7 +471,7 @@ function _build_mpdec_matrix_and_rhs!(M::AbstractSparseMatrix, rhs, P::AbstractS for idx_M in nzrange(M, j) if M_rows[idx_M] == i M_vals[idx_M] -= dt_th_P / σ[j] # M_ij <- P_ij - break + #break end end tmp[j] += dt_th_P / σ[j] # M_jj <- P_ij = D_ji @@ -492,7 +492,7 @@ function _build_mpdec_matrix_and_rhs!(M::AbstractSparseMatrix, rhs, P::AbstractS i = M_rows[idx_M] if i == j M_vals[idx_M] += tmp[j] - break + #break end end end @@ -506,14 +506,14 @@ function _build_mpdec_matrix_and_rhs!(M::AbstractSparseMatrix, rhs, P::AbstractS if M_rows[idx_M] == j M_vals[idx_M] += dt_th_P / σ[i] # M_ji <- P_ij end - break + #break end tmp[i] -= dt_th_P / σ[i] else for idx_M in nzrange(M, j) if i == M_rows[idx_M] M_vals[idx_M] -= dt_th_P / σ[i] # M_ij <- P_ij - break + #break end end end @@ -525,7 +525,7 @@ function _build_mpdec_matrix_and_rhs!(M::AbstractSparseMatrix, rhs, P::AbstractS i = M_rows[idx_M] if i == j M_vals[idx_M] += tmp[j] - break + #break end end end diff --git a/test/runtests.jl b/test/runtests.jl index cfba8af6..ccf2a8a3 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1098,9 +1098,6 @@ end @test_throws "MPDeC can only be applied to production-destruction systems" solve(prob_oop, MPDeC(2), dt = 0.1) - @test_throws "MPDeC requires the parameter K to be an integer." solve(prob_pds_linmod, - MPDeC(2.1), - dt = 0.1) @test_throws "MPDeC requires 2 ≤ K ≤ 10." solve(prob_pds_linmod, MPDeC(123)) @test_throws "MPDeC requires 2 ≤ K ≤ 10." solve(prob_pds_linmod, MPDeC(123, nodes = :lagrange)) From f1b32a4bb191a9ed0d73587b40d4a236a9f085fc Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Fri, 21 Mar 2025 18:30:46 +0100 Subject: [PATCH 27/38] fixed breaks in loops for sparse matrices --- src/mpdec.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/mpdec.jl b/src/mpdec.jl index 464c032f..b4a5ff06 100644 --- a/src/mpdec.jl +++ b/src/mpdec.jl @@ -471,7 +471,7 @@ function _build_mpdec_matrix_and_rhs!(M::AbstractSparseMatrix, rhs, P::AbstractS for idx_M in nzrange(M, j) if M_rows[idx_M] == i M_vals[idx_M] -= dt_th_P / σ[j] # M_ij <- P_ij - #break + break end end tmp[j] += dt_th_P / σ[j] # M_jj <- P_ij = D_ji @@ -492,7 +492,7 @@ function _build_mpdec_matrix_and_rhs!(M::AbstractSparseMatrix, rhs, P::AbstractS i = M_rows[idx_M] if i == j M_vals[idx_M] += tmp[j] - #break + break end end end @@ -505,15 +505,15 @@ function _build_mpdec_matrix_and_rhs!(M::AbstractSparseMatrix, rhs, P::AbstractS for idx_M in nzrange(M, i) if M_rows[idx_M] == j M_vals[idx_M] += dt_th_P / σ[i] # M_ji <- P_ij + break end - #break end tmp[i] -= dt_th_P / σ[i] else for idx_M in nzrange(M, j) if i == M_rows[idx_M] M_vals[idx_M] -= dt_th_P / σ[i] # M_ij <- P_ij - #break + break end end end @@ -525,7 +525,7 @@ function _build_mpdec_matrix_and_rhs!(M::AbstractSparseMatrix, rhs, P::AbstractS i = M_rows[idx_M] if i == j M_vals[idx_M] += tmp[j] - #break + break end end end From 3f7acd4562def03535053a8b820974d53d582dc5 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Sun, 23 Mar 2025 15:17:06 +0100 Subject: [PATCH 28/38] use algebraic expressions for integration nodes and weights --- src/mpdec.jl | 188 ++++++++++++++++++++++++++++----------------------- 1 file changed, 102 insertions(+), 86 deletions(-) diff --git a/src/mpdec.jl b/src/mpdec.jl index b4a5ff06..8962f399 100644 --- a/src/mpdec.jl +++ b/src/mpdec.jl @@ -85,113 +85,129 @@ end alg_order(alg::MPDeC) = alg.K isfsal(::MPDeC) = false -function get_constant_parameters(alg::MPDeC) +function get_constant_parameters(alg::MPDeC, type) + oneType = one(type) if alg.nodes == :lagrange - nodes = collect(0:(1 / alg.M):1) + nodes = collect(i * oneType / alg.M for i in 0:(alg.M)) # M is one less than the methods order if alg.M == 1 - theta = [0 0.5; 0 0.5] + theta = [0 oneType/2; 0 oneType/2] elseif alg.M == 2 - theta = [0.0 0.20833333333333337 0.16666666666666652; - 0.0 0.33333333333333337 0.6666666666666667; - 0.0 -0.04166666666666667 0.16666666666666663] + theta = [0 5oneType/24 oneType/6; + 0 oneType/3 2oneType/3; + 0 -oneType/24 oneType/6] elseif alg.M == 3 - theta = [0.0 0.125 0.11111111111111116 0.125; - 0.0 0.26388888888888895 0.4444444444444451 0.3750000000000009; - 0.0 -0.0694444444444445 0.11111111111111072 0.375; - 0.0 0.013888888888888888 5.551115123125783e-17 0.12499999999999978] + theta = [0 1oneType/8 1oneType/9 1oneType/8; + 0 19oneType/72 4oneType/9 3oneType/8; + 0 -5oneType/72 1oneType/9 3oneType/8; + 0 1oneType/72 0 1oneType/8] elseif alg.M == 4 - theta = [0.0 0.08715277777777783 0.0805555555555556 0.08437500000000009 0.07777777777777928; - 0.0 0.2243055555555556 0.34444444444444455 0.31874999999999964 0.35555555555555785; - 0.0 -0.09166666666666667 0.06666666666666643 0.22499999999999964 0.13333333333333286; - 0.0 0.036805555555555564 0.011111111111111016 0.1312500000000001 0.35555555555555607; - 0.0 -0.0065972222222222265 -0.0027777777777778234 -0.009374999999999911 0.0777777777777775] + theta = th2 = [0 251oneType/2880 29oneType/360 27oneType/320 7oneType/90; + 0 323oneType/1440 31oneType/90 51oneType/160 16oneType/45; + 0 -11oneType/120 1oneType/15 9oneType/40 2oneType/15; + 0 53oneType/1440 1oneType/90 21oneType/160 16oneType/45; + 0 -19oneType/2880 -1oneType/360 -3oneType/320 7oneType/90] elseif alg.M == 5 - theta = [0.0 0.06597222222222225 0.0622222222222224 0.06374999999999953 0.06222222222222573 0.06597222222222676; - 0.0 0.1981944444444445 0.28666666666666707 0.2737500000000015 0.28444444444444983 0.2604166666666785; - 0.0 -0.1108333333333335 0.031111111111109757 0.14249999999999297 0.10666666666664959 0.1736111111111054; - 0.0 0.06694444444444442 0.031111111111110645 0.1424999999999983 0.2844444444444356 0.17361111111109295; - 0.0 -0.024027777777777766 -0.013333333333332975 -0.026250000000000107 0.06222222222222484 0.2604166666666785; - 0.0 0.0037500000000000033 0.0022222222222222365 0.003750000000000364 0.0 0.06597222222222454] + theta = [0 19oneType/288 14oneType/225 51oneType/800 14oneType/225 19oneType/288; + 0 1427oneType/7200 43oneType/150 219oneType/800 64oneType/225 25oneType/96; + 0 -133oneType/1200 7oneType/225 57oneType/400 8oneType/75 25oneType/144; + 0 241oneType/3600 7oneType/225 57oneType/400 64oneType/225 25oneType/144; + 0 -173oneType/7200 -1oneType/75 -21oneType/800 14oneType/225 25oneType/96; + 0 3oneType/800 1oneType/450 3oneType/800 0 19oneType/288] elseif alg.M == 6 - theta = [0.0 0.05259865520282189 0.05022045855379187 0.05096726190476186 0.050440917107582806 0.05118772045855735 0.04880952380950987; - 0.0 0.179431216931217 0.2486772486772483 0.24107142857142755 0.24550264550265144 0.239748677248647 0.25714285714280294; - 0.0 -0.12803406084656088 0.0014550264550234893 0.08638392857142296 0.06772486772487873 0.08783895502639893 0.0321428571428477; - 0.0 0.1033509700176369 0.058553791887126866 0.16190476190476288 0.26525573192239804 0.22045855379185753 0.3238095238095582; - 0.0 -0.055696097883597806 -0.03558201058201127 -0.05424107142856682 0.030687830687839757 0.16017691798938927 0.03214285714290277; - 0.0 0.017394179894179934 0.011640211640212117 0.01607142857142918 0.008465608465618502 0.07771164021165333 0.2571428571429166; - 0.0 -0.00237819664902998 -0.0016313932980599258 -0.002157738095238018 -0.0014109347442683717 -0.003789131393296341 0.048809523809520694] + theta = [0 19087oneType/362880 1139oneType/22680 137oneType/2688 143oneType/2835 3715oneType/72576 41oneType/840; + 0 2713oneType/15120 47oneType/189 27oneType/112 232oneType/945 725oneType/3024 9oneType/35; + 0 -15487oneType/120960 11oneType/7560 387oneType/4480 64oneType/945 2125oneType/24192 9oneType/280; + 0 293oneType/2835 166oneType/2835 17oneType/105 752oneType/2835 125oneType/567 34oneType/105; + 0 -6737oneType/120960 -269oneType/7560 -243oneType/4480 29oneType/945 3875oneType/24192 9oneType/280; + 0 263oneType/15120 11oneType/945 9oneType/560 8oneType/945 235oneType/3024 9oneType/35; + 0 -863oneType/362880 -37oneType/22680 -29oneType/13440 -4oneType/2835 -275oneType/72576 41oneType/840] elseif alg.M == 7 - theta = [0.0 0.04346064814814822 0.0418367346938779 0.04225127551020513 0.04202569916855836 0.042251275510214015 0.041836734693902144 0.04346064814816941; - 0.0 0.1651655801209374 0.22161753590324995 0.21667729591836782 0.21889644746787695 0.21686626039310397 0.2204081632653896 0.20700231481467313; - 0.0 -0.14384566326530607 -0.02414965986394635 0.04390943877551834 0.032653061224550584 0.04118835034029189 0.02755102040850943 0.0765625000003638; - 0.0 0.1454235166288738 0.09251700680271657 0.18899872448978527 0.26969009826137835 0.24580144557796757 0.2775510204078273 0.17297453703622523; - 0.0 -0.10457648337112646 -0.07282690854119345 -0.09671556122447367 -0.016024187452693184 0.08045753023436752 0.02755102040857693 0.17297453703748644; - 0.0 0.04901147959183677 0.03537414965986363 0.04390943877552189 0.03265306122450795 0.1007121598638605 0.22040816326534696 0.07656249999990905; - 0.0 -0.013405848450491272 -0.009863945578231281 -0.011894132653058165 -0.009674981103547253 -0.014615221088391195 0.04183673469395899 0.2070023148147584; - 0.0 0.001623913454270594 0.0012093726379440242 0.0014349489795916492 0.0012093726379429626 0.0016239134542663791 -1.0658141036401503e-14 0.04346064814814099] + theta = [0 751oneType/17280 41oneType/980 265oneType/6272 278oneType/6615 265oneType/6272 41oneType/980 751oneType/17280; + 0 139849oneType/846720 1466oneType/6615 1359oneType/6272 1448oneType/6615 36725oneType/169344 54oneType/245 3577oneType/17280; + 0 -4511oneType/31360 -71oneType/2940 1377oneType/31360 8oneType/245 775oneType/18816 27oneType/980 49oneType/640; + 0 123133oneType/846720 68oneType/735 5927oneType/31360 1784oneType/6615 4625oneType/18816 68oneType/245 2989oneType/17280; + 0 -88547oneType/846720 -1927oneType/26460 -3033oneType/31360 -106oneType/6615 13625oneType/169344 27oneType/980 2989oneType/17280; + 0 1537oneType/31360 26oneType/735 1377oneType/31360 8oneType/245 1895oneType/18816 54oneType/245 49oneType/640; + 0 -11351oneType/846720 -29oneType/2940 -373oneType/31360 -64oneType/6615 -275oneType/18816 41oneType/980 3577oneType/17280; + 0 275oneType/169344 8oneType/6615 9oneType/6272 8oneType/6615 275oneType/169344 0 751oneType/17280] elseif alg.M == 8 - theta = [0.0 0.03685850005511468 0.035688932980599594 0.03594029017857134 0.0358289241622568 0.035914937444895045 0.035803571428594694 0.03605492862659787 0.034885361552085214; - 0.0 0.15387641920194012 0.2012610229276906 0.19782924107142996 0.1990828924162642 0.19819740685622378 0.1992857142857929 0.1969121334875581 0.2076895943564523; - 0.0 -0.15861283344356245 -0.046840828924162636 0.009592633928562577 0.0021516754850381403 0.0065018050044045594 0.0016071428572104196 0.011744309413188603 -0.03273368606460281; - 0.0 0.19274133322310427 0.13237213403880332 0.22303013392857451 0.2888183421517425 0.2741522679677928 0.2878571428567511 0.2618484760821502 0.3702292768929283; - 0.0 -0.17337411816578485 -0.12799823633157104 -0.15669642857142208 -0.08007054673731773 -0.003444664903213379 -0.0321428571432989 0.013233024690180173 -0.16014109348088823; - 0.0 0.10838080081569673 0.08237213403880084 0.09607700892858873 0.08141093474421268 0.14719914296711067 0.23785714285577342 0.1774879436707124 0.37022927689395146; - 0.0 -0.044477995480599525 -0.03434082892416224 -0.03923549107142321 -0.03488536155200972 -0.04232631999563363 0.01410714285702852 0.1258791473759011 -0.03273368606702931; - 0.0 0.010777460868606693 0.008403880070546516 0.009492187499996696 0.008606701940021111 0.009860353284821599 0.006428571428614305 0.053813175154459714 0.2076895943564523; - 0.0 -0.0011695670745149895 -0.0009182098765432678 -0.0010295758928574872 -0.0009435626102302086 -0.0010549286265453262 -0.0008035714285732354 -0.0019731385031036552 0.03488536155197153] + theta = [0 1070017oneType/29030400 32377oneType/907200 12881oneType/358400 4063oneType/113400 41705oneType/1161216 401oneType/11200 149527oneType/4147200 989oneType/28350; + 0 2233547oneType/14515200 22823oneType/113400 35451oneType/179200 2822oneType/14175 115075oneType/580608 279oneType/1400 408317oneType/2073600 2944oneType/14175; + 0 -2302297oneType/14515200 -21247oneType/453600 1719oneType/179200 61oneType/28350 3775oneType/580608 9oneType/5600 24353oneType/2073600 -464oneType/14175; + 0 2797679oneType/14515200 15011oneType/113400 39967oneType/179200 4094oneType/14175 159175oneType/580608 403oneType/1400 542969oneType/2073600 5248oneType/14175; + 0 -31457oneType/181440 -2903oneType/22680 -351oneType/2240 -227oneType/2835 -125oneType/36288 -9oneType/280 343oneType/25920 -454oneType/2835; + 0 1573169oneType/14515200 9341oneType/113400 17217oneType/179200 1154oneType/14175 85465oneType/580608 333oneType/1400 368039oneType/2073600 5248oneType/14175; + 0 -645607oneType/14515200 -15577oneType/453600 -7031oneType/179200 -989oneType/28350 -24575oneType/580608 79oneType/5600 261023oneType/2073600 -464oneType/14175; + 0 156437oneType/14515200 953oneType/113400 243oneType/25600 122oneType/14175 5725oneType/580608 9oneType/1400 111587oneType/2073600 2944oneType/14175; + 0 -33953oneType/29030400 -119oneType/129600 -369oneType/358400 -107oneType/113400 -175oneType/165888 -9oneType/11200 -8183oneType/4147200 989oneType/28350] elseif alg.M == 9 - theta = [0.0 0.03188616071428564 0.031009210268469034 0.03117187499999874 0.031111111111105316 0.031149339236719698 0.031111111111073342 0.03117187499992724 0.031009210268393872 0.031886160714066136; - 0.0 0.14467159330295903 0.18532725847540765 0.1828236607142859 0.18359396433471176 0.1831509191895293 0.18357142857178133 0.18292556155574857 0.18461297276007826 0.17568080357159488; - 0.0 -0.1725594013325496 -0.06735057809132261 -0.0193750000000108 -0.02461297276119012 -0.022122403488083364 -0.02428571428754367 -0.021130829907633597 -0.02909660984732909 0.012053571401338559; - 0.0 0.24498946698020774 0.1776641191456143 0.26335317460326735 0.3186204193613946 0.30879507152621954 0.31587301587512684 0.3064180384141082 0.3290926905915512 0.21589285713457684; - 0.0 -0.2646060834313152 -0.20377621007249935 -0.23694196428579062 -0.16401332549457948 -0.10071817435664343 -0.11857142856752034 -0.098733067548892 -0.14234763861895772 0.06448660725436639; - 0.0 0.20683424578679332 0.1632196747011676 0.18305803571432477 0.1652047815015294 0.22849993263764645 0.3014285714274365 0.26826281722094336 0.329092690586549 0.06448660718228894; - 0.0 -0.11319983343131501 -0.09052518126592918 -0.09998015873019028 -0.09290221438379831 -0.10272756221900181 -0.04746031746481094 0.03822873799435911 -0.029096609840053134 0.21589285708614625; - 0.0 0.04115018126592196 0.03318440133255063 0.036339285714291236 0.034175974916735186 0.03666654418941562 0.03142857142916 0.07940414951872299 0.18461297276189725 0.0120535714286234; - 0.0 -0.008932169189692364 -0.007244757985499331 -0.007890625000002531 -0.007470115618263051 -0.00791316076327453 -0.007142857142902415 -0.009646454903545987 0.03100921026901915 0.17568080356954852; - 0.0 0.0008769504458161903 0.0007142857142857367 0.000775049603174427 0.0007368214775622661 0.0007750496031686538 0.0007142857142810044 0.0008769504457859512 5.684341886080802e-14 0.03188616071417982] + theta = [0 2857oneType/89600 3956oneType/127575 399oneType/12800 7oneType/225 81385oneType/2612736 7oneType/225 399oneType/12800 3956oneType/127575 2857oneType/89600; + 0 9449717oneType/65318400 37829oneType/204120 16381oneType/89600 3346oneType/18225 478525oneType/2612736 257oneType/1400 341383oneType/1866240 23552oneType/127575 15741oneType/89600; + 0 -1408913oneType/8164800 -34369oneType/510300 -31oneType/1600 -628oneType/25515 -7225oneType/326592 -17oneType/700 -24647oneType/1166400 -3712oneType/127575 27oneType/2240; + 0 200029oneType/816480 45331oneType/255150 13273oneType/50400 40648oneType/127575 50425oneType/163296 199oneType/630 178703oneType/583200 41984oneType/127575 1209oneType/5600; + 0 -8641823oneType/32659200 -103987oneType/510300 -2123oneType/8960 -20924oneType/127575 -131575oneType/1306368 -83oneType/700 -460649oneType/4665600 -3632oneType/25515 2889oneType/44800; + 0 6755041oneType/32659200 83291oneType/510300 8201oneType/44800 21076oneType/127575 298505oneType/1306368 211oneType/700 1251607oneType/4665600 41984oneType/127575 2889oneType/44800; + 0 -462127oneType/4082400 -9239oneType/102060 -5039oneType/50400 -11852oneType/127575 -16775oneType/163296 -299oneType/6300 4459oneType/116640 -3712oneType/127575 1209oneType/5600; + 0 335983oneType/8164800 8467oneType/255150 407oneType/11200 872oneType/25515 11975oneType/326592 11oneType/350 92617oneType/1166400 23552oneType/127575 27oneType/2240; + 0 -116687oneType/13063680 -3697oneType/510300 -101oneType/12800 -953oneType/127575 -20675oneType/2612736 -1oneType/140 -90013oneType/9331200 3956oneType/127575 15741oneType/89600; + 0 8183oneType/9331200 1oneType/1400 25oneType/32256 94oneType/127575 25oneType/32256 1oneType/1400 8183oneType/9331200 0 2857oneType/89600] + else error("MPDeC requires 2 ≤ K ≤ 10.") end else # alg.nodes == :gausslobatto if alg.M == 1 - nodes = [0.0, 1.0] - theta = [0.0 0.5; 0.0 0.5] + nodes = [0, oneType] + theta = [0 oneType/2; 0 oneType/2] elseif alg.M == 2 - nodes = [0.0, 0.5, 1.0] - theta = [0.0 0.20833333333333337 0.16666666666666652; - 0.0 0.33333333333333337 0.6666666666666667; - 0.0 -0.04166666666666667 0.16666666666666663] + nodes = [0, oneType / 2, oneType] + theta = [0 5oneType/24 oneType/6; + 0 oneType/3 2oneType/3; + 0 -oneType/24 oneType/6] elseif alg.M == 3 - nodes = [0.0, 0.27639320225002106, 0.7236067977499789, 1.0] - theta = [0.0 0.11030056647916493 0.07303276685416865 0.08333333333333393; - 0.0 0.1896994335208352 0.45057403089581083 0.41666666666666785; - 0.0 -0.033907364229143935 0.22696723314583123 0.4166666666666661; - 0.0 0.010300566479164913 -0.026967233145831604 0.08333333333333326] + nodes = [ + 0, + (1 - sqrt(oneType / 5)) / 2, + (1 + sqrt(oneType / 5)) / 2, + oneType + ] + theta = [0 (11 + sqrt(5oneType))/120 (11 - sqrt(5oneType))/120 oneType/12; + 0 (25 - sqrt(5oneType))/120 (13 * sqrt(5oneType) + 25)/120 5oneType/12; + 0 (25 - 13 * sqrt(5oneType))/120 (sqrt(5oneType) + 25)/120 5oneType/12; + 0 (sqrt(5oneType) - 1)/120 (-sqrt(5oneType) - 1)/120 oneType/12] elseif alg.M == 4 - nodes = [0.0, 0.1726731646460114, 0.5, 0.8273268353539887, 1.0] - theta = [0.0 0.0677284321861569 0.04062499999999991 0.05370013924241501 0.05000000000000071; - 0.0 0.11974476934341162 0.30318418332304287 0.2615863979968083 0.27222222222222214; - 0.0 -0.021735721866558116 0.17777777777777748 0.3772912774221129 0.35555555555555296; - 0.0 0.010635824225415487 -0.0309619611008205 0.1524774528788102 0.27222222222222037; - 0.0 -0.0037001392424145354 0.009375000000000022 -0.017728432186157494 0.04999999999999982] + nodes = [ + 0, + (1 - sqrt(3oneType / 7)) / 2, + oneType / 2, + (1 + sqrt(3oneType / 7)) / 2, + oneType + ] + + theta = [0 (3sqrt(21oneType) + 119)/1960 13oneType/320 (119 - 3sqrt(21oneType))/1960 oneType/20; + 0 49oneType / 360-sqrt(21oneType) / 280 (7sqrt(21oneType)) / 192+49oneType / 360 (23sqrt(21oneType)) / 840+49oneType / 360 49oneType/180; + 0 8oneType / 45-(32sqrt(21oneType)) / 735 8oneType/45 (32sqrt(21oneType)) / 735+8oneType / 45 16oneType/45; + 0 49oneType / 360-(23sqrt(21oneType)) / 840 49oneType / 360-(7sqrt(21oneType)) / 192 sqrt(21oneType) / 280+49oneType / 360 49oneType/180; + 0 ((3sqrt(21oneType)) / 1960-3oneType / 280) 3oneType/320 (-(3sqrt(21oneType)) / 1960-3oneType / 280) 1oneType/20] + elseif alg.M == 5 nodes = [ - 0.0, - 0.11747233803526769, - 0.3573842417596774, - 0.6426157582403226, - 0.8825276619647323, - 1.0 + 0, + oneType / 2 - sqrt((2sqrt(7oneType)) / 21 + oneType / 3) / 2, + oneType / 2 - sqrt(oneType / 3 - (2sqrt(7oneType)) / 21) / 2, + oneType / 2 + sqrt(oneType / 3 - (2sqrt(7oneType)) / 21) / 2, + oneType / 2 + sqrt((2sqrt(7oneType)) / 21 + oneType / 3) / 2, + oneType ] - theta = [0.0 0.04567980513375505 0.025908385387879762 0.03746264288972734 0.03168990732937349 0.033333333333333215; - 0.0 0.08186781700897068 0.2138408086328255 0.177429781771262 0.19370925858950017 0.18923747814892522; - 0.0 -0.01487460578908985 0.13396073565086075 0.30143326325089315 0.2698015123994857 0.2774291885177327; - 0.0 0.007627676118250971 -0.024004074733154912 0.14346845286688126 0.2923037943068376 0.27742918851774334; - 0.0 -0.004471780440573713 0.011807696377659743 -0.02460333048390262 0.10736966113994595 0.18923747814891745; - 0.0 0.0016434260039545345 -0.004129309556393734 0.007424947945453564 -0.012346471800418257 0.03333333333333588] + + theta = [0 (sqrt(7oneType) / 756 + sqrt(6sqrt(7oneType) + 21) / 540 - sqrt(42sqrt(7oneType) + 147) / 3780+19oneType / 540) (19oneType / 540 - sqrt(21 - 6sqrt(7oneType)) / 540 - sqrt(147 - 42sqrt(7oneType)) / 3780-sqrt(7oneType) / 756) (sqrt(21 - 6sqrt(7oneType)) / 540-sqrt(7oneType) / 756+sqrt(147 - 42sqrt(7oneType))/3780+19oneType/540) (sqrt(7oneType) / 756-sqrt(6sqrt(7oneType) + 21) / 540+sqrt(42sqrt(7oneType) + 147)/3780+19oneType/540) oneType/30; + 0 (7oneType / 60 - sqrt(42sqrt(7oneType) + 147) / 1260-sqrt(7oneType) / 120) (-((sqrt(3oneType) - 2sqrt(21oneType)) * (63sqrt(21oneType) + 35sqrt(2sqrt(7oneType) + 7) + 70sqrt(6sqrt(7oneType) + 21) + 20sqrt(14sqrt(7oneType) + 49) - 23sqrt(42sqrt(7oneType) + 147)))/22680) (-((sqrt(3oneType) - 2sqrt(21oneType)) * (63sqrt(21oneType) + 35sqrt(2sqrt(7oneType) + 7) - 70sqrt(6sqrt(7oneType) + 21) + 20sqrt(14sqrt(7oneType) + 49) + 23sqrt(42sqrt(7oneType) + 147)))/22680) (sqrt(42sqrt(7oneType) + 147) / 60 - sqrt(6sqrt(7oneType) + 21) / 36 - sqrt(7oneType) / 120+7oneType / 60) (7oneType / 30-sqrt(7oneType) / 60); + 0 (-((sqrt(3oneType) + 2sqrt(21oneType)) * (21sqrt(2sqrt(7oneType) + 7) - 63sqrt(21oneType) + 70sqrt(6sqrt(7oneType) + 21) + 24 * sqrt(14sqrt(7oneType) + 49) - 25sqrt(42sqrt(7oneType) + 147)))/22680) (sqrt(7oneType) / 120 - sqrt(147 - 42sqrt(7oneType)) / 1260+7oneType / 60) (sqrt(7oneType)/120+sqrt(21 - 6sqrt(7oneType))/36+sqrt(147 - 42sqrt(7oneType))/60+7oneType/60) (((sqrt(3oneType) + 2sqrt(21oneType)) * (63sqrt(21oneType) + 21sqrt(2sqrt(7oneType) + 7) - 70sqrt(6sqrt(7oneType) + 21) + 24 * sqrt(14sqrt(7oneType) + 49) + 25sqrt(42sqrt(7oneType) + 147)))/22680) (sqrt(7oneType) / 60+7oneType / 30); + 0 (-((sqrt(3oneType) + 2sqrt(21oneType)) * (21sqrt(2sqrt(7oneType) + 7) - 63sqrt(21oneType) - 70sqrt(6sqrt(7oneType) + 21) + 24 * sqrt(14sqrt(7oneType) + 49) + 25sqrt(42sqrt(7oneType) + 147)))/22680) (sqrt(7oneType) / 120 - sqrt(21 - 6sqrt(7oneType)) / 36 - sqrt(147 - 42sqrt(7oneType)) / 60+7oneType / 60) (((sqrt(7oneType) + 14) * (32 * sqrt(2sqrt(7oneType) + 7) - 10sqrt(14sqrt(7oneType) + 49) + 567))/68040) (((sqrt(3oneType) + 2sqrt(21oneType)) * (63sqrt(21oneType) + 21sqrt(2sqrt(7oneType) + 7) + 70sqrt(6sqrt(7oneType) + 21) + 24 * sqrt(14sqrt(7oneType) + 49) - 25sqrt(42sqrt(7oneType) + 147)))/22680) (sqrt(7oneType) / 60+7oneType / 30); + 0 (sqrt(6sqrt(7oneType) + 21) / 36 - sqrt(7oneType) / 120 - sqrt(42sqrt(7oneType) + 147) / 60+7oneType / 60) (((sqrt(3oneType) - 2sqrt(21oneType)) * (35sqrt(2sqrt(7oneType) + 7) - 63sqrt(21oneType) - 70sqrt(6sqrt(7oneType) + 21) + 20sqrt(14sqrt(7oneType) + 49) + 23sqrt(42sqrt(7oneType) + 147)))/22680) (((sqrt(3oneType) - 2sqrt(21oneType)) * (35sqrt(2sqrt(7oneType) + 7) - 63sqrt(21oneType) + 70sqrt(6sqrt(7oneType) + 21) + 20sqrt(14sqrt(7oneType) + 49) - 23sqrt(42sqrt(7oneType) + 147)))/22680) (sqrt(42sqrt(7oneType) + 147) / 1260 - sqrt(7oneType) / 120+7oneType / 60) (7oneType / 30-sqrt(7oneType) / 60); + 0 (sqrt(6sqrt(7oneType) + 21) / 540 - sqrt(7oneType) / 756 - sqrt(42sqrt(7oneType) + 147) / 3780-oneType / 540) (sqrt(7oneType) / 756 - sqrt(21 - 6sqrt(7oneType)) / 540 - sqrt(147 - 42sqrt(7oneType)) / 3780-oneType / 540) (sqrt(7oneType) / 756 + sqrt(21 - 6sqrt(7oneType)) / 540 + sqrt(147 - 42sqrt(7oneType)) / 3780-oneType / 540) (sqrt(42sqrt(7oneType) + 147) / 3780 - sqrt(6sqrt(7oneType) + 21) / 540 - sqrt(7oneType) / 756-oneType / 540) oneType/30] else error("MPDeC requires 2 ≤ K ≤ 10.") end @@ -216,7 +232,7 @@ function alg_cache(alg::MPDeC, u, rate_prototype, ::Type{uEltypeNoUnits}, throw(ArgumentError("MPDeC can only be applied to production-destruction systems")) end - nodes, theta = get_constant_parameters(alg) + nodes, theta = get_constant_parameters(alg, uEltypeNoUnits) MPDeCConstantCache(alg.K, alg.M, nodes, theta, alg.small_constant_function(uEltypeNoUnits)) end @@ -610,7 +626,7 @@ function alg_cache(alg::MPDeC, u, rate_prototype, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck, ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} - nodes, theta = get_constant_parameters(alg) + nodes, theta = get_constant_parameters(alg, uEltypeNoUnits) tab = MPDeCConstantCache(alg.K, alg.M, nodes, theta, alg.small_constant_function(uEltypeNoUnits)) From f07407aaad33204870eb18ceea41655ddef11837 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Mon, 24 Mar 2025 07:40:20 +0100 Subject: [PATCH 29/38] fix strato benchmark --- docs/src/stratospheric_reaction_benchmark.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/stratospheric_reaction_benchmark.md b/docs/src/stratospheric_reaction_benchmark.md index cb6c3a60..4b746a91 100644 --- a/docs/src/stratospheric_reaction_benchmark.md +++ b/docs/src/stratospheric_reaction_benchmark.md @@ -23,7 +23,7 @@ function stratreac_plot(sols, labels = fill("", length(sols)), sol_ref = nothing labels = [labels] end - tspan = prob.tspan + tspan = prob_pds_stratreac.tspan layout = (3, 2) linewidth = 2 xticks = (range(first(tspan), last(tspan), 4), range(12.0, 84.0, 4)) From c71240284dbfdb16243d616493d9ca4ff2c18862 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Thu, 27 Mar 2025 14:36:43 +0100 Subject: [PATCH 30/38] Convergence tests for high order schemes --- src/mpdec.jl | 8 +-- test/runtests.jl | 177 +++++++++++++++++++++++++++++++++++++++++++++-- 2 files changed, 177 insertions(+), 8 deletions(-) diff --git a/src/mpdec.jl b/src/mpdec.jl index 8962f399..ccd4398b 100644 --- a/src/mpdec.jl +++ b/src/mpdec.jl @@ -559,11 +559,11 @@ end N = length(uprev) if uprev isa StaticArray - C = MMatrix{N, M + 1}(zeros(N, M + 1)) - C2 = MMatrix{N, M + 1}(zeros(N, M + 1)) + C = MMatrix{N, M + 1}(zeros(eltype(uprev), N, M + 1)) + C2 = MMatrix{N, M + 1}(zeros(eltype(uprev), N, M + 1)) else - C = zeros(N, M + 1) - C2 = zeros(N, M + 1) + C = zeros(eltype(uprev), N, M + 1) + C2 = zeros(eltype(uprev), N, M + 1) end for i in 1:(M + 1) diff --git a/test/runtests.jl b/test/runtests.jl index c180640a..8f958c7c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1911,7 +1911,6 @@ end end end - #TODO: Check order of MPDeC(K) for K ≥ 4 # Here we check the convergence order of pth-order schemes for which # no interpolation of order p is available @testset "Convergence tests (conservative)" begin @@ -1929,7 +1928,6 @@ end # Here we check the convergence order of pth-order schemes for which # no interpolation of order p is available - #TODO: Check order of MPDeC(K), K≥ 5 @testset "Convergence tests (nonconservative)" begin dts = 0.5 .^ (4:12) problems = (prob_pds_linmod_nonconservative, @@ -1968,7 +1966,6 @@ end end end - #TODO: Check MPDeC(K) order for K ≥ 4 @testset "Check convergence order (nonautonomous conservative PDS)" begin prod! = (P, u, p, t) -> begin fill!(P, zero(eltype(P))) @@ -1999,7 +1996,6 @@ end end end - #TODO: Check order of MPDeC(K), K≥ 6 @testset "Check convergence order (nonautonomous nonconservative PDS)" begin prod! = (P, u, p, t) -> begin fill!(P, zero(eltype(P))) @@ -2044,6 +2040,179 @@ end end end + @testset "Convergence of higher order schemes (autonomous)" begin + function f_analytic(u0, p, t) + u₁⁰, u₂⁰ = u0 + a, b = p + c = a + b + return ((u₁⁰ + u₂⁰) * [b; a] + + exp(-c * t) * (a * u₁⁰ - b * u₂⁰) * [1; -1]) / c + end + + # linear model problem - conservative + function linmodP!(P, u, p, t) + P[1, 2] = u[2] + P[2, 1] = 5u[1] + return nothing + end + function linmodP(u, p, t) + P = zeros(eltype(u), 2, 2) + linmodP!(P, u, p, t) + return SMatrix{2, 2}(P) + end + + # linear model problem - nonconservative - in-place + function linmodP_noncons!(P, u, p, t) + P[1, 1] = u[2] / 2 + P[1, 2] = u[2] / 2 + P[2, 1] = 3u[1] + P[2, 2] = 2u[1] + return nothing + end + function linmodD_noncons!(D, u, p, t) + D[1] = 2u[1] + D[2] = u[2] / 2 + return nothing + end + function linmodP_noncons(u, p, t) + P = zeros(eltype(u), 2, 2) + linmodP_noncons!(P, u, p, t) + return SMatrix{2, 2}(P) + end + function linmodD_noncons(u, p, t) + D = zeros(eltype(u), 2, 1) + linmodD_noncons!(D, u, p, t) + return SVector{2}(D) + end + + u0 = [Double64(9) / 10; Double64(1) / 10] + p = [Double64(5); Double64(1)] + tspan = (Double64(0), Double64(1)) + + prob_ip = ConservativePDSProblem(linmodP!, u0, tspan, p; analytic = f_analytic) + prob_op = ConservativePDSProblem(linmodP, SVector{2}(u0), tspan, SVector{2}(p); + analytic = f_analytic) + prob_ip_noncons = PDSProblem(linmodP_noncons!, linmodD_noncons!, u0, tspan, p; + analytic = f_analytic) + prob_op_noncons = PDSProblem(linmodP_noncons, linmodD_noncons, SVector{2}(u0), + tspan, + SVector{2}(p); analytic = f_analytic) + + dts = 0.5 .^ (8:13) + + algs = [] + for K in 4:10 + push!(algs, MPDeC(K)) + push!(algs, MPDeC(K, nodes = :lagrange)) + end + + atol = 0.3 + @testset "$alg" for alg in algs + orders = experimental_orders_of_convergence(prob_ip, alg, dts) + @test check_order(orders, PositiveIntegrators.alg_order(alg); atol) + + orders = experimental_orders_of_convergence(prob_op, alg, dts) + @test check_order(orders, PositiveIntegrators.alg_order(alg); atol) + + orders = experimental_orders_of_convergence(prob_ip_noncons, alg, dts) + @test check_order(orders, PositiveIntegrators.alg_order(alg); atol) + + orders = experimental_orders_of_convergence(prob_op_noncons, alg, dts) + @test check_order(orders, PositiveIntegrators.alg_order(alg); atol) + end + end + + @testset "Convergence of higher order schemes (nonautonomous)" begin + prod! = (P, u, p, t) -> begin + fill!(P, zero(eltype(P))) + P[1, 2] = sin(t)^2 * u[2] + P[2, 1] = cos(t)^2 * u[1] + return nothing + end + prod = (u, p, t) -> begin + P = similar(u, (length(u), length(u))) + prod!(P, u, p, t) + return SMatrix{2, 2}(P) + end + function fct!(du, u, p, t) + du[1] = sin(t)^2 * u[2] - cos(t)^2 * u[1] + du[2] = -sin(t)^2 * u[2] + cos(t)^2 * u[1] + return nothing + end + fct = (u, p, t) -> begin + f = similar(u) + fct!(f, u, p, t) + return SVector{2}(f) + end + + prod_noncons! = (P, u, p, t) -> begin + fill!(P, zero(eltype(P))) + P[1, 2] = sin(t)^2 * u[2] + P[2, 1] = cos(2 * t)^2 * u[1] + P[2, 2] = cos(t)^2 * u[2] + return nothing + end + dest_noncons! = (d, u, p, t) -> begin + fill!(d, zero(eltype(d))) + d[1] = sin(2 * t)^2 * u[1] + d[2] = sin(0.5 * t)^2 * u[2] + return nothing + end + prod_noncons = (u, p, t) -> begin + P = similar(u, (length(u), length(u))) + prod_noncons!(P, u, p, t) + return SMatrix{2, 2}(P) + end + dest_noncons = (u, p, t) -> begin + d = similar(u, (length(u),)) + dest_noncons!(d, u, p, t) + return SVector{2}(d) + end + fct_noncons! = (du, u, p, t) -> begin + du[1] = sin(t)^2 * u[2] - cos(2 * t)^2 * u[1] - sin(2 * t)^2 * u[1] + du[2] = -sin(t)^2 * u[2] + cos(2 * t)^2 * u[1] + cos(t)^2 * u[2] - + sin(0.5 * t)^2 * u[2] + return nothing + end + fct_noncons = (u, p, t) -> begin + f = similar(u) + fct_noncons!(f, u, p, t) + return SVector{2}(f) + end + + u0 = [Double64(11) / 10; Double64(9) / 10] # values close to zero may decrease the order + tspan = (Double64(0), Double64(1)) + prob_op = ConservativePDSProblem(prod, SVector{2}(u0), tspan; std_rhs = fct) #out-of-place + prob_ip = ConservativePDSProblem(prod!, u0, tspan; std_rhs = fct!) #in-place + prob_op_noncons = PDSProblem(prod_noncons, dest_noncons, SVector{2}(u0), tspan; + std_rhs = fct_noncons) #out-of-place + prob_ip_noncons = PDSProblem(prod_noncons!, dest_noncons!, u0, tspan; + std_rhs = fct_noncons!) #in-place + + dts = 0.5 .^ (5:10) + + algs = [] + for K in 4:10 + push!(algs, MPDeC(K)) + push!(algs, MPDeC(K, nodes = :lagrange)) + end + + atol = 0.3 + @testset "$alg" for alg in algs + orders = experimental_orders_of_convergence(prob_ip, alg, dts) + @test check_order(orders, PositiveIntegrators.alg_order(alg); atol) + + orders = experimental_orders_of_convergence(prob_op, alg, dts) + @test check_order(orders, PositiveIntegrators.alg_order(alg); atol) + + orders = experimental_orders_of_convergence(prob_ip_noncons, alg, dts) + @test check_order(orders, PositiveIntegrators.alg_order(alg); atol) + + orders = experimental_orders_of_convergence(prob_op_noncons, alg, dts) + @test check_order(orders, PositiveIntegrators.alg_order(alg); atol) + end + end + @testset "Interpolation tests (nonconservative)" begin algs = (MPE(), MPRK22(0.5), MPRK22(1.0), MPRK22(2.0), MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), MPRK43II(0.5), MPRK43II(2.0 / 3.0), From 736780a575cf007605f72a93c8a4a527def6ceb7 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Thu, 27 Mar 2025 14:41:12 +0100 Subject: [PATCH 31/38] Updated NEWS.md --- NEWS.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/NEWS.md b/NEWS.md index 00c80934..d0b4dc63 100644 --- a/NEWS.md +++ b/NEWS.md @@ -12,6 +12,9 @@ for human readability. - The minimum required Julia version was updated to v1.10 +### Added + +- MPDeC methods of Öffner and Torlo ## Changes when updating to v0.2 from v0.1.x From c63e353e0aa57d7fb1b4578a6ba72b2c704ee250 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Thu, 27 Mar 2025 15:56:43 +0100 Subject: [PATCH 32/38] Added DoubleFloats to test environment --- test/Project.toml | 2 ++ test/runtests.jl | 1 + 2 files changed, 3 insertions(+) diff --git a/test/Project.toml b/test/Project.toml index 464f1e9f..edcab80a 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,6 +1,7 @@ [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" +DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78" ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" @@ -20,6 +21,7 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" [compat] ADTypes = "1.12" Aqua = "0.8" +DoubleFloats = "1.4.3" ExplicitImports = "1.0.1" LinearAlgebra = "1" LinearSolve = "2.32" diff --git a/test/runtests.jl b/test/runtests.jl index 6907da2c..e5a789c4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,6 +3,7 @@ using LinearAlgebra using SparseArrays using Statistics: mean, median +using DoubleFloats using StaticArrays: MVector, @SVector, SA using Unitful: @u_str, ustrip From d078fd76c1e892c1af67c4c87c1f2af63aef7cbb Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Thu, 27 Mar 2025 20:23:08 +0100 Subject: [PATCH 33/38] made SVector available in tests --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index e5a789c4..cb737161 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,7 +4,7 @@ using SparseArrays using Statistics: mean, median using DoubleFloats -using StaticArrays: MVector, @SVector, SA +using StaticArrays: MVector, @SVector, SVector, SA using Unitful: @u_str, ustrip From 8a2eaa2e7b1fd11b96a0e3e13ac205d3934410fd Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Thu, 27 Mar 2025 21:14:17 +0100 Subject: [PATCH 34/38] made SMatrix available in tests --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index cb737161..f01aee15 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,7 +4,7 @@ using SparseArrays using Statistics: mean, median using DoubleFloats -using StaticArrays: MVector, @SVector, SVector, SA +using StaticArrays: SMatrix, MVector, @SVector, SVector, SA using Unitful: @u_str, ustrip From b99b3966fb17a838724baaf06412328fd7902c70 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Thu, 27 Mar 2025 21:24:34 +0100 Subject: [PATCH 35/38] import only Double64 from DoubleFloats --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index f01aee15..6cf2aef2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,7 +3,7 @@ using LinearAlgebra using SparseArrays using Statistics: mean, median -using DoubleFloats +using DoubleFloats: Double64 using StaticArrays: SMatrix, MVector, @SVector, SVector, SA using Unitful: @u_str, ustrip From 82e91686137abf93fce7aaa361e70e472deeb15e Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Sun, 30 Mar 2025 16:32:12 +0200 Subject: [PATCH 36/38] added suggested changes from reviews --- src/mpdec.jl | 18 +++++++++--------- test/runtests.jl | 2 +- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/src/mpdec.jl b/src/mpdec.jl index ccd4398b..7bd655f4 100644 --- a/src/mpdec.jl +++ b/src/mpdec.jl @@ -52,7 +52,7 @@ end function small_constant_function_MPDeC(type) if type == Float64 # small_constant is chosen such that - # the testet "Zero initial values" passes. + # the testset "Zero initial values" passes. small_constant = 1e-300 else small_constant = floatmin(type) @@ -72,8 +72,10 @@ function MPDeC(K::Integer; if nodes === :lagrange M = K - 1 - else # :gausslobatto + elseif nodes === :gausslobatto M = cld(K, 2) + else + throw(ArgumentError("Unknown node choice for MPDeC. Possible options are :gauslobatto or :lagrange")) end MPDeC{typeof(nodes), typeof(linsolve), typeof(small_constant_function)}(K, M, @@ -87,7 +89,7 @@ isfsal(::MPDeC) = false function get_constant_parameters(alg::MPDeC, type) oneType = one(type) - if alg.nodes == :lagrange + if alg.nodes === :lagrange nodes = collect(i * oneType / alg.M for i in 0:(alg.M)) # M is one less than the methods order if alg.M == 1 @@ -157,7 +159,7 @@ function get_constant_parameters(alg::MPDeC, type) else error("MPDeC requires 2 ≤ K ≤ 10.") end - else # alg.nodes == :gausslobatto + else # alg.nodes === :gausslobatto if alg.M == 1 nodes = [0, oneType] theta = [0 oneType/2; 0 oneType/2] @@ -267,7 +269,6 @@ end @muladd function _build_mpdec_matrix_and_rhs_oop(uprev, m, prod, C, p, t, dt, nodes, theta, small_constant, dest = nothing) N, M = size(C) - M = M - 1 # Create linear system matrix and rhs if uprev isa StaticArray @@ -286,7 +287,7 @@ end σ = add_small_constant(C[:, m], small_constant) - @fastmath @inbounds @simd for r in 1:(M + 1) + @fastmath @inbounds @simd for r in 1:M th = theta[r, m] dt_th = dt * th P = prod(C[:, r], p, t + nodes[r] * dt) @@ -306,12 +307,11 @@ end nodes, theta, small_constant, dest = nothing, d = nothing) N, M = size(C) - M = M - 1 oneMmat = one(eltype(Mmat)) zeroMmat = zero(eltype(Mmat)) - #Initialize Mmat as identity matrix + # Initialize Mmat as identity matrix if Mmat isa Tridiagonal Mmat.d .= oneMmat Mmat.du .= zeroMmat @@ -340,7 +340,7 @@ end σ .= C[:, m] .+ small_constant - @fastmath @inbounds @simd for r in 1:(M + 1) + @fastmath @inbounds @simd for r in 1:M th = theta[r, m] dt_th = dt * th diff --git a/test/runtests.jl b/test/runtests.jl index 6cf2aef2..7b921736 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2101,7 +2101,7 @@ end dts = 0.5 .^ (8:13) - algs = [] + algs = MPDeC[] for K in 4:10 push!(algs, MPDeC(K)) push!(algs, MPDeC(K, nodes = :lagrange)) From fb7bfb89c01a68e9f01eac1ab5876c44e54fc875 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Sun, 30 Mar 2025 19:50:50 +0200 Subject: [PATCH 37/38] format --- test/runtests.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 2f05bd2a..64b6c2fa 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1741,7 +1741,6 @@ end sol_sparse_ip = solve(prob_sparse_ip, alg) sol_sparse_op = solve(prob_sparse_op, alg) - @test isapprox(sol_tridiagonal_ip.t, sol_tridiagonal_op.t; rtol) @test isapprox(sol_dense_ip.t, sol_dense_op.t; rtol) @test isapprox(sol_sparse_ip.t, sol_sparse_op.t; rtol) @@ -1824,7 +1823,6 @@ end end #solve and test for alg in algs - for prob in (prob_default, prob_tridiagonal, prob_dense, prob_sparse, prob_default2, prob_tridiagonal2, prob_dense2, prob_sparse2) From ea54eef90c1c2a7fd627af79c709e1983b6cc80d Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Sun, 30 Mar 2025 20:35:50 +0200 Subject: [PATCH 38/38] use better variable name M_plus_1 --- src/mpdec.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/mpdec.jl b/src/mpdec.jl index 7bd655f4..c31d4ee0 100644 --- a/src/mpdec.jl +++ b/src/mpdec.jl @@ -268,7 +268,7 @@ end # out-of-place for dense arrays @muladd function _build_mpdec_matrix_and_rhs_oop(uprev, m, prod, C, p, t, dt, nodes, theta, small_constant, dest = nothing) - N, M = size(C) + N, M_plus_1 = size(C) # Create linear system matrix and rhs if uprev isa StaticArray @@ -287,7 +287,7 @@ end σ = add_small_constant(C[:, m], small_constant) - @fastmath @inbounds @simd for r in 1:M + @fastmath @inbounds @simd for r in 1:M_plus_1 th = theta[r, m] dt_th = dt * th P = prod(C[:, r], p, t + nodes[r] * dt) @@ -306,7 +306,7 @@ end @muladd function build_mpdec_matrix_and_rhs_ip!(Mmat, rhs, m, prod, P, C, p, t, dt, σ, tmp, nodes, theta, small_constant, dest = nothing, d = nothing) - N, M = size(C) + N, M_plus_1 = size(C) oneMmat = one(eltype(Mmat)) zeroMmat = zero(eltype(Mmat)) @@ -340,7 +340,7 @@ end σ .= C[:, m] .+ small_constant - @fastmath @inbounds @simd for r in 1:M + @fastmath @inbounds @simd for r in 1:M_plus_1 th = theta[r, m] dt_th = dt * th