From e185521e57c696f31bf9f446cf2dd8b74a74077c Mon Sep 17 00:00:00 2001 From: Chris Coutinho Date: Wed, 9 Nov 2016 06:22:34 +0100 Subject: [PATCH] Implement additional tableau - fehlbery_rk78 --- .gitignore | 1 + data.out | 198 +++++++++++++++++++++---------------------- plotter.py | 1 + src/main.f90 | 14 ++- src/misc.f90 | 4 +- src/rk_constants.f90 | 75 ++++++++++++---- src/runge_kutta.f90 | 54 +++++++----- 7 files changed, 205 insertions(+), 142 deletions(-) diff --git a/.gitignore b/.gitignore index 3b084d5..c36e7f4 100644 --- a/.gitignore +++ b/.gitignore @@ -1,5 +1,6 @@ .ipynb_checkpoints data.out +raw.out plot.png main diff --git a/data.out b/data.out index 055f01f..5f53a01 100644 --- a/data.out +++ b/data.out @@ -1,100 +1,100 @@ 0.0000000000000000 1.5000000000000000 3.0000000000000000 - 0.20202020202020202 1.9206248596079809 2.4378529935131867 - 0.40404040404040403 2.3324513659965129 1.7954446825547759 - 0.60606060606060608 2.4323743966638784 1.4098939240174651 - 0.80808080808080807 2.2439625918424486 1.3248863394318924 - 1.0101010101010102 1.9535462985358967 1.3927036539584361 - 1.2121212121212122 1.6551554698226734 1.5288445158419737 - 1.4141414141414141 1.3801655474948953 1.6997765020687188 - 1.6161616161616161 1.1386991324915792 1.8894455179547536 - 1.8181818181818181 0.93452394063678168 2.0868798701189748 - 2.0202020202020203 0.76883509136062878 2.2831835057149235 - 2.2222222222222223 0.64049756544921577 2.4717945496694602 - 2.4242424242424243 0.54582998647302927 2.6491801934849879 - 2.6262626262626263 0.47923804898358729 2.8146729922885356 - 2.8282828282828283 0.43448623094088068 2.9694686333537303 - 3.0303030303030303 0.40582998221767225 3.1154947020783408 - 3.2323232323232323 0.38862631305559836 3.2546318501945897 - 3.4343434343434343 0.37944413244255548 3.3883642968367993 - 3.6363636363636362 0.37589405690585731 3.5177175370662574 - 3.8383838383838382 0.37638816454364082 3.6433135712446219 - 4.0404040404040407 0.37991546615719302 3.7654567349270782 - 4.2424242424242422 0.38587045159306321 3.8842067676369543 - 4.4444444444444446 0.39393734842055028 3.9994245192929223 - 4.6464646464646462 0.40401657963564586 4.1107970921412615 - 4.8484848484848486 0.41619382893307366 4.2178262889240079 - 5.0505050505050502 0.43073595063159331 4.3197979141243916 - 5.2525252525252526 0.44812576630096862 4.4157055096252806 - 5.4545454545454541 0.46913981953040124 4.5041256293946068 - 5.6565656565656566 0.49501076529971422 4.5829800566007810 - 5.8585858585858590 0.52774916258514326 4.6490893419504102 - 6.0606060606060606 0.57079477784972987 4.6973125301748198 - 6.2626262626262630 0.63050350250814691 4.7186371751684977 - 6.4646464646464645 0.71988449040201608 4.6955426522082018 - 6.6666666666666670 0.86968691097161099 4.5886613618807219 - 6.8686868686868685 1.1685706770050810 4.2898652840507063 - 7.0707070707070709 1.9283238624518928 3.4323083638716234 - 7.2727272727272725 3.4739424266532057 1.5443774908359309 - 7.4747474747474749 3.6156622365325570 0.86126034682410690 - 7.6767676767676765 3.0893285893108224 0.91109911890430328 - 7.8787878787878789 2.5824091029827643 1.0482685855716802 - 8.0808080808080813 2.1443777121787724 1.2120227059149227 - 8.2828282828282820 1.7697147568959450 1.3943481014446395 - 8.4848484848484844 1.4503912908442944 1.5913176363225987 - 8.6868686868686869 1.1808973224083268 1.7978634657477370 - 8.8888888888888893 0.95820434862210935 2.0072845409532536 - 9.0909090909090917 0.78027041788365581 2.2123753298973790 - 9.2929292929292924 0.64414161929497971 2.4073173656722280 - 9.4949494949494948 0.54485918823617863 2.5890901458343887 - 9.6969696969696972 0.47579479033867622 2.7575304068440558 - 9.8989898989898997 0.42987872902400526 2.9143195370050679 - 10.101010101010100 0.40077868671288164 3.0617735612576609 - 10.303030303030303 0.38347341460130357 3.2020467883080310 - 10.505050505050505 0.37431238709501730 3.3367975669544059 - 10.707070707070708 0.37079242852503225 3.4671539002851515 - 10.909090909090908 0.37126507364401884 3.5938030419551850 - 11.111111111111111 0.37469070551565176 3.7170914801702328 - 11.313131313131313 0.38044518892229101 3.8371158243297301 - 11.515151515151516 0.38819638162477937 3.9537755040428388 - 11.717171717171718 0.39782585999772907 4.0668002770121863 - 11.919191919191919 0.40938746672019916 4.1757548961032134 - 12.121212121212121 0.42310068726159888 4.2800097184239387 - 12.323232323232324 0.43937066111732542 4.3786876642036985 - 12.525252525252526 0.45885057315802091 4.4705573015745195 - 12.727272727272727 0.48256241252345816 4.5538539435818004 - 12.929292929292929 0.51213939269131681 4.6259372029288572 - 13.131313131313131 0.55030347871379681 4.6826500462048699 - 13.333333333333334 0.60190272142903367 4.7169701251871379 - 13.535353535353535 0.67640159980253611 4.7158747002446546 - 13.737373737373737 0.79472320674162111 4.6520034788467930 - 13.939393939393939 1.0119749404433043 4.4567948585946260 - 14.141414141414142 1.5095312791003548 3.9143914746575619 - 14.343434343434344 2.8225321990354586 2.3841749020648453 - 14.545454545454545 3.7481381534620617 0.96112207892930535 - 14.747474747474747 3.3218265036555348 0.86862083704991366 - 14.949494949494950 2.7904932667385376 0.98537464458288837 - 15.151515151515152 2.3229635059729983 1.1396150506772385 - 15.353535353535353 1.9223132692756113 1.3145393802639307 - 15.555555555555555 1.5802165899411555 1.5057947090731012 - 15.757575757575758 1.2899587327013933 1.7090002440466507 - 15.959595959595960 1.0475651849438885 1.9180916586741106 - 16.161616161616163 0.85081575529077458 2.1258676910815164 - 16.363636363636363 0.69734413658792593 2.3256879814402112 - 16.565656565656564 0.58308589778044639 2.5132480548271525 - 16.767676767676768 0.50200858339959065 2.6872450059000279 - 16.969696969696969 0.44705821749630464 2.8487327984749697 - 17.171717171717173 0.41148360626831598 2.9998828400036008 - 17.373737373737374 0.38967114474781089 3.1429839739727425 - 17.575757575757574 0.37740468865260540 3.2799215912053152 - 17.777777777777779 0.37172394237646628 3.4120450807555018 - 17.979797979797979 0.37065177405380056 3.5402151815833696 - 18.181818181818183 0.37291573018018059 3.6649116005999534 - 18.383838383838384 0.37773259534751530 3.7863295806542410 - 18.585858585858585 0.38465850452108946 3.9044466320941358 - 18.787878787878789 0.39348957923200001 4.0190654624651199 - 18.989898989898990 0.40421007951934346 4.1298207698584255 - 19.191919191919190 0.41696810088968250 4.2361704331505727 - 19.393939393939394 0.43208411483903786 4.3373537507479538 - 19.595959595959595 0.45009453637804475 4.4323079947361164 - 19.797979797979799 0.47184398964825353 4.5195244763911635 - 20.000000000000000 0.49866591845763342 4.5967890211579494 + 0.20202020202020202 1.9228273596482419 2.4337073227227015 + 0.40404040404040403 2.3295379939372070 1.7950504631676178 + 0.60606060606060608 2.4260084390175374 1.4138421392890184 + 0.80808080808080807 2.2389134768006693 1.3287817289785369 + 1.0101010101010102 1.9501499547207590 1.3960246302359856 + 1.2121212121212122 1.6531650630809118 1.5316301779428794 + 1.4141414141414141 1.3792870041080729 1.7020748762093159 + 1.6161616161616161 1.1386899395675674 1.8913010133772596 + 1.8181818181818181 0.93509048492949876 2.0883310199995488 + 2.0202020202020203 0.76973204188308597 2.2843176596991053 + 2.2222222222222223 0.64162988561965573 2.4726833763191598 + 2.4242424242424243 0.54704122772270458 2.6499082925367090 + 2.6262626262626263 0.48035490068191056 2.8153424283053656 + 2.8282828282828283 0.43549086882672539 2.9701091612554205 + 3.0303030303030303 0.40671728136167717 3.1161198232215281 + 3.2323232323232323 0.38938090509460166 3.2552573696000682 + 3.4343434343434343 0.38009582896557559 3.3889745277175156 + 3.6363636363636362 0.37647119134566065 3.5182908610361623 + 3.8383838383838382 0.37691455991106404 3.6438243083256565 + 4.0404040404040407 0.38038154649432943 3.7659160989621996 + 4.2424242424242422 0.38629189791005092 3.8846010921751279 + 4.4444444444444446 0.39435323461674171 3.9996906739244662 + 4.6464646464646462 0.40441856966793516 4.1109530323192081 + 4.8484848484848486 0.41659226649230330 4.2178608105899231 + 5.0505050505050502 0.43114054403594065 4.3197021316738917 + 5.2525252525252526 0.44855609647799038 4.4154496317592651 + 5.4545454545454541 0.46965255590129362 4.5036014167113407 + 5.6565656565656566 0.49562770237026382 4.5821544963722030 + 5.8585858585858590 0.52850631874549547 4.6479120844663200 + 6.0606060606060606 0.57180653469716725 4.6956095368332349 + 6.2626262626262630 0.63189644880900575 4.7162407311256080 + 6.4646464646464645 0.72218995702309585 4.6917232898748678 + 6.6666666666666670 0.87397090869856675 4.5820707979270905 + 6.8686868686868685 1.1782095641881043 4.2764311463919649 + 7.0707070707070709 1.9549909248615882 3.3982085218696105 + 7.2727272727272725 3.4927584010647701 1.5112511142461154 + 7.4747474747474749 3.6014521596267524 0.86175709858044058 + 7.6767676767676765 3.0753099008548137 0.91470615338806938 + 7.8787878787878789 2.5711400845992336 1.0523043530061773 + 8.0808080808080813 2.1354241731711872 1.2162449988084660 + 8.2828282828282820 1.7625746350906142 1.3987041080686322 + 8.4848484848484844 1.4447028444892016 1.5957349627925383 + 8.6868686868686869 1.1764848714523968 1.8022394250479519 + 8.8888888888888893 0.95491482859665922 2.0114872074773231 + 9.0909090909090917 0.77798931399884763 2.2162934968831585 + 9.2929292929292924 0.64268702380879605 2.4109049490213335 + 9.4949494949494948 0.54410802353831678 2.5923331608443352 + 9.6969696969696972 0.47557051706286352 2.7604697535954306 + 9.8989898989898997 0.42992111793267701 2.9170678083836674 + 10.101010101010100 0.40102347313261461 3.0643541338256242 + 10.303030303030303 0.38382665352188139 3.2045022564565206 + 10.505050505050505 0.37473022610167511 3.3391412357774541 + 10.707070707070708 0.37122133500797094 3.4694115861465273 + 10.909090909090908 0.37170000456577573 3.5959654585400180 + 11.111111111111111 0.37511538344853684 3.7191666499367577 + 11.313131313131313 0.38088271752884140 3.8390613591672160 + 11.515151515151516 0.38865561609292282 3.9555626954474099 + 11.717171717171718 0.39831554438992395 4.0684015783227530 + 11.919191919191919 0.40991583595477665 4.1771484653327589 + 12.121212121212121 0.42366589320023684 4.2812044324572511 + 12.323232323232324 0.43999689990295954 4.3796417540125230 + 12.525252525252526 0.45956884285617883 4.4712180593723572 + 12.727272727272727 0.48340935262992762 4.5541642541999847 + 12.929292929292929 0.51320529648621338 4.6257472085247393 + 13.131313131313131 0.55170034140558755 4.6817919313649616 + 13.333333333333334 0.60382639488504730 4.7151679282352639 + 13.535353535353535 0.67943350666239688 4.7123168017210988 + 13.737373737373737 0.79979418320966067 4.6455117030896727 + 13.939393939393939 1.0225440882724239 4.4431032690230703 + 14.141414141414142 1.5368177664562332 3.8802774959603572 + 14.343434343434344 2.8796173081158947 2.3109971592952023 + 14.545454545454545 3.7367373854770385 0.95102344785598902 + 14.747474747474747 3.3002007207544355 0.87294377605358464 + 14.949494949494950 2.7720843934008390 0.99101409440307164 + 15.151515151515152 2.3078344664784458 1.1457541519462686 + 15.353535353535353 1.9098377341722368 1.3210820748744150 + 15.555555555555555 1.5699491230246023 1.5126501470035814 + 15.757575757575758 1.2816621331869149 1.7160056688348584 + 15.959595959595960 1.0411154988544584 1.9250118226341943 + 16.161616161616163 0.84599939080373476 2.1324740169202734 + 16.363636363636363 0.69393083686771306 2.3318429273305119 + 16.565656565656564 0.58093651523930878 2.5188685057586175 + 16.767676767676768 0.50078269028644173 2.6923926868888022 + 16.969696969696969 0.44647986346531887 2.8534987221981112 + 17.171717171717173 0.41131481489406863 3.0043689436169685 + 17.373737373737374 0.38976361314897068 3.1472572110695851 + 17.575757575757574 0.37767878047252879 3.2840112580710215 + 17.777777777777779 0.37208547244484813 3.4159969992226542 + 17.979797979797979 0.37109152272324525 3.5440118041481710 + 18.181818181818183 0.37339271637681176 3.6685692271543546 + 18.383838383838384 0.37825480172285553 3.7898132418446799 + 18.585858585858585 0.38522151208206001 3.9077390349872942 + 18.787878787878789 0.39408978816425594 4.0221592912290030 + 18.989898989898990 0.40484676902509076 4.1327168843298816 + 19.191919191919190 0.41765931977550447 4.2388375058919161 + 19.393939393939394 0.43286367021270983 4.3397272319250382 + 19.595959595959595 0.45099500426075234 4.4343300857781136 + 19.797979797979799 0.47292513140904857 4.5210892767632815 + 20.000000000000000 0.50002673252461749 4.5977273746884384 diff --git a/plotter.py b/plotter.py index 56dd602..4d08e77 100644 --- a/plotter.py +++ b/plotter.py @@ -6,6 +6,7 @@ from mpl_toolkits.mplot3d import Axes3D filename = 'data.out' +# filename = 'raw.out' fig = plt.figure() diff --git a/src/main.f90 b/src/main.f90 index 6dc6421..5ede2e4 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -13,6 +13,9 @@ program main real(wp), dimension(n) :: y0 real(wp), dimension(num_t,n) :: y + open(unit=22, file='raw.out', iostat=ios, status="replace", action="write") + if ( ios /= 0 ) stop "Error opening file 22" + ! Set time vector `t` t0 = 0.0_wp ! tend = 4_wp*pi @@ -20,29 +23,32 @@ program main call linspace(t0, tend, t) ! Set y vector `y` using y0 - ! y0 = [0._wp, 1._wp] + ! y0 = [0._wp, 0.1_wp] y0 = [1.5_wp, 3._wp] ! y0 = [0.1_wp, 0.0_wp, 0.0_wp] ! y0 = [0._wp, 1._wp, 0._wp] y(1,:) = y0 + write(22,*) t0, y0 + + call rk_wrapper(n, num_t, t, y, mysub) open(unit=21, file='data.out', iostat=ios, status="replace", action="write") if ( ios /= 0 ) stop "Error opening file 21" - open(unit=22, file='raw.out', iostat=ios, status="replace", action="write") - if ( ios /= 0 ) stop "Error opening file 22" do ii = 1, num_t write(21,*) t(ii), y(ii, :) end do + print*, t(num_t), y(num_t, :) + close(unit=21, iostat=ios) if ( ios /= 0 ) stop "Error closing file unit 21" - close(unit=22, iostat=ios, status="delete") + close(unit=22, iostat=ios) if ( ios /= 0 ) stop "Error closing file unit 22" diff --git a/src/misc.f90 b/src/misc.f90 index 7f8a91e..86d0b6f 100644 --- a/src/misc.f90 +++ b/src/misc.f90 @@ -52,12 +52,12 @@ subroutine vanderpol(n, t, y, dy) real(wp), intent(in), dimension(n) :: y real(wp), intent(out), dimension(n) :: dy - real(wp), parameter :: mu = 30.53_wp + real(wp), parameter :: mu = 0.2_wp real(wp), parameter :: A = 3.2_wp real(wp), parameter :: omega = 2._wp*pi/10_wp dy(1) = y(2) - dy(2) = mu * (1._wp-y(1)*y(1)) * y(2) - y(1) + A*dsin(t*omega) + dy(2) = mu * (1._wp-y(1)*y(1)) * y(2) - y(1) !+ A*dsin(t*omega) return end subroutine vanderpol diff --git a/src/rk_constants.f90 b/src/rk_constants.f90 index 92e1a71..d6cdb32 100644 --- a/src/rk_constants.f90 +++ b/src/rk_constants.f90 @@ -6,19 +6,19 @@ module rk_constants public :: midpoint, heun, ralston public :: heun_euler, fehlbery_rk12, bogacki_shampine public :: fehlbery_rk45, cash_karp_rk45, dormand_prince_rk45 + public :: fehlbery_rk78 private :: two_stage contains - subroutine midpoint(a, b, c, m, p) - integer, intent(out) :: m, p + subroutine midpoint(a, b, c, m) + integer, intent(out) :: m real(wp) :: alpha real(wp), intent(out), dimension(:), allocatable :: c, b real(wp), intent(out), dimension(:,:), allocatable :: a - p = 2 m = 2 allocate(a(m,m), b(m), c(m)) @@ -28,13 +28,12 @@ subroutine midpoint(a, b, c, m, p) return end subroutine midpoint - subroutine heun(a, b, c, m, p) - integer, intent(out) :: m, p + subroutine heun(a, b, c, m) + integer, intent(out) :: m real(wp) :: alpha real(wp), intent(out), dimension(:), allocatable :: c, b real(wp), intent(out), dimension(:,:), allocatable :: a - p = 2 m = 2 allocate(a(m,m), b(m), c(m)) @@ -44,13 +43,12 @@ subroutine heun(a, b, c, m, p) return end subroutine heun - subroutine ralston(a, b, c, m, p) - integer, intent(out) :: m, p + subroutine ralston(a, b, c, m) + integer, intent(out) :: m real(wp) :: alpha real(wp), intent(out), dimension(:), allocatable :: c, b real(wp), intent(out), dimension(:,:), allocatable :: a - p = 2 m = 2 allocate(a(m,m), b(m), c(m)) @@ -77,12 +75,11 @@ end subroutine two_stage - subroutine classic_rk4(a, b, c, m, p) - integer, intent(out) :: m, p + subroutine classic_rk4(a, b, c, m) + integer, intent(out) :: m real(wp), intent(out), dimension(:), allocatable :: c, b real(wp), intent(out), dimension(:,:), allocatable :: a - p = 4 m = 4 allocate(a(m,m), b(m), c(m)) @@ -97,12 +94,11 @@ subroutine classic_rk4(a, b, c, m, p) return end subroutine classic_rk4 - subroutine three_eighths_rk4(a, b, c, m, p) - integer, intent(out) :: m, p + subroutine three_eighths_rk4(a, b, c, m) + integer, intent(out) :: m real(wp), intent(out), dimension(:), allocatable :: c, b real(wp), intent(out), dimension(:,:), allocatable :: a - p = 4 m = 4 allocate(a(m,m), b(m), c(m)) @@ -131,7 +127,7 @@ subroutine heun_euler(a, b, bstar, c, m, p) ! There is already a routine to calculate the heun coefficients - just ! reuse it for the first row of `b`, then set bstar - call heun(a, b, c, m, p) + call heun(a, b, c, m) bstar = [1._wp, 0._wp] return @@ -262,4 +258,51 @@ subroutine dormand_prince_rk45(a, b, bstar, c, m, p) return end subroutine dormand_prince_rk45 + subroutine fehlbery_rk78(a, b, bstar, c, m, p) + integer, intent(out) :: m, p + real(wp), intent(out), dimension(:), allocatable :: c, b, bstar + real(wp), intent(out), dimension(:,:), allocatable :: a + + p = 7 + m = 13 + allocate(a(m,m), b(m), bstar(m), c(m)) + + c = [0._wp, 2._wp/27._wp, 1._wp/9._wp, 1._wp/6._wp, 5._wp/12._wp, & + & 0.5_wp, 5._wp/6._wp, 1._wp/6._wp, 2._wp/3._wp, & + & 1._wp/3._wp, 1._wp, 0._wp, 1._wp] + b = [41._wp/840._wp, 0._wp, 0._wp, 0._wp, 0._wp, 34._wp/105._wp, & + & 9._wp/35._wp, 9._wp/35._wp, 9._wp/280._wp, 9._wp/280._wp, & + & 41._wp/840._wp, 0._wp, 0._wp] + bstar = [0._wp, 0._wp, 0._wp, 0._wp, 0._wp, 34._wp/105._wp, & + & 9._wp/35._wp, 9._wp/35._wp, 9._wp/280._wp, & + & 9._wp/280._wp, 0._wp, 41._wp/840._wp, 41._wp/840._wp] + + a = 0._wp + a(2,1) = 2._wp/27._wp + a(3, 1:2) = [1._wp/36._wp, 1._wp/12._wp] + a(4, 1:3) = [1._wp/24._wp, 0._wp, 1._wp/8._wp] + a(5, 1:4) = [5._wp/12._wp, 0._wp, -25._wp/16._wp, 25._wp/16._wp] + a(6, 1:5) = [1._wp/20._wp, 0._wp, 0._wp, 1._wp/4._wp, 1._wp/5._wp] + a(7, 1:6) = [-25._wp/108._wp, 0._wp, 0._wp, 125._wp/108._wp, & + & -65._wp/27._wp, 125._wp/54._wp] + a(8, 1:7) = [31._wp/300._wp, 0._wp, 0._wp, 0._wp, 61._wp/225._wp, & + & -2._wp/9._wp, 13._wp/900._wp] + a(9, 1:8) = [2._wp, 0._wp, 0._wp, -53._wp/6._wp, 704._wp/45._wp, & + & -107._wp/9._wp, 67._wp/90._wp, 3._wp] + a(10, 1:9) = [-91._wp/108._wp, 0._wp, 0._wp, 23._wp/108._wp, & + & -976._wp/135._wp, 311._wp/54._wp, -19._wp/60._wp, & + & 17._wp/6._wp, -1._wp/12._wp] + a(11, 1:10) = [2383._wp/4100._wp, 0._wp, 0._wp, -341._wp/164._wp, & + & 4496._wp/1025._wp, -301._wp/82._wp, 2133._wp/4100._wp, & + & 45._wp/82._wp, 45._wp/164._wp, 18._wp/41._wp] + a(12, 1:11) = [3._wp/205._wp, 0._wp, 0._wp, 0._wp, 0._wp, -6._wp/41._wp, & + & -3._wp/205._wp, -3._wp/41._wp, 3._wp/41._wp, 6._wp/41._wp, & + & 0._wp] + a(13, 1:12) = [-1777._wp/4100._wp, 0._wp, 0._wp, -341._wp/164._wp, & + & 4496._wp/1025._wp, -289._wp/82._wp, 2193._wp/4100._wp, & + & 51._wp/82._wp, 33._wp/164, 19._wp/41._wp, 0._wp, 1._wp] + + return + end subroutine fehlbery_rk78 + end module rk_constants diff --git a/src/runge_kutta.f90 b/src/runge_kutta.f90 index a0e6035..a1478aa 100644 --- a/src/runge_kutta.f90 +++ b/src/runge_kutta.f90 @@ -4,6 +4,7 @@ module runge_kutta use rk_constants, only: midpoint, heun, ralston use rk_constants, only: heun_euler, fehlbery_rk12, bogacki_shampine use rk_constants, only: fehlbery_rk45, cash_karp_rk45, dormand_prince_rk45 + use rk_constants, only: fehlbery_rk78 implicit none contains @@ -33,8 +34,12 @@ subroutine dysub(n, t, y, dy) real(wp), dimension(:,:), allocatable :: a ! Get the tableau coefficients associated with a runge kutta implementation - call dormand_prince_rk45(a, b, bstar, c, m, p) + ! call heun_euler(a, b, bstar, c, m, p) + ! call dormand_prince_rk45(a, b, bstar, c, m, p) + call fehlbery_rk78(a, b, bstar, c, m, p) + ! call heun(a, b, c, m, p) + ! call three_eighths_rk4(a, b, c, m, p) ! Initial guess for dt is just the difference between t(1:2), this way new ! versions of dt will be saved in each call to rk_adaptive @@ -45,14 +50,17 @@ subroutine dysub(n, t, y, dy) tspan = t(ii:ii+1) call rk_adaptive(n, tspan, dt, yy, dysub, a, b, bstar, c, m, p) - ! call rk_explicit(n, t(ii), t(ii+1)-t(ii), yy, dysub, a, b, c, m, p) + ! call rk_explicit(n, t(ii), t(ii+1)-t(ii), yy, dysub, a, b, c, m) y(ii+1, :) = yy ! if (ii == 1) stop end do - deallocate(a, b, bstar, c) + if (allocated(a)) deallocate(a) + if (allocated(b)) deallocate(b) + if (allocated(bstar)) deallocate(bstar) + if (allocated(c)) deallocate(c) return end subroutine rk_wrapper @@ -77,14 +85,15 @@ subroutine dysub(n, t, y, dy) end interface ! Local arguments - integer :: ii, jj, kk real(wp) :: t, error - real(wp), parameter :: eps1 = sqrt(epsilon(1._wp)) + real(wp), parameter :: eps = sqrt(epsilon(1._wp)) + real(wp), parameter :: eps_abs = 1d-5 ! sqrt(epsilon(1._wp)) + real(wp), parameter :: eps_rel = 1d-5 real(wp), dimension(n) :: dummy_y, y_star integer :: eval, num_dt - integer, parameter :: max_eval = 100 - real(wp) :: s, new_s + integer, parameter :: max_eval = 500 + real(wp) :: s, gamma, sc, maxy logical :: last num_dt = 1 @@ -112,38 +121,39 @@ subroutine dysub(n, t, y, dy) ! error = maxval(abs(dummy_y - y_star)) error = norm2(dummy_y - y_star) - - s = eps1*dt / & - & (2._wp * (tspan(2) - tspan(1)) * error) + maxy = maxval([ maxval(abs(dummy_y)), maxval(abs(y_star)) ]) + + sc = eps_abs + maxy * eps_rel + gamma = 0.25_wp ** (1._wp / real(p,wp)) + ! gamma = 0.9_wp + ! s = gamma * (eps/error) ** ( 1._wp/real(p,wp) ) + s = gamma * (sc/error) ** ( 1._wp/real(p,wp) ) + ! s = eps*dt / & + ! & (2._wp * (tspan(2) - tspan(1)) * error) ! print*, 'error = ', error + ! print*, 'maxy = ', maxy + ! print*, 'sc = ', sc ! print*, 's = ', s ! print*, 'sqrt(s) = ', sqrt(s) ! print*, - - new_s = 0.9_wp * (eps1/error) ** (1._wp/real(p, wp)) - ! print*, 'new_s = ', new_s ! stop ! Adjust dt based on estimation of error - if ( new_s >= 2._wp ) then + if ( s >= 2._wp ) then t = t + dt y = y_star dt = dt * 2._wp ! print*, 'Adjusted dt with factor of 2' exit - else if ( new_s < 2._wp .and. new_s >= 1._wp ) then + else if ( s < 2._wp .and. s >= 1._wp ) then t = t + dt y = y_star ! dt = dt * s ! print 121, 'Adjusted dt with factor of ', s - ! dt = dt * new_s - ! print 121, 'Adjusted dt with factor of ', new_s exit - else if ( new_s <= 1._wp ) then - ! dt = dt * s + else if ( s <= 1._wp ) then + dt = dt * s ! print 121, 'Adjusted dt with factor of ', s - dt = dt * new_s - ! print 121, 'Adjusted dt with factor of ', new_s end if if ( eval == max_eval ) then @@ -152,6 +162,8 @@ subroutine dysub(n, t, y, dy) end do + write(22,*) t, y + ! num_dt = num_dt + 1 ! print 120, t, tspan, dt, y