diff --git a/Project.toml b/Project.toml index 43081c23..6e217c52 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "SummationByPartsOperators" uuid = "9f78cca6-572e-554e-b819-917d2f1cf240" author = ["Hendrik Ranocha"] -version = "0.5.66" +version = "0.5.67" [deps] ArgCheck = "dce04be8-c92d-5529-be00-80e4d2c0e197" diff --git a/src/SBP_coefficients/WilliamsDuru2024.jl b/src/SBP_coefficients/WilliamsDuru2024.jl new file mode 100644 index 00000000..b8ccc521 --- /dev/null +++ b/src/SBP_coefficients/WilliamsDuru2024.jl @@ -0,0 +1,766 @@ +""" + WilliamsDuru2024(version::Symbol) + +Coefficients of the upwind SBP operators given in +- Williams, Duru (2024) + Full-Spectrum Dispersion Relation Preserving Summation-By-Parts Operators. + SIAM Journal on Numerical Analysis 62.4, pp. 1565-1588. + +You can choose between the different versions `:central`, `:plus`, and `:minus`. +""" +struct WilliamsDuru2024 <: SourceOfCoefficients + kind::Symbol + + function WilliamsDuru2024(kind::Symbol) + if (kind !== :plus) && (kind !== :minus) && (kind !== :central) + throw(ArgumentError("The only choices are :plus, :minus, and :central, not :$kind.")) + end + new(kind) + end +end + +function Base.show(io::IO, source::WilliamsDuru2024) + if get(io, :compact, false) + summary(io, source) + else + print(io, + "Williams, Duru (2024) \n", + " Full-Spectrum Dispersion Relation Preserving Summation-By-Parts Operators. \n", + " SIAM Journal on Numerical Analysis 62.4, pp. 1565-1588. \n", + " (upwind coefficients ", source.kind, ")") + end +end + + +function first_derivative_coefficients(source::WilliamsDuru2024, order::Int, + T=Float64, mode=FastMode()) + if order == 4 + left_boundary_plus = ( + DerivativeCoefficientRow{T,1, 6}(SVector( T(-7244216288712856142349841908536610983621117531423//4383417301651166629768772548289373010602184474830), + T(334203102386451445389405940943896029825001960741//146113910055038887658959084942979100353406149161), + T(-154761287365716184513371737750091076641062676551//292227820110077775317918169885958200706812298322), + T(-47551546840572027745239855308325563115118731109//438341730165116662976877254828937301060218447483), + T(-2027691859629510742263865985535929607908586//25292350710583155211867593031500623221984793), + T(10572883047809310846529020642420628641879384//126461753552915776059337965157503116109923965), )), + DerivativeCoefficientRow{T,1, 6}(SVector( T(-103891878110009340509786809946415352778485344626//214324208707778398746264664339045603364518299681), + T(-30656114770770165599996417750517210097119869413//1285945252246670392477587986034273620187109798086), + T(105617810201506085617651923044209572660186559252//214324208707778398746264664339045603364518299681), + T(1724489859413945488654948596921605417066392689//142882805805185599164176442892697068909678866454), + T(1841279993087236048336359167378795063646952//111298706270267473816651201837828771004596659), + T(-475853544315554201937773124125473316948296//37099568756755824605550400612609590334865553), )), + DerivativeCoefficientRow{T,1, 7}(SVector( T(24402651831679141629765081741317534139059529921//64879085639584867376981030570117731701177761942), + T(-115008374618384762928474480582274121580026720269//97318628459377301065471545855176597551766642913), + T(-25337437961186700206692845540531216985830930791//194637256918754602130943091710353195103533285826), + T(2952784101484936236780961439554019034921446585//2495349447676341052960808868081451219276067767), + T(-120564914274524309894860197017988707048800//1295836652765972504566803982039874270006613), + T(-247565206189901598995271496383878380330120//1295836652765972504566803982039874270006613), + T(205497436417156398146302285161289402323048//5615292161985880853122817255506121836695323), )), + DerivativeCoefficientRow{T,1, 8}(SVector( T(-2045943276599844567924465770331924307140274804//69169663884034527267421986449825308966044065547), + T(70833766889864192195958289226074523178045259801//1152827731400575454457033107497088482767401092450), + T(-283405962527858709629721527199101966603555673374//576413865700287727228516553748544241383700546225), + T(-88930074506280496946385186439556417902700092547//691696638840345272674219864498253089660440655470), + T(12899980754583921740520921771537796210010688//19955473972660125574814490349612056132376685), + T(-1031661509381888342683499028256291845530976//19955473972660125574814490349612056132376685), + T(-2484650822134709177587109447859226409905944//99777369863300627874072451748060280661883425), + T(1849476927754407583316720566451604620907432//99777369863300627874072451748060280661883425), )), + DerivativeCoefficientRow{T,1, 9}(SVector( T(-30848740023524146777443615//1274365888916385703504957136), + T(-2847973112561150368136731//477887208343644638814358926), + T(35193984920688246728384925//159295736114548212938119642), + T(-95097865163352317713472025//159295736114548212938119642), + T(-2501516650966057622106220669//3823097666749157110514871408), + T(448480699319675746249380855//318591472229096425876239284), + T(-215606731781823666495914871//637182944458192851752478568), + T(-10949100926682912426100297//318591472229096425876239284), + T(8150082644673746843488191//318591472229096425876239284), )), + DerivativeCoefficientRow{T,1,10}(SVector( T(239243576443645949444669436//15024124021976694625393539385), + T(-43208769313625564285244989//3004824804395338925078707877), + T(-65895740807931831205699290//3004824804395338925078707877), + T(-179675720290104294481700384//3004824804395338925078707877), + T(-86853447171538637656723943//3004824804395338925078707877), + T(-29791922049731222907256317667//30048248043953389250787078770), + T(4231334563340413936350007050//3004824804395338925078707877), + T(-1801929304290820178748076173//6009649608790677850157415754), + T(-91506910068987813582853811//3004824804395338925078707877), + T(68114166141577395072951333//3004824804395338925078707877), )), + ) + right_boundary_plus = ( + DerivativeCoefficientRow{T,1, 6}(SVector( T(6645540386415520249117061217571584879855787531777//4383417301651166629768772548289373010602184474830), + T(-311675634330028021529360429839246058335456033878//146113910055038887658959084942979100353406149161), + T(219623866485112274667885735671857807251535769289//292227820110077775317918169885958200706812298322), + T(-51148581914996114198111644258298107678506870100//438341730165116662976877254828937301060218447483), + T(-1750105962455727357832249604331334378691370//25292350710583155211867593031500623221984793), + T(6496085907097095583419573712929848231230944//126461753552915776059337965157503116109923965), )), + DerivativeCoefficientRow{T,1, 6}(SVector( T(334203102386451445389405940943896029825001960741//642972626123335196238793993017136810093554899043), + T(-30656114770770165599996417750517210097119869413//1285945252246670392477587986034273620187109798086), + T(-115008374618384762928474480582274121580026720269//214324208707778398746264664339045603364518299681), + T(70833766889864192195958289226074523178045259801//1285945252246670392477587986034273620187109798086), + T(-430855390649649532271220388941725427561008//111298706270267473816651201837828771004596659), + T(-391076856445741420231527017941029013349752//37099568756755824605550400612609590334865553), )), + DerivativeCoefficientRow{T,1, 6}(SVector( T(-154761287365716184513371737750091076641062676551//583911770756263806392829275131059585310599857478), + T(8124446938577391201357840234169967127706658404//7486048343029023158882426604244353657828203301), + T(-25337437961186700206692845540531216985830930791//194637256918754602130943091710353195103533285826), + T(-283405962527858709629721527199101966603555673374//291955885378131903196414637565529792655299928739), + T(1774773089750676192133205715016194910450800//5615292161985880853122817255506121836695323), + T(-15292657306572795914969129703142596942480//431945550921990834855601327346624756668871), )), + DerivativeCoefficientRow{T,1, 7}(SVector( T(-47551546840572027745239855308325563115118731109//1729241597100863181685549661245632724151101638675), + T(15520408734725509397894537372294448753597534201//1152827731400575454457033107497088482767401092450), + T(69095147974747507940674497685564045417161850089//115282773140057545445703310749708848276740109245), + T(-88930074506280496946385186439556417902700092547//691696638840345272674219864498253089660440655470), + T(-1726424775466658030439702201256805556456624//3991094794532025114962898069922411226475337), + T(-4878663543549720465771003262056834494382336//99777369863300627874072451748060280661883425), + T(2465969237005876777755627421935472827876576//99777369863300627874072451748060280661883425), )), + DerivativeCoefficientRow{T,1, 8}(SVector( T(-35741686713503424823053247//1274365888916385703504957136), + T(24341883735526123464906853//955774416687289277628717852), + T(-5180106064713036907224775//79647868057274106469059821), + T(71057867556242896847948430//79647868057274106469059821), + T(-2501516650966057622106220669//3823097666749157110514871408), + T(-10392298878790214848122661//318591472229096425876239284), + T(-109408685199711207626220261//637182944458192851752478568), + T(2716694214891248947829397//79647868057274106469059821), )), + DerivativeCoefficientRow{T,1, 9}(SVector( T(389387454207582946160914171//15024124021976694625393539385), + T(-52575461023873067818412747//3004824804395338925078707877), + T(-355584433607894620879408295//3004824804395338925078707877), + T(-189974696081312018416760220//3004824804395338925078707877), + T(3748169214543475593647140365//3004824804395338925078707877), + T(-29791922049731222907256317667//30048248043953389250787078770), + T(117651741517270046035097757//3004824804395338925078707877), + T(-914381078809660182342952743//6009649608790677850157415754), + T(90818888188769860097268444//3004824804395338925078707877), )), + ) + upper_coef_plus = SVector(T(205//143), + T(-873//2860), + T(-133//4290), + T(3//130), ) + central_coef_plus = T(-889//858) + lower_coef_plus = SVector(T(57//1430), + T(-443//2860), + T(2//65), ) + left_weights = SVector( T(25292350710583155211867593031500623221984793//80144000202690995277057891212902866905988720), + T(12366522918918608201850133537536530111621851//8904888911410110586339765690322540767332080), + T(431945550921990834855601327346624756668871//684991454723854660487674283870964674410160), + T(19955473972660125574814490349612056132376685//16028800040538199055411578242580573381197744), + T(159295736114548212938119642//176585123967931181608910805), + T(3004824804395338925078707877//2951613866135020453161224430), ) + right_weights = left_weights + left_boundary_derivatives = Tuple{}() + right_boundary_derivatives = left_boundary_derivatives + + left_boundary_minus = .- right_boundary_plus + right_boundary_minus = .- left_boundary_plus + + upper_coef_minus = .- lower_coef_plus + central_coef_minus = .- central_coef_plus + lower_coef_minus = .- upper_coef_plus + + left_boundary_central = (left_boundary_plus .+ left_boundary_minus) ./ 2 + right_boundary_central = (right_boundary_plus .+ right_boundary_minus) ./ 2 + upper_coef_central = widening_plus(upper_coef_plus, upper_coef_minus) / 2 + central_coef_central = (central_coef_plus + central_coef_minus) / 2 + lower_coef_central = widening_plus(lower_coef_plus, lower_coef_minus) / 2 + + if source.kind === :plus + left_boundary = left_boundary_plus + right_boundary = right_boundary_plus + upper_coef = upper_coef_plus + central_coef = central_coef_plus + lower_coef = lower_coef_plus + elseif source.kind === :minus + left_boundary = left_boundary_minus + right_boundary = right_boundary_minus + upper_coef = upper_coef_minus + central_coef = central_coef_minus + lower_coef = lower_coef_minus + elseif source.kind === :central + left_boundary = left_boundary_central + right_boundary = right_boundary_central + upper_coef = upper_coef_central + central_coef = central_coef_central + lower_coef = lower_coef_central + end + DerivativeCoefficients( left_boundary, right_boundary, + left_boundary_derivatives, right_boundary_derivatives, + lower_coef, central_coef, upper_coef, + left_weights, right_weights, mode, 1, order, source) + elseif order == 5 + left_boundary_plus = ( + DerivativeCoefficientRow{T,1, 6}(SVector( T(-96882324580117737612379334671496968887381092020605//59481851337402660666131943933658024760415232778698), + T(67150135620386012386529380149468175242733452110302//29740925668701330333065971966829012380207616389349), + T(-34128689447023583810334708985949827478930669509869//59481851337402660666131943933658024760415232778698), + T(-30549247576972892988497781381892473298203298542//803808801856792711704485728833216550816422064577), + T(-81720290341560256302271750655855050130279607//906210599613069573511257867906670294043316871), + T(2443823858232594607650863782592729811365949406//33529792185683574219916541112546800879602724227), )), + DerivativeCoefficientRow{T,1, 6}(SVector( T(-62974207502890835520031233971938871364620362475167//129398528148882141931797718376357864964546168868731), + T(-2014745257661048877151542919979509115504239523621//86265685432588094621198478917571909976364112579154), + T(65893029085576154023117249860090254268607930432409//129398528148882141931797718376357864964546168868731), + T(-145900185403143985727975021530756459242897398569//6994515035074710374691768560884208917002495614526), + T(53629508462843620401302448983625062085423008//1314264380885890713019873837069561991169202483), + T(-2793405615514967063407214012901800983632581620//145883346278333869145205995914721381019781475613), )), + DerivativeCoefficientRow{T,1, 7}(SVector( T(42828397538990176037716367499252080214733846099611//118214961210965459175614120008276252890416216134262), + T(-69817923295633568701918335049529759021086075531298//59107480605482729587807060004138126445208108067131), + T(-9889850321422102066251815583314096451021466796487//118214961210965459175614120008276252890416216134262), + T(1954632344730448306884683933225358654714502299200//1597499475823857556427217837949679093113732650463), + T(-524225905911016261508838953095736472108133200//1801014065190369285712759682017676542405561049), + T(-1441269847050379650669223756859345207454840300//66637520412043663571372108234654032069005758813), + T(-14923389126717029828607189853264529939286156//1801014065190369285712759682017676542405561049), )), + DerivativeCoefficientRow{T,1, 8}(SVector( T(-3168368585418325543534203816307743770727899447319//116339210945378926973099672836792750309398968236081), + T(17959939046155037792234208501951268102746904488043//232678421890757853946199345673585500618797936472162), + T(-62111250282499028018134580845475927427441764682551//116339210945378926973099672836792750309398968236081), + T(-933649263413230315057268285331840500550789509115//6288605997047509566113495829015824341048592877626), + T(2798763298980180796353446689888325092994405280//3544873730015507083491260332026958478606873099), + T(-29494146476930109297920914510358254243074591216//131160328010573762089176632284997463708454304663), + T(255054286892981964343468335673975602598708848//3544873730015507083491260332026958478606873099), + T(-14923389126717029828607189853264529939286156//3544873730015507083491260332026958478606873099), )), + DerivativeCoefficientRow{T,1, 9}(SVector( T(-3644478718032016407646992871//196481820033750264422977074804), + T(-9640303890867238936682396837//294722730050625396634465612206), + T(3850603396841089396709474000//16373485002812522035248089567), + T(-9149018907977344036501850080//16373485002812522035248089567), + T(-419036161755782772084170213183//589445460101250793268931224412), + T(49760345160398052799281005525//32746970005625044070496179134), + T(-8626411949233527849912615766//16373485002812522035248089567), + T(4855585169029650406537640012//49120455008437566105744268701), + T(-25827580686327927694349149//4465495909857960555067660791), )), + DerivativeCoefficientRow{T,1,10}(SVector( T(34725115704461833309317221//2663315948188497036110182846), + T(-14585214695327843265197651//2663315948188497036110182846), + T(-47294554476800957543403165//2663315948188497036110182846), + T(-87169805540447359252082175//1331657974094248518055091423), + T(-135369399420288825414436800//1331657974094248518055091423), + T(-4735428907200742461567397635//5326631896376994072220365692), + T(1934909297289556650463662225//1331657974094248518055091423), + T(-1250825236054281137912638419//2663315948188497036110182846), + T(117342886416269887189409193//1331657974094248518055091423), + T(-27463228735722739554968109//5326631896376994072220365692), )), + ) + right_boundary_plus = ( + DerivativeCoefficientRow{T,1, 6}(SVector( T(90121036560686812221188786049958681469456715780195//59481851337402660666131943933658024760415232778698), + T(-62974207502890835520031233971938871364620362475167//29740925668701330333065971966829012380207616389349), + T(2254126186262640844090335131539583169196518215769//3130623754600140035059575996508317092653433304142), + T(-3168368585418325543534203816307743770727899447319//29740925668701330333065971966829012380207616389349), + T(-47859315003533657615625122233305055473892971//906210599613069573511257867906670294043316871), + T(37738928595375789737653856004633685025463128//906210599613069573511257867906670294043316871), )), + DerivativeCoefficientRow{T,1, 6}(SVector( T(67150135620386012386529380149468175242733452110302//129398528148882141931797718376357864964546168868731), + T(-2014745257661048877151542919979509115504239523621//86265685432588094621198478917571909976364112579154), + T(-69817923295633568701918335049529759021086075531298//129398528148882141931797718376357864964546168868731), + T(17959939046155037792234208501951268102746904488043//258797056297764283863595436752715729929092337737462), + T(-28132561080836897573754311064403217603464786//1314264380885890713019873837069561991169202483), + T(-15851073920668845503571308281082269544475368//3942793142657672139059621511208685973507607449), )), + DerivativeCoefficientRow{T,1, 6}(SVector( T(-34128689447023583810334708985949827478930669509869//118214961210965459175614120008276252890416216134262), + T(65893029085576154023117249860090254268607930432409//59107480605482729587807060004138126445208108067131), + T(-9889850321422102066251815583314096451021466796487//118214961210965459175614120008276252890416216134262), + T(-62111250282499028018134580845475927427441764682551//59107480605482729587807060004138126445208108067131), + T(606793745326435362345644842761043224946488000//1801014065190369285712759682017676542405561049), + T(-51399276233967068573421640651971189728781720//1801014065190369285712759682017676542405561049), )), + DerivativeCoefficientRow{T,1, 7}(SVector( T(-30549247576972892988497781381892473298203298542//3144302998523754783056747914507912170524296438813), + T(-145900185403143985727975021530756459242897398569//6288605997047509566113495829015824341048592877626), + T(1954632344730448306884683933225358654714502299200//3144302998523754783056747914507912170524296438813), + T(-933649263413230315057268285331840500550789509115//6288605997047509566113495829015824341048592877626), + T(-1441739612495089163036064889179753501085176960//3544873730015507083491260332026958478606873099), + T(-4621235288646505744746116943180163298758800//86460334878427002036372203220169718990411539), + T(70546930417207777371597624760886868803898192//3544873730015507083491260332026958478606873099), )), + DerivativeCoefficientRow{T,1, 8}(SVector( T(-6222986245399120158716767907//196481820033750264422977074804), + T(9188725434809734501177200968//147361365025312698317232806103), + T(-3326642816542419513745601100//16373485002812522035248089567), + T(17760445866510407196023236440//16373485002812522035248089567), + T(-419036161755782772084170213183//589445460101250793268931224412), + T(-5601512526578427503269971200//49120455008437566105744268701), + T(-5630412589619488237368114482//49120455008437566105744268701), + T(1343034195689052240106155748//49120455008437566105744268701), )), + DerivativeCoefficientRow{T,1, 9}(SVector( T(8994644988289303196582066593//394170760331897561344307061208), + T(-5140651142922821739320165555//197085380165948780672153530604), + T(-2652341444918970761354282325//197085380165948780672153530604), + T(-13569378982598338979535684081//49271345041487195168038382651), + T(14430448220578584527721041325//10653263792753988144440731384), + T(-4735428907200742461567397635//5326631896376994072220365692), + T(-127329515047441792482124869//1331657974094248518055091423), + T(-272135630199434419226502171//2663315948188497036110182846), + T(32456543051308692201325947//1331657974094248518055091423), )), + ) + upper_coef_plus = SVector(T(31//21), + T(-167//350), + T(47//525), + T(-11//2100), ) + central_coef_plus = T(-127//140) + lower_coef_plus = SVector(T(-17//175), + T(-109//1050), + T(13//525), ) + left_weights = SVector( T(906210599613069573511257867906670294043316871//2849010651464160240006827153805046624772811600), + T(1314264380885890713019873837069561991169202483//949670217154720080002275717935015541590937200), + T(1801014065190369285712759682017676542405561049//2849010651464160240006827153805046624772811600), + T(3544873730015507083491260332026958478606873099//2849010651464160240006827153805046624772811600), + T(16373485002812522035248089567//18079306480429549386044404300), + T(1331657974094248518055091423//1310745007841312569668932475), ) + right_weights = left_weights + left_boundary_derivatives = Tuple{}() + right_boundary_derivatives = left_boundary_derivatives + + left_boundary_minus = .- right_boundary_plus + right_boundary_minus = .- left_boundary_plus + + upper_coef_minus = .- lower_coef_plus + central_coef_minus = .- central_coef_plus + lower_coef_minus = .- upper_coef_plus + + left_boundary_central = (left_boundary_plus .+ left_boundary_minus) ./ 2 + right_boundary_central = (right_boundary_plus .+ right_boundary_minus) ./ 2 + upper_coef_central = widening_plus(upper_coef_plus, upper_coef_minus) / 2 + central_coef_central = (central_coef_plus + central_coef_minus) / 2 + lower_coef_central = widening_plus(lower_coef_plus, lower_coef_minus) / 2 + + if source.kind === :plus + left_boundary = left_boundary_plus + right_boundary = right_boundary_plus + upper_coef = upper_coef_plus + central_coef = central_coef_plus + lower_coef = lower_coef_plus + elseif source.kind === :minus + left_boundary = left_boundary_minus + right_boundary = right_boundary_minus + upper_coef = upper_coef_minus + central_coef = central_coef_minus + lower_coef = lower_coef_minus + elseif source.kind === :central + left_boundary = left_boundary_central + right_boundary = right_boundary_central + upper_coef = upper_coef_central + central_coef = central_coef_central + lower_coef = lower_coef_central + end + DerivativeCoefficients( left_boundary, right_boundary, + left_boundary_derivatives, right_boundary_derivatives, + lower_coef, central_coef, upper_coef, + left_weights, right_weights, mode, 1, order, source) + elseif order == 6 + left_boundary_plus = ( + DerivativeCoefficientRow{T,1, 8}(SVector( T(-10190236575107321762019143736566294294054131524341443320227005023241702591640711//5894080668304428053004462044979693842898335662835353883451657852358878379957980), + T(33849616183326212095840203384095875150272728489801320422924473195174503084450//14033525400724828697629671535665937721186513482941318770122994886568758047519), + T(-14756990225446231802868362677209096386626752839972063725488994830567228382400//42100576202174486092889014606997813163559540448823956310368984659706274142557), + T(-577782220691921029128956803673992150878007652982133817852810//1139747145076559125170571984701210219412867300185306608134863), + T(-22130255370700075839343660989324534016694595891384661639248252556105//732554278358794216898758240608410768840028581733452373404775329129436), + T(12061110918173060345373662673278029589584871274329967646371824523331625596//45791359802234594401663057001302820495496563464024316195746122101050983405), + T(-241539699722156267393936459940614547277635169144350000//9332474760507988611613909984697980130625228656933410369), + T(-167915699824312121476321137082881057458113179282336204181029080//5087546717798786806696653843204678962742836995430687455965287991), )), + DerivativeCoefficientRow{T,1, 8}(SVector( T(-10778063091600688787531674067163947116264383321006560097506204249638414692838//24281561194595593666082827936996758400653537085704287941380525629605621345085), + T(-2504257501903645414594428064616987136329388735336980697085512957205781619077//124876600429348767425568829390269043203361047869336337984242703237971766917580), + T(522174380355188005949548246569724709878524331906166304223209404967913615700//2081276673822479457092813823171150720056017464488938966404045053966196115293), + T(4696073686577061086523934757926059890686315713230063178083//18781446408335918352353026449007753604108126697727327759429), + T(2564688071343521044815934066063373168615206761482482809668807301802//27160818260575318641757169860666699573888719356902084360327395436773), + T(-1540813932637229597064931258811873730020094387805535901477397853726176561//9054934408625100966251093422541443202332031605346699875588623249798547380), + T(69000445907020468653359916483762295249825374916568832//2306792896943663360970918770911561768172655518156590405), + T(6660985033230050960964546840872892339668707897859107381230928//754521191816122906512807605209891715234598208551437274291904677), )), + DerivativeCoefficientRow{T,1, 8}(SVector( T(5203451144742395947914637709323581713288592179367902564153666915140812645507//11968082974104968031966928413200852399227233029868857778065568216063925273130), + T(-261773961779589948581045603565717026651865712020098158151784273116590706964//170972613915785257599527548760012177131817614712412253972365260229484646759), + T(-46899447741653031874756366546060913088167964891146847363109386493631626051//1025835683494711545597165292560073062790905688274473523834191561376907880554), + T(11008175217168535988181702767728168513246655218544761410351//4628572009585186296497557478226045354577086887767955746781), + T(-1104336572340638906512036152633011635494172715077410417703085286230//743735187879244750179757393088731980958119486834708405920349829333), + T(-12348258906737384419291729230900546969770755415585403131503608023106467//557883230092838560798980472351573342827336136760100893971172265269147205), + T(40620079991051268625810390040772639140768866884918204//113698985717781971959457551849125386390652850420495609), + T(-1788293891060940633801974356182748663603745646055806138996622//20660789927996421809384072615701623801530227662077330956756117), )), + DerivativeCoefficientRow{T,1, 9}(SVector( T(10429356677491351621591110982620542691792538079628633156915830817823187326106//129237912644109975667974451674311753872090000276178466946435650008603134551615), + T(-14803989657513543667646839581829040039059381211705452899636976489520271674//72402191957484580206148152198493979760274509958643398849543781517424725239), + T(-1219750844773001309344516354781628097986388790164928195024119194257418211531//3692511789831713590513555762123192967774000007890813341326732857388660987189), + T(-2473823990023511562322244850271871930741663486004654315887//199927419064865975463071722294123710675234037336631469184302), + T(4525628116417329186104717270483447143619279124211502949994529099377//16062519527890687547780648054508813430753424886418832451772728847743), + T(3271516476790879957448755516432433094390179267678901374712403119431581597//8032438089692655189283349493415690597724603019122935265013558532496543370), + T(-245060249614038620470090633473868095785711660651967600//818522599690756243349430192726111796225380289930283513), + T(35457539706339958784099650800861784243083718000252993846218344//446213043417238885766472603383338095804567255715701754805334207), + T(-2156274690904298987224821447963907322348744126884560//818522599690756243349430192726111796225380289930283513), )), + DerivativeCoefficientRow{T,1,10}(SVector( T(-33603711929248584617691521623986535174769364561215487822300933//330299784504341369600960246454453096407859191925957125021168636), + T(-43752736219233495153747386791335506748884218588263813620457200//247724838378256027200720184840839822305894393944467843765876477), + T(84693281169574411235593761171609288090326307000610569980448519//82574946126085342400240061613613274101964797981489281255292159), + T(-423772881352053074370563780387889980031833032628943000//365089939389748562274173538583206470381976103350108947), + T(-370248443994833842319600811852409627252571620164210333689107141//990899353513024108802880739363359289223577575777871375063505908), + T(53734365297335130761642550483260570472930703641963235200//438116831900366317376868591995910769494234299047041714667), + T(425670744750746818340677394608350308506482251216572800//365089939389748562274173538583206470381976103350108947), + T(-867772845280365241571824239535093296580891931611212237600//1678318451374674140774375756867000144345944147100450829359), + T(9990739401189918640808339375566103926882514454565128//365089939389748562274173538583206470381976103350108947), + T(-4312549381808597974449642895927814644697488253769120//365089939389748562274173538583206470381976103350108947), )), + DerivativeCoefficientRow{T,1,11}(SVector( T(-9761074702163836995498860772533361649227383750020861883827090961214//193423055963891860585731178825703288009963888171893085918701622261945), + T(27385466332874798864294075507095630295260701920718960901215575167755//154738444771113488468584943060562630407971110537514468734961297809556), + T(240690444104463026571148494576149913895645505949172376676845483630//116053833578335116351438707295421972805978332903135851551220973357167), + T(-18557115617391487513955535863322182806865526830381211816000//33838344148329054577844146082124429844515166950405543430021), + T(2999917163720605965976532839479786062217690266528144100000//24039862219996182353922487627824820937218382750071942851449), + T(-613672798393387563604786645445144099644446331577325763250117126295157//2321076671566702327028774145908439456119566658062717031024419467143340), + T(578838417322931331335373799884612173351651729994672000//831225138134787260258030069078690949041125229074787969), + T(-953328245843170660415497649807533680070786119297970968000//4693928355047143658677095800087367789235234168585327660943), + T(356504082229510765887837146063366010628325695644913920//5818575966943510821806210483550836643287876603523515783), + T(49953697005949593204041696877830519634412572272825640//5818575966943510821806210483550836643287876603523515783), + T(-21562746909042989872248214479639073223487441268845600//5818575966943510821806210483550836643287876603523515783), )), + DerivativeCoefficientRow{T,1,12}(SVector( T(3763886388709655470705217623826631133//67563692944975151075524114453069695540), + T(-6961648168952882024293793425163113747//64346374233309667690975347098161614800), + T(-22192573574159038134860743266589290//160865935583274169227438367745404037), + T(83236646051757970783660064721433590//160865935583274169227438367745404037), + T(-78323897995521960295235711956272306//160865935583274169227438367745404037), + T(-134437121123763992247478559012971935//2573854969332386707639013883926464592), + T(-3640315397886009911900969493378722343//4021648389581854230685959193635100925), + T(42182341844716490888368602504391516733//27025477177990060430209645781227878216), + T(-431000890918868728712438178135066926//804329677916370846137191838727020185), + T(206347916887798155831437583354240536//2412989033749112538411575516181060555), + T(115654679223402998631329935831610623//9651956134996450153646302064724242220), + T(-832048051966928047707409610299357//160865935583274169227438367745404037), )), + DerivativeCoefficientRow{T,1,13}(SVector( T(-4253142829266554424171900//328734601968698139864969403), + T(2499764893488171595737913//140886257986584917084986887), + T(1824467274137859818131350//46962085995528305694995629), + T(-10365457883474672990469675//93924171991056611389991258), + T(186680298876494887250325//1092141534779728039418503), + T(-11678769280707117659824950//46962085995528305694995629), + T(-1538817743362627278958050//46962085995528305694995629), + T(-3522043255326282016490092819//3944815223624377678379632836), + T(69241856625979803037911030//46962085995528305694995629), + T(-22942824135772412349874998//46962085995528305694995629), + T(3661403079582727324819176//46962085995528305694995629), + T(2052157371217738298991393//187848343982113222779982516), + T(-221455831426374636581805//46962085995528305694995629), )), + ) + right_boundary_plus = ( + DerivativeCoefficientRow{T,1, 8}(SVector( T(9828684694031336445928090771707221772051618911454050836644257938873070465799289//5894080668304428053004462044979693842898335662835353883451657852358878379957980), + T(-32334189274802066362595022201491841348793149963019680292518612748915244078514//14033525400724828697629671535665937721186513482941318770122994886568758047519), + T(5203451144742395947914637709323581713288592179367902564153666915140812645507//14033525400724828697629671535665937721186513482941318770122994886568758047519), + T(20858713354982703243182221965241085383585076159257266313831661635646374652212//42100576202174486092889014606997813163559540448823956310368984659706274142557), + T(-1176129917523700461619203256839528731116927759642542073780532655//8443164463639964234740133935069220391641927800254167944916327572), + T(-68327522915146858968492025407733531544591686250146033186789636728498//310233259501735021657168600918022130278494092017264663968524502219135), + T(1625701973226982433951287830313673690261040994259512440//9332474760507988611613909984697980130625228656933410369), + T(-2898844851515363135435056563456103204204126526597136000//65327323323555920281297369892885860914376600598533872583), )), + DerivativeCoefficientRow{T,1, 8}(SVector( T(967131890952177488452577239545596432864935099708609154940699234147842945270//2081276673822479457092813823171150720056017464488938966404045053966196115293), + T(-2504257501903645414594428064616987136329388735336980697085512957205781619077//124876600429348767425568829390269043203361047869336337984242703237971766917580), + T(-174515974519726632387363735710478017767910474680065438767856182077727137976//693758891274159819030937941057050240018672488162979655468015017988732038431), + T(-503335648355460484699992545782187361328018961197985398587657200643689236916//2081276673822479457092813823171150720056017464488938966404045053966196115293), + T(-43752736219233495153747386791335506748884218588263813620457200//939138282236966862893992941484274387949542524701845868411444813), + T(1825697755524986590952938367139708686350713461381264060081038344517//12269309443067012142369861618722442229943878654697668577994516710092), + T(-150344138961359886330663274805973355513616893456406298//2306792896943663360970918770911561768172655518156590405), + T(48679501806947432707588766885451198992903547326908192//4152227214498594049747653787640811182710779932681862729), )), + DerivativeCoefficientRow{T,1, 8}(SVector( T(-210814146077803311469548038245844234094667897713886624649842783293817548320//512917841747355772798582646280036531395452844137236761917095780688453940277), + T(261087190177594002974774123284862354939262165953083152111604702483956807850//170972613915785257599527548760012177131817614712412253972365260229484646759), + T(-46899447741653031874756366546060913088167964891146847363109386493631626051//1025835683494711545597165292560073062790905688274473523834191561376907880554), + T(-1219750844773001309344516354781628097986388790164928195024119194257418211531//512917841747355772798582646280036531395452844137236761917095780688453940277), + T(84693281169574411235593761171609288090326307000610569980448519//51432190303187631837056629652413953940605061155195768190612346), + T(24069044410446302657114849457614991389564550594917237667684548363//2267771915582360361777120271342344028891023096116342732754099572241), + T(-57512645586835349315904484993635049994803967599063200//113698985717781971959457551849125386390652850420495609), + T(17764502213682011513243886170744751674321219241679200//113698985717781971959457551849125386390652850420495609), )), + DerivativeCoefficientRow{T,1, 8}(SVector( T(-8254031724170300416127954338199887869685823614030483112183//99963709532432987731535861147061855337617018668315734592151), + T(14088221059731183259571804273778179672058947139690189534249//66642473021621991821023907431374570225078012445543823061434), + T(11008175217168535988181702767728168513246655218544761410351//33321236510810995910511953715687285112539006222771911530717), + T(-2473823990023511562322244850271871930741663486004654315887//199927419064865975463071722294123710675234037336631469184302), + T(-211886440676026537185281890193944990015916516314471500//818522599690756243349430192726111796225380289930283513), + T(-12989980932174041259768875104325527964805868781266848271200//33321236510810995910511953715687285112539006222771911530717), + T(215709985514513177997480213811684530984364475133007200//818522599690756243349430192726111796225380289930283513), + T(-2968428133822149873770967110901931999527863006141400//48148388217103308432319423101535988013257664113546089), )), + DerivativeCoefficientRow{T,1, 9}(SVector( T(-632293010591430738266961742552129543334131311182418903978521501603//28657800202950170250688113863127714003635087069071818038211654365268), + T(2564688071343521044815934066063373168615206761482482809668807301802//7164450050737542562672028465781928500908771767267954509552913591317), + T(-32965270816138474821254810526358556283408140748579415453823441380//35644030103171853545631982416825514929894386901830619450512007917), + T(9051256232834658372209434540966894287238558248423005899989058198754//7164450050737542562672028465781928500908771767267954509552913591317), + T(-370248443994833842319600811852409627252571620164210333689107141//990899353513024108802880739363359289223577575777871375063505908), + T(4199884029208848352367145975271700487104766373139401740000//10558766137090918169531372909364914329917130884988500856187), + T(-405956936120383711129131210063334564894745960745744960//365089939389748562274173538583206470381976103350108947), + T(156319768597662826013897313951478637584954262987186400//365089939389748562274173538583206470381976103350108947), + T(-5390686727260747468062053619909768305871860317211400//365089939389748562274173538583206470381976103350108947), )), + DerivativeCoefficientRow{T,1,10}(SVector( T(1723015845453294335053380381896861369940695896332852520910260646190232228//28549823329438330298035679188210282420134699785835935160858115550729867835), + T(-220116276091032799580704465544553390002870626829362271639628264818025223//1087612317311936392306121111931820282671798087079464196604118687646852108), + T(-24696517813474768838583458461801093939541510831170806263007216046212934//5709964665887666059607135837642056484026939957167187032171623110145973567), + T(3271516476790879957448755516432433094390179267678901374712403119431581597//5709964665887666059607135837642056484026939957167187032171623110145973567), + T(1573649269421957400876674692724059563850113463800351888000//40897108021369667991955337428740673383772402395708642862769), + T(-613672798393387563604786645445144099644446331577325763250117126295157//2321076671566702327028774145908439456119566658062717031024419467143340), + T(-31106908977700033041336860796603304679530489988360250//831225138134787260258030069078690949041125229074787969), + T(-162448611462872476636001170341953910213162079624872000//831225138134787260258030069078690949041125229074787969), + T(214189952629827032730998930497748127353308583270532960//5818575966943510821806210483550836643287876603523515783), + T(-26953433636303737340310268099548841529359301586057000//5818575966943510821806210483550836643287876603523515783), )), + DerivativeCoefficientRow{T,1,11}(SVector( T(-7988882895681781717303491340527875//965195613499645015364630206472424222), + T(199690536335522891947980787509571528//4021648389581854230685959193635100925), + T(940451378570533970096091476278074353//9651956134996450153646302064724242220), + T(-94562118716135158524736058218555095//160865935583274169227438367745404037), + T(82127410631674873637636848409605580//160865935583274169227438367745404037), + T(156350656851900249410278189778575380//160865935583274169227438367745404037), + T(-3640315397886009911900969493378722343//4021648389581854230685959193635100925), + T(-5781605738039420622199779572229570//160865935583274169227438367745404037), + T(-331155124682837362987549024899144086//2412989033749112538411575516181060555), + T(123975159743072279108404031934604193//2412989033749112538411575516181060555), + T(-4160240259834640238537048051496785//643463742333096676909753470981616148), )), + DerivativeCoefficientRow{T,1,12}(SVector( T(-4197892495607803036908028427072026436452829482058405104525727//436228443028051414848261247307437193320249950506000384943419316), + T(416311564576878185060284177554555771229294243616194211326933//31159174502003672489161517664816942380017853607571456067387094), + T(-894146945530470316900987178091374331801872823027903069498311//41545566002671563318882023553089256506690471476761941423182792), + T(4432192463292494848012456350107723030385464750031624230777293//31159174502003672489161517664816942380017853607571456067387094), + T(-579298558474153488722063336325//2806501221178767076638633784669), + T(-17532658174026079978181501850//67840555715912366624559150851), + T(2277430469806238926656543465270947291800379554839515055//1600414172171538614368076667869805394163618902864752088), + T(-3522043255326282016490092819//3944815223624377678379632836), + T(-4960610623950791859432432//46962085995528305694995629), + T(-5875961393846473690637226//46962085995528305694995629), + T(2199794592168654723379263//46962085995528305694995629), + T(-1107279157131873182909025//187848343982113222779982516), )), + ) + upper_coef_plus = SVector(T(67//45), + T(-37//75), + T(124//1575), + T(139//12600), + T(-1//210), ) + central_coef_plus = T(-8//9) + lower_coef_plus = SVector(T(-8//75), + T(-199//1575), + T(149//3150), + T(-1//168), ) + left_weights = SVector( T(9332474760507988611613909984697980130625228656933410369//31697237956293195112204875285069437638526538665203032000), + T(461358579388732672194183754182312353634531103631318081//301878456726601858211475002714947025128824177763838400), + T(113698985717781971959457551849125386390652850420495609//452817685089902787317212504072420537693236266645757600), + T(818522599690756243349430192726111796225380289930283513//452817685089902787317212504072420537693236266645757600), + T(365089939389748562274173538583206470381976103350108947//905635370179805574634425008144841075386472533291515200), + T(831225138134787260258030069078690949041125229074787969//646882407271289696167446434389172196704623238065368000), + T(160865935583274169227438367745404037//174730090913054890018556018162864970), + T(46962085995528305694995629//46505724599538673682179050), ) + right_weights = left_weights + left_boundary_derivatives = Tuple{}() + right_boundary_derivatives = left_boundary_derivatives + + left_boundary_minus = .- right_boundary_plus + right_boundary_minus = .- left_boundary_plus + + upper_coef_minus = .- lower_coef_plus + central_coef_minus = .- central_coef_plus + lower_coef_minus = .- upper_coef_plus + + left_boundary_central = (left_boundary_plus .+ left_boundary_minus) ./ 2 + right_boundary_central = (right_boundary_plus .+ right_boundary_minus) ./ 2 + upper_coef_central = widening_plus(upper_coef_plus, upper_coef_minus) / 2 + central_coef_central = (central_coef_plus + central_coef_minus) / 2 + lower_coef_central = widening_plus(lower_coef_plus, lower_coef_minus) / 2 + + if source.kind === :plus + left_boundary = left_boundary_plus + right_boundary = right_boundary_plus + upper_coef = upper_coef_plus + central_coef = central_coef_plus + lower_coef = lower_coef_plus + elseif source.kind === :minus + left_boundary = left_boundary_minus + right_boundary = right_boundary_minus + upper_coef = upper_coef_minus + central_coef = central_coef_minus + lower_coef = lower_coef_minus + elseif source.kind === :central + left_boundary = left_boundary_central + right_boundary = right_boundary_central + upper_coef = upper_coef_central + central_coef = central_coef_central + lower_coef = lower_coef_central + end + DerivativeCoefficients( left_boundary, right_boundary, + left_boundary_derivatives, right_boundary_derivatives, + lower_coef, central_coef, upper_coef, + left_weights, right_weights, mode, 1, order, source) + elseif order == 7 + left_boundary_plus = ( + DerivativeCoefficientRow{T,1, 8}(SVector( T(-13033386346315743757800580520477868798421453633825229808543341244803220668398489275683839//9303484962498927924217104006028014429616200818506730195522770559309308246047703285543740), + T(94534600595942768184202885243545592573105350049495265538253571674514882396997487988876//66453464017849485172979314328771531640115720132190929968019789709352201757483594896741), + T(25029185149587713148452074085971367459020298148603520885647890266467793543611830183300//66453464017849485172979314328771531640115720132190929968019789709352201757483594896741), + T(-3988245245869599382342757973690628281240162671728753604781287877710//21358550461055418572524515226933712412739104587742670505709499725579), + T(-52157423959894608944973301855104987466410153171243153529022947293888251801017//191609505828893758073804733504642375625829589620470614068084997381987797976732), + T(-106605000205582851041349790209086779318888396869523693703320808170537131183838//428552095421252129104334762915725513795316459914038297848258582625572975770845), + T(6424530485059236141989612161728634325383288251158263495972//14314108189450995827895179823054016375253568054036851156531), + T(-30211820855294239685121615490232423847886239838993797556058293910816//216680925024219527013296082822218682938851358246995609423580833054417), )), + DerivativeCoefficientRow{T,1, 8}(SVector( T(-135146918110229697743023630544897352321734834549107741835094115073161535903039839554558//329813515401639007346851235924351255641537079124235100530522534418861169624635783183735), + T(-8611401568934025643746009711155941752158662729292510906734188887942690704148392943017//113078919566276231090348995174063287648526998556880605896179154657895258157017982805852), + T(2388663556720989147127940850279183797138183469218141745696071258249908195955797325726//9423243297189685924195749597838607304043916546406717158014929554824604846418165233821), + T(654816496895358115744162237611885257011113997683883652050775214915//3028688127014230968848197157022190352663681738030253696509108630099), + T(2688494260040757369918491293737609729641396466810264490285611829592842708466//20377978238218356085253882695579636943762406326817445914925018097791823370669), + T(4750792864453494801708115725276626651301942702876254016217004028705884509439//243078414066734139263321598716241892706648267554842503887485349839990983357780), + T(-1472614277041425963426205781584773607163380248787594349190//6089313448668155990224990402356718810001747782821362176433), + T(9825267112423038163423866576734043031591416885048180829483855439730//92177455511494366958428324629004026919950025868273329645345181579131), )), + DerivativeCoefficientRow{T,1, 8}(SVector( T(-8179893666037024061449426919024537106302008587900670151735602015520149963002921298257//71376980709683814954392615120030484348451483831792609051503942425515918046301263166466), + T(-986265724112681035018927869386854795232109613919161668705028438472813393923317316627//5098355764977415353885186794287891739175105987985186360821710173251137003307233083319), + T(-1797343214906560683369889884974996500110898522265708169777303743204384365008913245523//30590134589864492123311120765727350435050635927911118164930261039506822019843398499914), + T(190896110009019417896364217622484020902396605298822966938107732445//1638642777830723707545643914852887655810037103866870612260858145561), + T(618293072006416119852416591997676308764190166168717775369575850397535001426//3675103484465532682939557887290052396733906668232602048760748175739677419597), + T(3791719511567593760224740815940444617084570443077195546144867687511775649248//19727288067606729774788221700893733040806326495720571692081052849094060829113), + T(-324318072235400109801833123693924239038220654795438909200//3294564869722780234528662350831441953388554779509986524787), + T(-583325182747644756612052829972381162288735910491552396053945831638//49871731726180917348261310854476248628279119939258672401486094123409), )), + DerivativeCoefficientRow{T,1, 9}(SVector( T(495092963247565489562008647927384516075798052996069445979534829934850995204791408951//10019791574268785829163929830415717685921930685986937500537653579870832917619540237875), + T(-2453709997248680546377294894770001173560462416138912851447690857245076685448829220229//10019791574268785829163929830415717685921930685986937500537653579870832917619540237875), + T(-922703647373417149751432119850770162696591969285493578270021964997659455910022607897//10019791574268785829163929830415717685921930685986937500537653579870832917619540237875), + T(-118406792849597143876793929463055206753672470185553418794360217681//1288168959203059103203232055635200807114301232338714804971008932850), + T(911615481086333559914009872946299095322387070176645429640544748940560863653//7222675824462932822572891522916071600833463298659991609928949672607879798625), + T(234408699218018902887694451228449634383531789393200068286430398539497849289//1033866989671742667945979499496639616523320707196594848054191975376762276270), + T(10542246789883503272525229721000100441805781611243219072//86330717442394841356459447413383751091710043106405202365), + T(-458551277701112222664673116287163272826446107992979752897746297002//4667278370661829311461936327194202442721095932269744357961598846625), + T(5761000960280699594156050291875609244322038396434676491//2158267936059871033911486185334593777292751077660130059125), )), + DerivativeCoefficientRow{T,1,10}(SVector( T(6345423310918556759173882941128648675858332170190557131447849199145431//106810293391173183530634001855866678189870557771787363707764464381192940), + T(-53408486147166906456673797139595152225113700340327839599509668270031146//560754040303659213535828509743300060496820428301883659465763438001262935), + T(-28876661672685538927936577885807710391592089915077022782477600207353147//186918013434553071178609503247766686832273476100627886488587812667087645), + T(-71922806105245269610067949091932816379415113995274047253632//997789275075513785566394128324999387459122710238571202286931), + T(-63772660870741854683965488506866672617212415988898155409761516905436601//320430880173519550591902005567600034569611673315362091123293393143578820), + T(1945099026636814992982183820424467030173358352418202377851840//6888094657195160129770624721272134052988593526619786967131323), + T(30491132878264481254975353440431648676978285293973448960//138948513448755575207686134009887117039287384798575574751), + T(-132058088335222471369931948595617222011283666801138588077440//15002132048548690699598664202913502139614819649317406230290719), + T(-165185172717302230425652177432720996677087266476875961476//4863197970706445132269014690346049096375058467950145116285), + T(11522001920561399188312100583751218488644076792869352982//4863197970706445132269014690346049096375058467950145116285), )), + DerivativeCoefficientRow{T,1,11}(SVector( T(91881697306003381275551832387043092083226489538991032867696422749630274//1471182889624041302596275267240610955593379422438865495875392405143288785), + T(-26224531532713018683536394498099616885960700497981132492549732057178877//1176946311699233042077020213792488764474703537951092396700313924114631028), + T(-140810052918767147552030196779153306150442817706268161427791177052898340//882709733774424781557765160344366573356027653463319297525235443085973271), + T(-59945783600832685027379103121600599295987206016689875964424800//506168474353967852150086907368394302434842122150276061697718271), + T(11015661683542626722793535883927010803141707490135071102373760//212097766934369128870421037387247272563510429972242931127050191), + T(-7065987610069312558859776643548695751570451137162287647684932490218558253//17654194675488495631155303206887331467120553069266385950504708861719465420), + T(5139603843778379578002936126553397633111890399184071564800//7307920164502950379713366550227311875530111634642970441617), + T(-33290037069214624152083342387720539524081530396534436953515200//146370332974829593155279018634502829554992605930264054975146893), + T(959218385126140038803361252939809717217430089149649961416//7307920164502950379713366550227311875530111634642970441617), + T(-165185172717302230425652177432720996677087266476875961476//7307920164502950379713366550227311875530111634642970441617), + T(1047454720051036289846554598522838044422188799351759362//664356378591177307246669686384301079593646512240270040147), )), + DerivativeCoefficientRow{T,1,12}(SVector( T(-9318116360191386282001991582603093067667//103364409360928627901203664747100758788890), + T(553620088164233504155108503997100337457//2953268838883675082891533278488593108254), + T(5328529535893303353380441940208658610//492211473147279180481922213081432184709), + T(-42267899301897733168814765063322289200//492211473147279180481922213081432184709), + T(-186826146815920769295286382160941543760//492211473147279180481922213081432184709), + T(1004788743830377524045808400049722725636//2461057365736395902409611065407160923545), + T(-1896846955281589729183330783433054950817//1476634419441837541445766639244296554127), + T(73868049498686808818621373133748430219265//41345763744371451160481465898840303515556), + T(-362218101993197645420308737862239054324//492211473147279180481922213081432184709), + T(107201965471823313326813994946146578187//492211473147279180481922213081432184709), + T(-36922092938760190645124084705684358339//984422946294558360963844426162864369418), + T(5150782225225489400994974527368800321//1968845892589116721927688852325728738836), )), + DerivativeCoefficientRow{T,1,13}(SVector( T(693724230826827761508080759//105321124336545997237887616458), + T(-819119549882808222227494693//15045874905220856748269659494), + T(49456512091726657157585800//835881939178936486014981083), + T(41956469511136512069192000//835881939178936486014981083), + T(57191554007203266783566000//835881939178936486014981083), + T(-1001164524483276495132590800//2507645817536809458044943249), + T(60249054510843418297714000//835881939178936486014981083), + T(-84866567058497733601671303623//105321124336545997237887616458), + T(1192370714658039515438925995//835881939178936486014981083), + T(-1413095281525109975559793208//2507645817536809458044943249), + T(418219273815575787510837154//2507645817536809458044943249), + T(-24006915091897417100416723//835881939178936486014981083), + T(30141561197152293701886073//15045874905220856748269659494), )), + ) + right_boundary_plus = ( + DerivativeCoefficientRow{T,1, 8}(SVector( T(12060529458203493529220036537008447034706973719099929496295699684328606444868532418940161//9303484962498927924217104006028014429616200818506730195522770559309308246047703285543740), + T(-135146918110229697743023630544897352321734834549107741835094115073161535903039839554558//110755773363082475288298857214619219400192866886984883280032982848920336262472658161235), + T(-8179893666037024061449426919024537106302008587900670151735602015520149963002921298257//22151154672616495057659771442923843880038573377396976656006596569784067252494531632247), + T(6931301485465916853868121070983383225061172741944972243713487619087913932867079725314//66453464017849485172979314328771531640115720132190929968019789709352201757483594896741), + T(310925742235009281199520264115303785117058276339337299440944610758126119//2200662759752538309545357515357272687475790345823089895003790067440625228), + T(643171881142023668928862826709301644582585426772937230073874959247411918//2881623030152516686531880680449206313889391805445425923037800029757952755), + T(-13896058230060338262526894852341601050111485231623688205676//71570540947254979139475899115270081876267840270184255782655), + T(1856295467693724053830361573348614098445679152256045267142//100198757326156970795266258761378114626774976378257958095717), )), + DerivativeCoefficientRow{T,1, 8}(SVector( T(13504942942277538312028983606220798939015050007070752219750510239216411770999641141268//28269729891569057772587248793515821912131749639220151474044788664473814539254495701463), + T(-8611401568934025643746009711155941752158662729292510906734188887942690704148392943017//113078919566276231090348995174063287648526998556880605896179154657895258157017982805852), + T(-1972531448225362070037855738773709590464219227838323337410056876945626787846634633254//9423243297189685924195749597838607304043916546406717158014929554824604846418165233821), + T(-4907419994497361092754589789540002347120924832277825702895381714490153370897658440458//28269729891569057772587248793515821912131749639220151474044788664473814539254495701463), + T(-53408486147166906456673797139595152225113700340327839599509668270031146//702132041422952695629462243585419734133701075933481925194673813795673203), + T(-26224531532713018683536394498099616885960700497981132492549732057178877//980688738090227096457029426164411646285874437919199983408247835878365172), + T(825610744176553348016712569736017075868971604620735829796//6089313448668155990224990402356718810001747782821362176433), + T(-313119050641967950037848339416935276305830295905002956062//6089313448668155990224990402356718810001747782821362176433), )), + DerivativeCoefficientRow{T,1, 8}(SVector( T(1787798939256265224889433863283669104215735582043108634689135019033413824543702155950//15295067294932246061655560382863675217525317963955559082465130519753411009921699249957), + T(1194331778360494573563970425139591898569091734609070872848035629124954097977898662863//5098355764977415353885186794287891739175105987985186360821710173251137003307233083319), + T(-1797343214906560683369889884974996500110898522265708169777303743204384365008913245523//30590134589864492123311120765727350435050635927911118164930261039506822019843398499914), + T(-922703647373417149751432119850770162696591969285493578270021964997659455910022607897//15295067294932246061655560382863675217525317963955559082465130519753411009921699249957), + T(-28876661672685538927936577885807710391592089915077022782477600207353147//253254555660375059982741817681842152550315726715542986511439077679059878), + T(-70405026459383573776015098389576653075221408853134080713895588526449170//397944204861653113081480275571253162826666259772871758660582430942126981), + T(23839224747516106082271205029652749854321305130986713240//3294564869722780234528662350831441953388554779509986524787), + T(170148527217986853422934989259276778577587271319523354800//3294564869722780234528662350831441953388554779509986524787), )), + DerivativeCoefficientRow{T,1, 8}(SVector( T(-56974932083851419747753685338437546874859466738982194354018398253//644084479601529551601616027817600403557150616169357402485504466425), + T(130963299379071623148832447522377051402222799536776730410155042983//429389653067686367734410685211733602371433744112904934990336310950), + T(38179222001803883579272843524496804180479321059764593387621546489//214694826533843183867205342605866801185716872056452467495168155475), + T(-118406792849597143876793929463055206753672470185553418794360217681//1288168959203059103203232055635200807114301232338714804971008932850), + T(-251729821368358443635237821821764857327952898983459165387712//3099704409769186778903676459377543582947849097735478790915325), + T(-1198915672016653700547582062432011985919744120333797519288496//5979524482212593896872450706193198751865112715678943531406995), + T(-7564067679654094355695130013158961937152699842549723712//86330717442394841356459447413383751091710043106405202365), + T(1154765056030785849992278517000867683585006663731014016//17266143488478968271291889482676750218342008621281040473), )), + DerivativeCoefficientRow{T,1, 9}(SVector( T(-1064437223671318549897414323573571172783880676964145990388223414160984730633//9299865435276057916828771907588455803313839594631753970671344149206088092860), + T(2688494260040757369918491293737609729641396466810264490285611829592842708466//16274764511733101354450350838279797655799219290605569448674852261110654162505), + T(176655163430404605672119026285050373932625761762490792962735957256438571836//774988786273004826402397658965704650276153299552646164222612012433840674405), + T(260461566024666731404002820841799741520682020050470122754441356840160246758//2324966358819014479207192976897113950828459898657938492667836037301522023215), + T(-63772660870741854683965488506866672617212415988898155409761516905436601//320430880173519550591902005567600034569611673315362091123293393143578820), + T(314733190958360763508386739540771737232620214003859174353536//4032702705823233059252674667368953797831237769009058905998273), + T(-47762203948239308024810696875085858934982476516241938048//138948513448755575207686134009887117039287384798575574751), + T(11243423215456684356063579923590162005400655723693611200//138948513448755575207686134009887117039287384798575574751), + T(-6722470591372322457224156378579408344799122145093380980//972639594141289026453802938069209819275011693590029023257), )), + DerivativeCoefficientRow{T,1,10}(SVector( T(-15229285743654693005907112887012397045555485267074813386188686881505304454834//218792848161997798480815461468756420704891794325685637681092483100504764816415), + T(4750792864453494801708115725276626651301942702876254016217004028705884509439//291723797549330397974420615291675227606522392434247516908123310800673019755220), + T(7583439023135187520449481631880889234169140886154391092289735375023551298496//43758569632399559696163092293751284140978358865137127536218496620100952963283), + T(5860217480450472572192361280711240859588294734830001707160759963487446232225//43758569632399559696163092293751284140978358865137127536218496620100952963283), + T(68078465932288524754376433714856346056067542334637083224814400//362275526314904759173530719994418531605654224064155973702279541), + T(-7065987610069312558859776643548695751570451137162287647684932490218558253//17654194675488495631155303206887331467120553069266385950504708861719465420), + T(8990617214972829373695625271405697409213659851621680234848//36539600822514751898566832751136559377650558173214852208085), + T(-208749871950906571407751335531720907695078784247196811200//664356378591177307246669686384301079593646512240270040147), + T(227891753047521731299898901233841942888690240718665615222//7307920164502950379713366550227311875530111634642970441617), + T(-33612352956861612286120781892897041723995610725466904900//7307920164502950379713366550227311875530111634642970441617), )), + DerivativeCoefficientRow{T,1,11}(SVector( T(615431694249176941010392012747358985607//2953268838883675082891533278488593108254), + T(-658315759283321089689698541071957971945//1968845892589116721927688852325728738836), + T(-72491385321176112918608140222878276300//492211473147279180481922213081432184709), + T(58909920508660469900873050630201615200//492211473147279180481922213081432184709), + T(119268802458777332440377289778948935200//492211473147279180481922213081432184709), + T(574400618611084392128548068149859873600//492211473147279180481922213081432184709), + T(-1896846955281589729183330783433054950817//1476634419441837541445766639244296554127), + T(46330842209555457760405678840766601000//492211473147279180481922213081432184709), + T(-53842098050389647686927115005222846232//492211473147279180481922213081432184709), + T(101876462004222468111674007714321577041//1968845892589116721927688852325728738836), + T(-7513013422140300008235546291616635475//984422946294558360963844426162864369418), )), + DerivativeCoefficientRow{T,1,12}(SVector( T(-7552955213823559921280403872558105961971559959748449389014573477704//152361093569235137761695327078965325398585446284303870298717681636101), + T(1637544518737173027237311096122340505265236147508030138247309239955//14510580339927155977304316864663364323674804408028940028449303012962), + T(-97220863791274126102008804995396860381455985081925399342324305273//7255290169963577988652158432331682161837402204014470014224651506481), + T(-1604929471953892779326355907005071454892561377975429135142112039507//21765870509890733965956475296995046485512206612043410042673954519443), + T(-671735568997268648025425199200//90249337091210593458551492550427), + T(-14514458171639073281379921338800//50225638079444756635182168334221), + T(27539759042933016809783936661418096897440573456978268846105//20130181528155525492450616832227809669414692740616280464286), + T(-84866567058497733601671303623//105321124336545997237887616458), + T(-481528821621806114038000846//2507645817536809458044943249), + T(-70016760881308991611436848//835881939178936486014981083), + T(198721636169556302710263811//5015291635073618916089886498), + T(-43964963754326615643863675//7522937452610428374134829747), )), + ) + upper_coef_plus = SVector(T(119//80), + T(-617//1050), + T(5113//29400), + T(-587//19600), + T(737//352800), ) + central_coef_plus = T(-1111//1400) + lower_coef_plus = SVector(T(-841//4200), + T(-107//1225), + T(4859//117600), + T(-43//7056), ) + left_weights = SVector( T(14314108189450995827895179823054016375253568054036851156531//38608868308030436140373899052417922903060260580285148065600), + T(6089313448668155990224990402356718810001747782821362176433//5515552615432919448624842721773988986151465797183592580800), + T(3294564869722780234528662350831441953388554779509986524787//2757776307716459724312421360886994493075732898591796290400), + T(86330717442394841356459447413383751091710043106405202365//110311052308658388972496854435479779723029315943671851616), + T(138948513448755575207686134009887117039287384798575574751//157587217583797698532138363479256828175756165633816930880), + T(7307920164502950379713366550227311875530111634642970441617//5515552615432919448624842721773988986151465797183592580800), + T(492211473147279180481922213081432184709//616416543100255312303604821321476510600), + T(835881939178936486014981083//801593757753303875925328400), ) + right_weights = left_weights + left_boundary_derivatives = Tuple{}() + right_boundary_derivatives = left_boundary_derivatives + + left_boundary_minus = .- right_boundary_plus + right_boundary_minus = .- left_boundary_plus + + upper_coef_minus = .- lower_coef_plus + central_coef_minus = .- central_coef_plus + lower_coef_minus = .- upper_coef_plus + + left_boundary_central = (left_boundary_plus .+ left_boundary_minus) ./ 2 + right_boundary_central = (right_boundary_plus .+ right_boundary_minus) ./ 2 + upper_coef_central = widening_plus(upper_coef_plus, upper_coef_minus) / 2 + central_coef_central = (central_coef_plus + central_coef_minus) / 2 + lower_coef_central = widening_plus(lower_coef_plus, lower_coef_minus) / 2 + + if source.kind === :plus + left_boundary = left_boundary_plus + right_boundary = right_boundary_plus + upper_coef = upper_coef_plus + central_coef = central_coef_plus + lower_coef = lower_coef_plus + elseif source.kind === :minus + left_boundary = left_boundary_minus + right_boundary = right_boundary_minus + upper_coef = upper_coef_minus + central_coef = central_coef_minus + lower_coef = lower_coef_minus + elseif source.kind === :central + left_boundary = left_boundary_central + right_boundary = right_boundary_central + upper_coef = upper_coef_central + central_coef = central_coef_central + lower_coef = lower_coef_central + end + DerivativeCoefficients( left_boundary, right_boundary, + left_boundary_derivatives, right_boundary_derivatives, + lower_coef, central_coef, upper_coef, + left_weights, right_weights, mode, 1, order, source) + else + throw(ArgumentError("Order $order not implemented/derived.")) + end +end diff --git a/src/SummationByPartsOperators.jl b/src/SummationByPartsOperators.jl index 2f985c7a..303029bf 100644 --- a/src/SummationByPartsOperators.jl +++ b/src/SummationByPartsOperators.jl @@ -123,6 +123,7 @@ include("SBP_coefficients/MattssonAlmquistVanDerWeide2018Minimal.jl") include("SBP_coefficients/MattssonAlmquistVanDerWeide2018Accurate.jl") include("SBP_coefficients/DienerDorbandSchnetterTiglio2007.jl") include("SBP_coefficients/SharanBradyLivescu2022.jl") +include("SBP_coefficients/WilliamsDuru2024.jl") include("conservation_laws/general_laws.jl") include("conservation_laws/burgers.jl") @@ -168,7 +169,8 @@ export MattssonNordström2004, MattssonSvärdNordström2004, MattssonSvärdShoey Mattsson2017, MattssonAlmquistVanDerWeide2018Minimal, MattssonAlmquistVanDerWeide2018Accurate, DienerDorbandSchnetterTiglio2007, - SharanBradyLivescu2022 + SharanBradyLivescu2022, + WilliamsDuru2024 export Tadmor1989, MadayTadmor1989, Tadmor1993, TadmorWaagan2012Standard, TadmorWaagan2012Convergent export GlaubitzNordströmÖffner2023 diff --git a/test/upwind_operators_test.jl b/test/upwind_operators_test.jl index 6398421b..288ebcb4 100644 --- a/test/upwind_operators_test.jl +++ b/test/upwind_operators_test.jl @@ -2,18 +2,22 @@ using Test using LinearAlgebra using SummationByPartsOperators +D_test_list = vcat( + [(Mattsson2017, acc_order) for acc_order in 2:9], + [(WilliamsDuru2024, acc_order) for acc_order in 4:7] +) + # check construction of interior part of upwind operators @testset "Check interior parts" begin N = 21 xmin = 0.0 xmax = Float64(N + 1) - interior = 10:N-10 - - for acc_order in 2:9 - Dp_bounded = derivative_operator(Mattsson2017(:plus ), 1, acc_order, xmin, xmax, N) - Dm_bounded = derivative_operator(Mattsson2017(:minus ), 1, acc_order, xmin, xmax, N) - Dc_bounded = derivative_operator(Mattsson2017(:central), 1, acc_order, xmin, xmax, N) + # bounded derivative operators + for (source, acc_order) in D_test_list + Dp_bounded = derivative_operator(source(:plus ), 1, acc_order, xmin, xmax, N) + Dm_bounded = derivative_operator(source(:minus ), 1, acc_order, xmin, xmax, N) + Dc_bounded = derivative_operator(source(:central), 1, acc_order, xmin, xmax, N) for compact in (true, false) show(IOContext(devnull, :compact=>compact), Dp_bounded) show(IOContext(devnull, :compact=>compact), Dm_bounded) @@ -25,26 +29,21 @@ using SummationByPartsOperators M = mass_matrix(Dp_bounded) @test M == mass_matrix(Dm_bounded) @test M == mass_matrix(Dc_bounded) - Dp_periodic = periodic_derivative_operator(1, acc_order, xmin, xmax, N-1, -(acc_order - 1) ÷ 2) - Dm_periodic = periodic_derivative_operator(1, acc_order, xmin, xmax, N-1, -acc_order + (acc_order - 1) ÷ 2) Dp = Matrix(Dp_bounded) Dm = Matrix(Dm_bounded) - @test Dp[interior,interior] ≈ Matrix(Dp_periodic)[interior,interior] - @test Dm[interior,interior] ≈ Matrix(Dm_periodic)[interior,interior] - - D_periodic = upwind_operators(periodic_derivative_operator; - accuracy_order = acc_order, - xmin, xmax, N = N - 1) - @test D_periodic.minus == Dm_periodic - @test D_periodic.plus == Dp_periodic - @test Matrix(D_periodic.central) ≈ (Matrix(Dm_periodic) + Matrix(Dp_periodic)) / 2 + Dc = Matrix(Dc_bounded) + @test Dc ≈ (Dm + Dp) / 2 + # upwind summation by parts res = M * Dp + Dm' * M res[1,1] += 1 res[end,end] -= 1 @test norm(res) < N * eps() - x = grid(Dp_bounded) + + # derivative accuracy order for D in (Dp_bounded, Dm_bounded, Dc_bounded) + x = grid(D) + interior = length(D.coefficients.left_boundary)+1:N-length(D.coefficients.right_boundary)-1 @test norm(D * x.^0) < N * eps() for k in 1:acc_order÷2 @test D * x.^k ≈ k .* x.^(k-1) @@ -53,9 +52,32 @@ using SummationByPartsOperators @test (D * x.^k)[interior] ≈ (k .* x.^(k-1))[interior] end end + + # dissipation matrix diss = M * (Dp - Dm) @test diss ≈ diss' @test maximum(eigvals(Symmetric(diss))) < N * eps() + end + + # periodic derivative operators + for acc_order in 2:9 + Dp_bounded = derivative_operator(Mattsson2017(:plus ), 1, acc_order, xmin, xmax, N) + Dm_bounded = derivative_operator(Mattsson2017(:minus ), 1, acc_order, xmin, xmax, N) + Dc_bounded = derivative_operator(Mattsson2017(:central), 1, acc_order, xmin, xmax, N) + + Dp_periodic = periodic_derivative_operator(1, acc_order, xmin, xmax, N-1, -(acc_order - 1) ÷ 2) + Dm_periodic = periodic_derivative_operator(1, acc_order, xmin, xmax, N-1, -acc_order + (acc_order - 1) ÷ 2) + + interior = length(Dc_bounded.coefficients.left_boundary)+1:N-length(Dc_bounded.coefficients.right_boundary)-1 + @test Matrix(Dp_bounded)[interior,interior] ≈ Matrix(Dp_periodic)[interior,interior] + @test Matrix(Dm_bounded)[interior,interior] ≈ Matrix(Dm_periodic)[interior,interior] + + D_periodic = upwind_operators(periodic_derivative_operator; + accuracy_order = acc_order, + xmin, xmax, N = N - 1) + @test D_periodic.minus == Dm_periodic + @test D_periodic.plus == Dp_periodic + @test Matrix(D_periodic.central) ≈ (Matrix(Dm_periodic) + Matrix(Dp_periodic)) / 2 # mass matrix scaling x1 = grid(D_periodic) @@ -69,6 +91,44 @@ using SummationByPartsOperators end end +@testset "Rational Tests" begin + N = 21 + xmin = big(0 // 1) + xmax = big((N + 1) // 1) + + # bounded derivative operators + for (source, acc_order) in D_test_list + Dp_bounded = derivative_operator(source(:plus ), 1, acc_order, xmin, xmax, N) + Dm_bounded = derivative_operator(source(:minus ), 1, acc_order, xmin, xmax, N) + Dc_bounded = derivative_operator(source(:central), 1, acc_order, xmin, xmax, N) + + M = mass_matrix(Dp_bounded) + Dp = Matrix(Dp_bounded) + Dm = Matrix(Dm_bounded) + Dc = Matrix(Dc_bounded) + @test Dc == (Dm + Dp) / 2 + + # upwind summation by parts + res = M * Dp + Dm' * M + res[1,1] += 1 + res[end,end] -= 1 + @test all(isequal(0), res) + + # derivative accuracy order + for D in (Dp_bounded, Dm_bounded, Dc_bounded) + x = grid(D) + interior = length(D.coefficients.left_boundary)+1:N-length(D.coefficients.right_boundary)-1 + @test all(isequal(0), D * x.^0) + for k in 1:acc_order÷2 + @test D * x.^k == k .* x.^(k-1) + end + for k in acc_order÷2+1:acc_order + @test (D * x.^k)[interior] == (k .* x.^(k-1))[interior] + end + end + end +end + @testset "UpwindOperators" begin N = 21 xmin = 0.