diff --git a/sources/driver.F90 b/sources/driver.F90 index 8809f08..8be0daa 100644 --- a/sources/driver.F90 +++ b/sources/driver.F90 @@ -48,7 +48,7 @@ program amun use equations , only : nv, nf use evolution , only : initialize_evolution, finalize_evolution use evolution , only : print_evolution - use evolution , only : advance, new_time_step + use evolution , only : advance, initialize_time_step, new_time_step use evolution , only : registers, step, time, dt, dtp, errtol use forcing , only : initialize_forcing, finalize_forcing use forcing , only : print_forcing @@ -62,7 +62,7 @@ program amun use io , only : restart_snapshot_number, restart_from_snapshot use io , only : read_snapshot_parameter use io , only : read_restart_snapshot, write_restart_snapshot - use io , only : write_snapshot, next_tout, precise_snapshots + use io , only : write_snapshot, update_dtp use mesh , only : initialize_mesh, finalize_mesh use mesh , only : generate_mesh, store_mesh_stats use mpitools , only : initialize_mpitools, finalize_mpitools @@ -326,10 +326,6 @@ program amun ! trun = trun - tsav / 6.0d+01 -! initialize dtp -! - dtp = 2.0d+00 * tmax - ! get integral calculation interval ! call get_parameter("ndat", ndat) @@ -553,9 +549,9 @@ program amun ! call boundary_variables(0.0d+00, 0.0d+00) -! calculate new timestep +! estimate the initial timestep ! - call new_time_step() + call initialize_time_step() end if @@ -621,10 +617,6 @@ program amun ! do while((nsteps <= nmax) .and. (time < tmax) .and. proceed) -! get the next snapshot time -! - if (precise_snapshots) dtp = next_tout() - time - ! performe one step evolution ! call advance(status) @@ -643,6 +635,14 @@ program amun step = step + 1 nsteps = nsteps + 1 +! update time step for precise snapshots +! + call update_dtp() + +! estimate the next time step +! + call new_time_step() + ! get current time in seconds ! tm_curr = get_timer_total() diff --git a/sources/evolution.F90 b/sources/evolution.F90 index a9b4d3c..86abe1e 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -68,8 +68,11 @@ module evolution ! evolution parameters ! logical , save :: error_control = .false. + logical , save :: fsal = .false. + logical , save :: fsal_update = .true. character(len=255), save :: name_int = "" character(len=255), save :: name_norm = "" + integer , save :: order = 2 integer , save :: stages = 2 integer , save :: registers = 1 real(kind=8) , save :: cfl = 5.0d-01 @@ -83,10 +86,9 @@ module evolution integer , save :: step = 0 real(kind=8) , save :: time = 0.0d+00 real(kind=8) , save :: dt = 1.0d+00 - real(kind=8) , save :: dtn = 1.0d+00 real(kind=8) , save :: dth = 1.0d+00 real(kind=8) , save :: dte = 0.0d+00 - real(kind=8) , save :: dtp = 1.0d+00 ! dt for the precise snapshots + real(kind=8) , save :: dtp = huge(dt) ! dt for precise snapshots ! the absolute and relative tolerances, limiting factors, the maximum error, ! the maximum number of passes for the adaptive step, @@ -103,6 +105,11 @@ module evolution integer , save :: niterations = 0 integer , save :: nrejections = 0 +! coefficients of Runge-Kutta methods +! + real(kind=8), dimension(:) , allocatable, save :: c, bt, dl, bh + real(kind=8), dimension(:,:), allocatable, save :: gm + ! errors in three recent steps ! real(kind=8), dimension(3), save :: betas, errs @@ -118,11 +125,11 @@ module evolution ! declare public subroutines ! public :: initialize_evolution, finalize_evolution, print_evolution - public :: advance, new_time_step + public :: initialize_time_step, new_time_step, advance ! declare public variables ! - public :: step, time, dt, dtn, dth, dte, dtp, cfl, glm_alpha, registers + public :: step, time, dt, dth, dte, dtp, cfl, glm_alpha, registers public :: atol, rtol, mrej, niterations, nrejections, errs, errtol !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - @@ -222,12 +229,14 @@ module evolution case ("euler", "EULER", "rk1", "RK1") evolve => evolve_euler + order = 1 registers = 1 name_int = "1st order Euler" case ("rk2", "RK2") evolve => evolve_rk2 + order = 2 registers = 2 name_int = "2nd order Runge-Kutta" @@ -235,6 +244,7 @@ module evolution error_control = .true. evolve => evolve_ssprk2_m + order = 2 registers = 3 stages = max(2, min(9, stages)) cfl = (stages - 1) * cfl @@ -243,12 +253,14 @@ module evolution case ("rk3", "RK3") evolve => evolve_rk3 + order = 3 registers = 2 name_int = "3rd order Runge-Kutta" case ("rk3.4", "RK3.4", "ssprk(4,3)", "SSPRK(4,3)") evolve => evolve_ssprk34 + order = 3 registers = 2 cfl = 2.0d+00 * cfl name_int = "3rd order SSPRK(4,3)" @@ -256,6 +268,7 @@ module evolution case ("rk3.5", "RK3.5", "ssprk(5,3)", "SSPRK(5,3)") evolve => evolve_ssprk35 + order = 3 registers = 2 cfl = 2.65062919143939d+00 * cfl name_int = "3rd order SSPRK(5,3)" @@ -264,6 +277,7 @@ module evolution error_control = .true. evolve => evolve_ssp324 + order = 3 registers = 3 stages = 4 cfl = 2.0d+00 * cfl @@ -281,6 +295,7 @@ module evolution error_control = .true. evolve => evolve_ssprk32m + order = 3 registers = 4 n = 2 do while(n**2 <= stages) @@ -312,10 +327,524 @@ module evolution case ("rk4.10", "RK4.10", "ssprk(10,4)", "SSPRK(10,4)") evolve => evolve_ssprk4_10 + order = 4 registers = 2 cfl = 6.0d+00 * cfl name_int = "Optimal 4th order SSPRK(10,4)" + case ("RK3(2)5", "rk3(2)5") + +! References: +! +! [1] Ranocha, H., Dalcin, L., Parsani, M., Ketcheson, D. I. +! "Coefficients of Optimized Runge-Kutta Methods with Automatic +! Step Size Control for Compressible Computational Fluid Dynamics", +! 2021, DOI:10.5281/zenodo.4671927, +! https://github.com/ranocha/Optimized-RK-CFD +! + error_control = .true. + evolve => evolve_3sstarplus + registers = 4 + order = 3 + stages = 5 + cfl = 2.5d+00 * cfl + name_int = "Optimized 3ʳᵈ-order 5-step embedded RK3(2)5" + + allocate(c(stages), dl(stages), bt(stages), bh(stages+1), gm(3,stages)) + c (:) = [ 0.000000000000000000000000000000000000d+00, & + 2.300285062878154351930669430512780706d-01, & + 4.050049049262914975700372321130661410d-01, & + 8.947823877926760224705450466361360720d-01, & + 7.235108137218888081489570284485201518d-01 ] + dl(:) = [ 1.000000000000000000000000000000000000d+00, & + 3.407687209321455242558804921815861422d-01, & + 3.414399280584625023244387687873774697d-01, & + 7.229302732875589702087936723400941329d-01, & + 0.000000000000000000000000000000000000d+00 ] + bt(:) = [ 2.300285062878154351930669430512780706d-01, & + 3.021457892454169700189445968126242994d-01, & + 8.025601039472704213300183888573974531d-01, & + 4.362158997637629844305216319994356355d-01, & + 1.129268494470295369172265188216779157d-01 ] + bh(:) = [ 1.046363371354093758897668305991705199d-01, & + 9.520431574956758809511173383346476348d-02, & + 4.482446645568668405072421350300379357d-01, & + 2.449030295461310135957132640369862245d-01, & + 1.070116530120251819121660365003405564d-01, & + 0.000000000000000000000000000000000000d+00 ] + gm(1,:) = [ 0.000000000000000000000000000000000000d+00, & + 2.587669070352079020144955303389306026d-01, & + -1.324366873994502973977035353758550057d-01, & + 5.055601231460399101814291350373559483d-02, & + 5.670552807902877312521811889846000976d-01 ] + gm(2,:) = [ 1.000000000000000000000000000000000000d+00, & + 5.528418745102160639901976698795928733d-01, & + 6.731844400389673824374042790213570079d-01, & + 2.803103804507635075215805236096803381d-01, & + 5.521508873507393276457754945308880998d-01 ] + gm(3,:) = [ 0.000000000000000000000000000000000000d+00, & + 0.000000000000000000000000000000000000d+00, & + 0.000000000000000000000000000000000000d+00, & + 2.752585813446636957256614568573008811d-01, & + -8.950548709279785077579454232514633376d-01 ] + + betas(:) = [ 6.4d-01, -3.1d-01, 4.0d-02 ] + + call get_parameter("beta(1)", betas(1)) + call get_parameter("beta(2)", betas(2)) + call get_parameter("beta(3)", betas(3)) + + betas(:) = - betas(:) / 3.0d+00 + + case ("RK3(2)5F", "rk3(2)5f") + +! References: +! +! [1] Ranocha, H., Dalcin, L., Parsani, M., Ketcheson, D. I. +! "Coefficients of Optimized Runge-Kutta Methods with Automatic +! Step Size Control for Compressible Computational Fluid Dynamics", +! 2021, DOI:10.5281/zenodo.4671927, +! https://github.com/ranocha/Optimized-RK-CFD +! + error_control = .true. + fsal = .true. + evolve => evolve_3sstarplus + registers = 4 + order = 3 + stages = 5 + cfl = 2.5d+00 * cfl + name_int = "Optimized 3ʳᵈ-order 5-step embedded RK3(2)5 FSAL" + + allocate(c(stages), dl(stages), bt(stages), bh(stages+1), gm(3,stages)) + c (:) = [ 0.000000000000000000000000000000000000d+00, & + 2.300298624518076223899418286314123354d-01, & + 4.050046072094990912268498160116125481d-01, & + 8.947822893693433545220710894560512805d-01, & + 7.235136928826589010272834603680114769d-01 ] + dl(:) = [ 1.000000000000000000000000000000000000d+00, & + 3.407655879334525365094815965895763636d-01, & + 3.414382655003386206551709871126405331d-01, & + 7.229275366787987419692007421895451953d-01, & + 0.000000000000000000000000000000000000d+00 ] + bt(:) = [ 2.300298624518076223899418286314123354d-01, & + 3.021434166948288809034402119555380003d-01, & + 8.025606185416310937583009085873554681d-01, & + 4.362158943603440930655148245148766471d-01, & + 1.129272530455059129782111662594436580d-01 ] + bh(:) = [ 9.484166705035703392326247283838082847d-02, & + 1.726371339430353766966762629176676070d-01, & + 3.998243189084371024483169698618455770d-01, & + 1.718016807580178450618829007973835152d-01, & + 5.881914422155740300718268359027168467d-02, & + 1.020760551185952388626787099944507877d-01 ] + gm(1,:) = [ 0.000000000000000000000000000000000000d+00, & + 2.587771979725733308135192812685323706d-01, & + -1.324380360140723382965420909764953437d-01, & + 5.056033948190826045833606441415585735d-02, & + 5.670532000739313812633197158607642990d-01 ] + gm(2,:) = [ 1.000000000000000000000000000000000000d+00, & + 5.528354909301389892439698870483746541d-01, & + 6.731871608203061824849561782794643600d-01, & + 2.803103963297672407841316576323901761d-01, & + 5.521525447020610386070346724931300367d-01 ] + gm(3,:) = [ 0.000000000000000000000000000000000000d+00, & + 0.000000000000000000000000000000000000d+00, & + 0.000000000000000000000000000000000000d+00, & + 2.752563273304676380891217287572780582d-01, & + -8.950526174674033822276061734289327568d-01 ] + + betas(:) = [ 7.0d-01, -2.3d-01, 0.0d+00 ] + + call get_parameter("beta(1)", betas(1)) + call get_parameter("beta(2)", betas(2)) + call get_parameter("beta(3)", betas(3)) + + betas(:) = - betas(:) / 3.0d+00 + + case ("RK4(3)9", "rk4(3)9") + +! References: +! +! [1] Ranocha, H., Dalcin, L., Parsani, M., Ketcheson, D. I. +! "Coefficients of Optimized Runge-Kutta Methods with Automatic +! Step Size Control for Compressible Computational Fluid Dynamics", +! 2021, DOI:10.5281/zenodo.4671927, +! https://github.com/ranocha/Optimized-RK-CFD +! + error_control = .true. + evolve => evolve_3sstarplus + registers = 4 + order = 4 + stages = 9 + cfl = 3.5d+00 * cfl + name_int = "Optimized 4ᵗʰ-order 9-step embedded RK4(3)9" + + allocate(c(stages), dl(stages), bt(stages), bh(stages+1), gm(3,stages)) + c (:) = [ 0.000000000000000000000000000000000000d+00, & + 2.836343531977826022543660465926414772d-01, & + 5.484073767552486705240014599676811834d-01, & + 3.687229456675706936558667052479014150d-01, & + -6.806119916032093175251948474173648331d-01, & + 3.518526451892056368706593492732753284d-01, & + 1.665941920204672094647868254892387293d+00, & + 9.715276989307335935187466054546761665d-01, & + 9.051569554420045339601721625247585643d-01 ] + dl(:) = [ 1.000000000000000000000000000000000000d+00, & + 1.262923854387806460989545005598562667d+00, & + 7.574967177560872438940839460448329992d-01, & + 5.163591158111222863455531895152351544d-01, & + -2.746333792042827389548936599648122146d-02, & + -4.382674653941770848797864513655752318d-01, & + 1.273587103668392811985704533534301656d+00, & + -6.294740045442794829622796613103492913d-01, & + 0.000000000000000000000000000000000000d+00 ] + bt(:) = [ 2.836343531977826022543660465926414772d-01, & + 9.736497978646965372894268287659773644d-01, & + 3.382358566377620380505126936670933370d-01, & + -3.584937820217850715182820651063453804d-01, & + -4.113955814725134294322006403954822487d-03, & + 1.427968962196019024010757034274849198d+00, & + 1.808467712038743032991177525728915926d-02, & + 1.605771316794521018947553625079465692d-01, & + 2.952226811394310028003810072027839487d-01 ] + bh(:) = [ 4.550655927970944948340364817140593012d-02, & + 1.175968310492638562142460384341959193d-01, & + 3.658257330515213200375475084421083608d-02, & + -5.311555834355629559010061596928357525d-03, & + 5.178250012713127329531367677410650996d-03, & + 4.954639022118682638697706200022961443d-01, & + -5.999303132737865921441409466809521699d-03, & + 9.405093434568315929035250835218733824d-02, & + 2.169318087627035072893925375820310602d-01, & + 0.000000000000000000000000000000000000d+00 ] + gm(1,:) = [ 0.000000000000000000000000000000000000d+00, & + -4.655641301259180308677051498071354582d+00, & + -7.720264924836063859141482018013692338d-01, & + -4.024423213419724605695005429153112050d+00, & + -2.129685246739018613087466942802498152d-02, & + -2.435022519234470128602335652131234586d+00, & + 1.985627480986167686791439120784668251d-02, & + -2.810790112885283952929218377438668784d-01, & + 1.689434895835535695524003319503844110d-01 ] + gm(2,:) = [ 1.000000000000000000000000000000000000d+00, & + 2.499262752607825957145627300817258023d+00, & + 5.866820365436136799319929406678132638d-01, & + 1.205141365412670762568835277881144391d+00, & + 3.474793796700868848597960521248007941d-01, & + 1.321346140128723105871355808477092220d+00, & + 3.119636324379370564023292317172847140d-01, & + 4.351419055894087609560896967082486864d-01, & + 2.359698299440788299161958168555704234d-01 ] + gm(3,:) = [ 0.000000000000000000000000000000000000d+00, & + 0.000000000000000000000000000000000000d+00, & + 0.000000000000000000000000000000000000d+00, & + 7.621037111138170045618771082985664430d-01, & + -1.981182159087218433914909510116664154d-01, & + -6.228960706317566993192689455719570179d-01, & + -3.752246993432626328289874575355102038d-01, & + -3.355436539000946543242869676125143358d-01, & + -4.560963110717484359015342341157302403d-02 ] + + betas(:) = [ 2.5d-01, -1.2d-01, 0.0d+00 ] + + call get_parameter("beta(1)", betas(1)) + call get_parameter("beta(2)", betas(2)) + call get_parameter("beta(3)", betas(3)) + + betas(:) = - betas(:) / 4.0d+00 + + case ("RK4(3)9F", "rk4(3)9f") + +! References: +! +! [1] Ranocha, H., Dalcin, L., Parsani, M., Ketcheson, D. I. +! "Coefficients of Optimized Runge-Kutta Methods with Automatic +! Step Size Control for Compressible Computational Fluid Dynamics", +! 2021, DOI:10.5281/zenodo.4671927, +! https://github.com/ranocha/Optimized-RK-CFD +! + error_control = .true. + fsal = .true. + evolve => evolve_3sstarplus + registers = 4 + order = 4 + stages = 9 + cfl = 3.5d+00 * cfl + name_int = "Optimized 4ᵗʰ-order 9-step embedded RK4(3)9 FSAL" + + allocate(c(stages), dl(stages), bt(stages), bh(stages+1), gm(3,stages)) + c (:) = [ 0.000000000000000000000000000000000000d+00, & + 2.836343005184365275160654678626695428d-01, & + 5.484076570002894365286665352032296535d-01, & + 3.687228761669438493478872632332010073d-01, & + -6.806126440140844191258463830024463902d-01, & + 3.518526124230705801739919476290327750d-01, & + 1.665941994879593315477304663913129942d+00, & + 9.715279295934715835299192116436237065d-01, & + 9.051569840159589594903399929316959062d-01 ] + dl(:) = [ 1.000000000000000000000000000000000000d+00, & + 1.262923876648114432874834923838556100d+00, & + 7.574967189685911558308119415539596711d-01, & + 5.163589453140728104667573195005629833d-01, & + -2.746327421802609557034437892013640319d-02, & + -4.382673178127944142238606608356542890d-01, & + 1.273587294602656522645691372699677063d+00, & + -6.294740283927400326554066998751383342d-01, & + 0.000000000000000000000000000000000000d+00 ] + bt(:) = [ 2.836343005184365275160654678626695428d-01, & + 9.736500104654741223716056170419660217d-01, & + 3.382359225242515288768487569778320563d-01, & + -3.584943611106183357043212309791897386d-01, & + -4.113944068471528211627210454497620358d-03, & + 1.427968894048586363415504654313371031d+00, & + 1.808470948394314017665968411915568633d-02, & + 1.605770645946802213926893453819236685d-01, & + 2.952227015964591648775833803635147962d-01 ] + bh(:) = [ 2.483675912451591196775756814283216443d-02, & + 1.866327774562103796990092260942180726d-01, & + 5.671080795936984495604436622517631183d-02, & + -3.447695439149287702616943808570747099d-03, & + 3.602245056516636472203469198006404016d-03, & + 4.545570622145088936800484247980581766d-01, & + -2.434665289427612407531544765622888855d-04, & + 6.642755361103549971517945063138312147d-02, & + 1.613697079523505006226025497715177578d-01, & + 4.955424859358438183052504342394102722d-02 ] + gm(1,:) = [ 0.000000000000000000000000000000000000d+00, & + -4.655641447335068552684422206224169103d+00, & + -7.720265099645871829248487209517314217d-01, & + -4.024436690519806086742256154738379161d+00, & + -2.129676284018530966221583708648634733d-02, & + -2.435022509790109546199372365866450709d+00, & + 1.985627297131987000579523283542615256d-02, & + -2.810791146791038566946663374735713961d-01, & + 1.689434168754859644351230590422137972d-01 ] + gm(2,:) = [ 1.000000000000000000000000000000000000d+00, & + 2.499262792574495009336242992898153462d+00, & + 5.866820377718875577451517985847920081d-01, & + 1.205146086523094569925592464380295241d+00, & + 3.474793722186732780030762737753849272d-01, & + 1.321346060965113109321230804210670518d+00, & + 3.119636464694193615946633676950358444d-01, & + 4.351419539684379261368971206040518552d-01, & + 2.359698130028753572503744518147537768d-01 ] + gm(3,:) = [ 0.000000000000000000000000000000000000d+00, & + 0.000000000000000000000000000000000000d+00, & + 0.000000000000000000000000000000000000d+00, & + 7.621006678721315291614677352949377871d-01, & + -1.981182504339400567765766904309673119d-01, & + -6.228959218699007450469629366684127462d-01, & + -3.752248380775956442989480369774937099d-01, & + -3.355438309135169811915662336248989661d-01, & + -4.560955005031121479972862973705108039d-02 ] + + betas(:) = [ 3.8d-01, -1.8d-01, 1.0d-02 ] + + call get_parameter("beta(1)", betas(1)) + call get_parameter("beta(2)", betas(2)) + call get_parameter("beta(3)", betas(3)) + + betas(:) = - betas(:) / 4.0d+00 + + case ("RK5(4)10", "rk5(4)10") + +! References: +! +! [1] Ranocha, H., Dalcin, L., Parsani, M., Ketcheson, D. I. +! "Coefficients of Optimized Runge-Kutta Methods with Automatic +! Step Size Control for Compressible Computational Fluid Dynamics", +! 2021, DOI:10.5281/zenodo.4671927, +! https://github.com/ranocha/Optimized-RK-CFD +! + error_control = .true. + evolve => evolve_3sstarplus + registers = 4 + order = 5 + stages = 10 + cfl = 4.2d+00 * cfl + name_int = "Optimized 5ᵗʰ-order 10-step embedded RK5(4)10" + + allocate(c(stages), dl(stages), bt(stages), bh(stages+1), gm(3,stages)) + c (:) = [ 0.000000000000000000000000000000000000d+00, & + 2.597883575710995826783320802193635406d-01, & + 9.904573115730917688557891428202061598d-02, & + 2.155511882303785204133426661931565216d-01, & + 5.007950078421880417512789524851012021d-01, & + 5.592251914858131230054392022144328176d-01, & + 5.449986973408778242805929551952000165d-01, & + 7.615224662599497796472095353126697300d-01, & + 8.427062083059167761623893618875787414d-01, & + 9.152209807185253394871325258038753352d-01 ] + dl(:) = [ 1.000000000000000000000000000000000000d+00, & + -1.331778409133849616712007380176762548d-01, & + 8.260422785246030254485064732649153253d-01, & + 1.513700430513332405798616943654007796d+00, & + -1.305810063177048110528482211982726539d+00, & + 3.036678789342507704281817524408221954d+00, & + -1.449458267074592489788800461540171106d+00, & + 3.834313873320957483471400258279635203d+00, & + 4.122293971923324492772059928094971199d+00, & + 0.000000000000000000000000000000000000d+00 ] + bt(:) = [ 2.597883575710995826783320802193635406d-01, & + 1.777008800169541694837687556103565007d-02, & + 2.481636637328140606807905234325691851d-01, & + 7.941736827560429420202759490815682546d-01, & + 3.885391296871822541486945325814526190d-01, & + 1.455051664264339366757555740296587660d-01, & + 1.587517379462528932413419955691782412d-01, & + 1.650605631567659573994022720500446501d-01, & + 2.118093299943235065178000892467421832d-01, & + 1.559392340339606299335442956580114440d-01 ] + bh(:) = [ 5.734588484676193812418453938089759359d-02, & + 1.971447518039733870541652912891291496d-02, & + 7.215296605683716720707226840456658773d-02, & + 1.739659489807939956977075317768151880d-01, & + 3.703693600445487815015171515640585668d-01, & + -1.215599039055065009827765147821222534d-01, & + 1.180372945491121604465067725859678821d-01, & + 4.155688823364870056536983972605056553d-02, & + 1.227886627910379901351569893551486490d-01, & + 1.456284232223684285998448928597043056d-01, & + 0.000000000000000000000000000000000000d+00 ] + gm(1,:) = [ 0.000000000000000000000000000000000000d+00, & + 4.043660078504695837542588769963326988d-01, & + -8.503427464263185087039788184485627962d-01, & + -6.950894167072419998080989313353063399d+00, & + 9.238765225328278557805080247596562995d-01, & + -2.563178039957404359875124580586147888d+00, & + 2.545744869966347362604059848503340890d-01, & + 3.125831733863168874151935287174374515d-01, & + -7.007114800567584871263283872289072079d-01, & + 4.839620970980726631935174740648996010d-01 ] + gm(2,:) = [ 1.000000000000000000000000000000000000d+00, & + 6.871467069752345566001768382316915820d-01, & + 1.093024760468898686510433898645775908d+00, & + 3.225975382330161123625348062949430509d+00, & + 1.041153700841396427100436517666787823d+00, & + 1.292821488864702752767390075072674807d+00, & + 7.391462769297006312785029455392854586d-01, & + 1.239129257039300081860496157739352186d-01, & + 1.842753479366766790220633908793933781d-01, & + 5.712788942697077644959290025755003720d-02 ] + gm(3,:) = [ 0.000000000000000000000000000000000000d+00, & + 0.000000000000000000000000000000000000d+00, & + 0.000000000000000000000000000000000000d+00, & + -2.393405159342139386425044844626597490d+00, & + -1.902854422095986544338294743445530533d+00, & + -2.820042210583207174321941694153843259d+00, & + -1.832698464130564949123807896975136336d+00, & + -2.199094510750697865007677774395365522d-01, & + -4.082430660384876496971887725512427800d-01, & + -1.377669791121207993339861855818881150d-01 ] + + betas(:) = [ 4.7d-01, -2.0d-01, 6.0d-02 ] + + call get_parameter("beta(1)", betas(1)) + call get_parameter("beta(2)", betas(2)) + call get_parameter("beta(3)", betas(3)) + + betas(:) = - betas(:) / 5.0d+00 + + case ("RK5(4)10F", "rk5(4)10f") + +! References: +! +! [1] Ranocha, H., Dalcin, L., Parsani, M., Ketcheson, D. I. +! "Coefficients of Optimized Runge-Kutta Methods with Automatic +! Step Size Control for Compressible Computational Fluid Dynamics", +! 2021, DOI:10.5281/zenodo.4671927, +! https://github.com/ranocha/Optimized-RK-CFD +! + error_control = .true. + fsal = .true. + evolve => evolve_3sstarplus + registers = 4 + order = 5 + stages = 10 + cfl = 4.2d+00 * cfl + name_int = "Optimized 5ᵗʰ-order 10-step embedded RK5(4)10 FSAL" + + allocate(c(stages), dl(stages), bt(stages), bh(stages+1), gm(3,stages)) + c (:) = [ 0.000000000000000000000000000000000000d+00, & + 2.597883554788674084039539165398464630d-01, & + 9.904573247592460887087003212056568980d-02, & + 2.155511890524058691860390281856497503d-01, & + 5.007950088969676776844289399972611534d-01, & + 5.592251911688643533787800688765883636d-01, & + 5.449986978853637084972622392134732553d-01, & + 7.615224694532590139829150720490417596d-01, & + 8.427062083267360939805493320684741215d-01, & + 9.152209805057669959657927210873423883d-01 ] + dl(:) = [ 1.000000000000000000000000000000000000d+00, & + -1.331778419508803397033287009506932673d-01, & + 8.260422814750207498262063505871077303d-01, & + 1.513700425755728332485300719652378197d+00, & + -1.305810059935023735972298885749903694d+00, & + 3.036678802924163246003321318996156380d+00, & + -1.449458274398895177922690618003584514d+00, & + 3.834313899176362315089976408899373409d+00, & + 4.122293760012985409330881631526514714d+00, & + 0.000000000000000000000000000000000000d+00 ] + bt(:) = [ 2.597883554788674084039539165398464630d-01, & + 1.777008889438867858759149597539211023d-02, & + 2.481636629715501931294746189266601496d-01, & + 7.941736871152005775821844297293296135d-01, & + 3.885391285642019129575902994397298066d-01, & + 1.455051657916305055730603387469193768d-01, & + 1.587517385964749337690916959584348979d-01, & + 1.650605617880053419242434594242509601d-01, & + 2.118093284937153836908655490906875007d-01, & + 1.559392342362059886106995325687547506d-01 ] + bh(:) = [ -2.019255440012066080909442770590267512d-02, & + 2.737903480959184339932730854141598275d-02, & + 3.028818636145965534365173822296811090d-01, & + -3.656843880622222190071445247906780540d-02, & + 3.982664774676767729863101188528827405d-01, & + -5.715959421140685436681459970502471634d-02, & + 9.849855103848558320961101178888983150d-02, & + 6.654601552456084978615342374581437947d-02, & + 9.073479542748112726465375642050504556d-02, & + 8.432289325330803924891866923939606351d-02, & + 4.529095628204896774513180907141004447d-02 ] + gm(1,:) = [ 0.000000000000000000000000000000000000d+00, & + 4.043660121685749695640462197806189975d-01, & + -8.503427289575839690883191973980814832d-01, & + -6.950894175262117526410215315179482885d+00, & + 9.238765192731084931855438934978371889d-01, & + -2.563178056509891340215942413817786020d+00, & + 2.545744879365226143946122067064118430d-01, & + 3.125831707411998258746812355492206137d-01, & + -7.007114414440507927791249989236719346d-01, & + 4.839621016023833375810172323297465039d-01 ] + gm(2,:) = [ 1.000000000000000000000000000000000000d+00, & + 6.871467028161416909922221357014564412d-01, & + 1.093024748914750833700799552463885117d+00, & + 3.225975379607193001678365742708874597d+00, & + 1.041153702510101386914019859778740444d+00, & + 1.292821487912164945157744726076279306d+00, & + 7.391462755788122847651304143259254381d-01, & + 1.239129251371800313941948224441873274d-01, & + 1.842753472370123193132193302369345580d-01, & + 5.712788998796583446479387686662738843d-02 ] + gm(3,:) = [ 0.000000000000000000000000000000000000d+00, & + 0.000000000000000000000000000000000000d+00, & + 0.000000000000000000000000000000000000d+00, & + -2.393405133244194727221124311276648940d+00, & + -1.902854422421760920850597670305403139d+00, & + -2.820042207399977261483046412236557428d+00, & + -1.832698465277380999601896111079977378d+00, & + -2.199094483084671192328083958346519535d-01, & + -4.082430635847870963724591602173546218d-01, & + -1.377669797880289713535665985132703979d-01 ] + + betas(:) = [ 4.5d-01, -1.3d-01, 0.0d+00 ] + + call get_parameter("beta(1)", betas(1)) + call get_parameter("beta(2)", betas(2)) + call get_parameter("beta(3)", betas(3)) + + betas(:) = - betas(:) / 5.0d+00 + case default if (verbose) then @@ -327,7 +856,8 @@ module evolution " 'ssprk(m,2)', 'ssprk(4,3)', 'ssprk(5,3)'," // & " 'ssprk(m,3)', 'ssprk(10,4)'." write(*,"(1x,a)") "Available embedded methods: 'ssprk3(2)4'," // & - " 'ssprk3(2)m'." + " 'ssprk3(2)m', 'rk3(2)5', 'rk4(3)9', 'rk5(4)10',"// & + "'rk3(2)5f', 'rk4(3)9f', 'rk5(4)10f'." end if status = 1 @@ -411,6 +941,11 @@ module evolution if (allocated(adecay)) deallocate(adecay) if (allocated(aglm)) deallocate(aglm) + if (allocated(c)) deallocate(c) + if (allocated(dl)) deallocate(dl) + if (allocated(bt)) deallocate(bt) + if (allocated(bh)) deallocate(bh) + if (allocated(gm)) deallocate(gm) ! reset the status flag ! @@ -534,10 +1069,6 @@ module evolution call start_timer(ima) #endif /* PROFILE */ -! find new time step -! - call new_time_step() - ! advance the solution using the selected method ! call evolve() @@ -591,105 +1122,78 @@ module evolution ! !=============================================================================== ! -! subroutine NEW_TIME_STEP: -! ------------------------ +! subroutine INITIALIZE_TIME_STEP: +! ------------------------------- ! -! Subroutine estimates the new time step from the maximum speed in the system -! and source term constraints. +! Subroutine estimates the initial time step. ! -! Arguments: +! References: ! +! [1] Courant, R., Friedrichs, K., Lewy, H., +! "Über die partiellen Differenzengleichungen der mathematischen Physik", +! Mathematische Annalen, 1928, volume 100, pages 32–74, +! DOI: 10.1007/BF01448839 +! [2] Hairer, E., Nørsett, S. P., Wanner, G., +! "Solving Ordinary Differential Equations I: Nonstiff Problems", +! Springer-Verlag Berlin Heidelberg, 1993, +! ISBN: 978-3-540-56670-0, DOI: 10.1007/978-3-540-78862-1 ! !=============================================================================== ! - subroutine new_time_step() + subroutine initialize_time_step() -! include external procedures -! - use equations , only : maxspeed, cmax, cmax2 -#ifdef MPI - use mpitools , only : reduce_maximum -#endif /* MPI */ - -! include external variables -! - use blocks , only : block_data, list_data - use coordinates , only : adx, ady + use blocks , only : block_data, list_data, get_nleafs + use coordinates, only : nc => ncells, nb, ne, adx, ady #if NDIMS == 3 - use coordinates , only : adz + use coordinates, only : adz #endif /* NDIMS == 3 */ - use sources , only : viscosity, resistivity + use equations , only : nf, maxspeed, cmax, cmax2 +#ifdef MPI + use mpitools , only : reduce_maximum, reduce_sum +#endif /* MPI */ + use sources , only : viscosity, resistivity, update_sources -! local variables are not implicit by default -! implicit none -! local pointers -! type(block_data), pointer :: pdata -! local variables -! - integer :: mlev - real(kind=8) :: cm, dx_min + integer :: mlev + real(kind=8) :: dx_min, fnorm, h0, h1 + real(kind=8), dimension(3) :: d -! local parameters -! - real(kind=8), parameter :: eps = tiny(cmax) +#if NDIMS == 3 + real(kind=8), dimension(nf,nc,nc,nc) :: sc, df +#else /* NDIMS == 3 */ + real(kind=8), dimension(nf,nc,nc, 1) :: sc, df +#endif /* NDIMS == 3 */ + + real(kind=8), parameter :: eps = tiny(cmax) ! !------------------------------------------------------------------------------- ! #ifdef PROFILE -! start accounting time for new time step estimation -! - call start_timer(imt) + call start_timer(imi) #endif /* PROFILE */ -! reset the maximum speed, and the highest level -! cmax = eps mlev = 1 -! assign pdata with the first block on the data block list -! pdata => list_data - -! iterate over all data blocks in order to find the maximum speed among them -! and the highest level which is required to obtain the minimum spacial step -! do while (associated(pdata)) -! update the maximum level -! mlev = max(mlev, pdata%meta%level) + cmax = max(cmax, maxspeed(pdata%q(:,:,:,:))) -! get the maximum speed for the current block -! - cm = maxspeed(pdata%q(:,:,:,:)) - -! update the maximum speed -! - cmax = max(cmax, cm) - -! assign pdata to the next block -! pdata => pdata%next - - end do ! over data blocks + end do #ifdef MPI -! reduce maximum speed and level over all processors -! call reduce_maximum(cmax) call reduce_maximum(mlev) #endif /* MPI */ -! calculate the squared cmax -! cmax2 = cmax * cmax -! find the smallest spacial step -! #if NDIMS == 2 dx_min = min(adx(mlev), ady(mlev)) #endif /* NDIMS == 2 */ @@ -697,26 +1201,197 @@ module evolution dx_min = min(adx(mlev), ady(mlev), adz(mlev)) #endif /* NDIMS == 3 */ -! calculate the new time step +! determine the time step due to the stability ! - dth = cfl * dx_min / max(cmax & + dth = dx_min / max(cmax & + 2.0d+00 * max(viscosity, resistivity) / dx_min, eps) + dt = cfl * dth -! round the time +! determine the time step due to the error ! - dtn = dth if (error_control) then - if (dte > 0.0d+00) then - dtn = min(dtn, dte) + + fnorm = 1.0d+00 / (get_nleafs() * nc**NDIMS) + + d(:) = 0.0d+00 + + pdata => list_data + do while (associated(pdata)) + + call update_increment(pdata) + call update_sources(pdata, 0.0d+00, 0.0d+00, pdata%du(:,:,:,:)) + +#if NDIMS == 3 + sc(:,:,:,:) = atol + rtol * abs(pdata%uu(:,nb:ne,nb:ne,nb:ne,1)) + + d(1) = d(1) + sum((pdata%uu(:,nb:ne,nb:ne,nb:ne,1) / sc(:,:,:,:))**2) + d(2) = d(2) + sum((pdata%du(:,nb:ne,nb:ne,nb:ne) / sc(:,:,:,:))**2) +#else /* NDIMS == 3 */ + sc(:,:,:,:) = atol + rtol * abs(pdata%uu(:,nb:ne,nb:ne, : ,1)) + + d(1) = d(1) + sum((pdata%uu(:,nb:ne,nb:ne, : ,1) / sc(:,:,:,:))**2) + d(2) = d(2) + sum((pdata%du(:,nb:ne,nb:ne, : ) / sc(:,:,:,:))**2) +#endif /* NDIMS == 3 */ + + pdata => pdata%next + end do + +#ifdef MPI + call reduce_sum(d(1:2)) +#endif /* MPI */ + + d(1:2) = sqrt(fnorm * d(1:2)) + + if (minval(d(1:2)) >= 1.0d-05) then + h0 = 1.0d-02 * d(1) / d(2) else - dte = dth + h0 = 1.0d-06 end if + + pdata => list_data + do while (associated(pdata)) + + pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) + h0 * pdata%du(:,:,:,:) + + pdata%u => pdata%uu(:,:,:,:,2) + + pdata => pdata%next + end do + + call update_variables(time + h0, h0) + + pdata => list_data + do while (associated(pdata)) + +#if NDIMS == 3 + df(:,:,:,:) = pdata%du(:,nb:ne,nb:ne,nb:ne) + sc(:,:,:,:) = atol + rtol * abs(pdata%uu(:,nb:ne,nb:ne,nb:ne,1)) +#else /* NDIMS == 3 */ + df(:,:,:,:) = pdata%du(:,nb:ne,nb:ne, : ) + sc(:,:,:,:) = atol + rtol * abs(pdata%uu(:,nb:ne,nb:ne, : ,1)) +#endif /* NDIMS == 3 */ + + call update_increment(pdata) + call update_sources(pdata, time + h0, 0.0d+00, pdata%du(:,:,:,:)) + +#if NDIMS == 3 + df(:,:,:,:) = df(:,:,:,:) - pdata%du(:,nb:ne,nb:ne,nb:ne) +#else /* NDIMS == 3 */ + df(:,:,:,:) = df(:,:,:,:) - pdata%du(:,nb:ne,nb:ne, : ) +#endif /* NDIMS == 3 */ + d(3) = d(3) + sum((df(:,:,:,:) / sc(:,:,:,:))**2) + + pdata%u => pdata%uu(:,:,:,:,1) + + pdata => pdata%next + end do + +#ifdef MPI + call reduce_sum(d(3:3)) +#endif /* MPI */ + + d(3) = sqrt(fnorm * d(3)) / h0 + + if (max(d(2), d(3)) > 1.0d-15) then + h1 = (1.0d-02 / max(d(2), d(3)))**(1.0d+00 / (order + 1)) + else + h1 = max(1.0d-06, 1.0d-03 * h0) + end if + + dte = min(1.0d+02 * h0, h1) + dt = min(dt, dte) + end if - dt = min(dtn, dtp) #ifdef PROFILE -! stop accounting time for new time step estimation + call stop_timer(imi) +#endif /* PROFILE */ + +!------------------------------------------------------------------------------- ! + end subroutine initialize_time_step +! +!=============================================================================== +! +! subroutine NEW_TIME_STEP: +! ------------------------ +! +! Subroutine estimates the new time step from the maximum speed in the system +! and source term constraints. +! +! References: +! +! [1] Courant, R., Friedrichs, K., Lewy, H., +! "Über die partiellen Differenzengleichungen der mathematischen Physik", +! Mathematische Annalen, 1928, volume 100, pages 32–74, +! DOI: 10.1007/BF01448839 +! +!=============================================================================== +! + subroutine new_time_step() + + use blocks , only : block_data, list_data + use coordinates , only : adx, ady +#if NDIMS == 3 + use coordinates , only : adz +#endif /* NDIMS == 3 */ + use equations , only : maxspeed, cmax, cmax2 +#ifdef MPI + use mpitools , only : reduce_maximum +#endif /* MPI */ + use sources , only : viscosity, resistivity + + implicit none + + type(block_data), pointer :: pdata + + integer :: mlev + real(kind=8) :: cm, dx_min + + real(kind=8), parameter :: eps = tiny(cmax) +! +!------------------------------------------------------------------------------- +! +#ifdef PROFILE + call start_timer(imt) +#endif /* PROFILE */ + + cmax = eps + mlev = 1 + + pdata => list_data + + do while (associated(pdata)) + + mlev = max(mlev, pdata%meta%level) + cm = maxspeed(pdata%q(:,:,:,:)) + cmax = max(cmax, cm) + + pdata => pdata%next + + end do + +#ifdef MPI + call reduce_maximum(cmax) + call reduce_maximum(mlev) +#endif /* MPI */ + + cmax2 = cmax * cmax + +#if NDIMS == 2 + dx_min = min(adx(mlev), ady(mlev)) +#endif /* NDIMS == 2 */ +#if NDIMS == 3 + dx_min = min(adx(mlev), ady(mlev), adz(mlev)) +#endif /* NDIMS == 3 */ + + dth = dx_min / max(cmax & + + 2.0d+00 * max(viscosity, resistivity) / dx_min, eps) + dt = cfl * dth + dt = min(dt, dtp) + if (error_control) dt = min(dt, dte) + +#ifdef PROFILE call stop_timer(imt) #endif /* PROFILE */ @@ -2448,6 +3123,256 @@ module evolution ! !=============================================================================== ! +! subroutine EVOLVE_3SSTARPLUS: +! ---------------------------- +! +! Subroutine advances the solution by one time step using the minimum +! storage implementation of 3S*+ methods. +! +! References: +! +! [1] Ranocha, H., Dalcin, L., Parsani, M., Ketcheson, D. I. +! "Optimized Runge-Kutta Methods with Automatic Step Size Control +! for Compressible Computational Fluid Dynamics", +! 2021, arXiv:2104.06836 +! [2] Gottlieb, S., Ketcheson, D., Shu C. - W., +! "Strong stability preserving Runge-Kutta and multistep +! time discretizations", +! World Scientific Publishing, 2011 +! +!=============================================================================== +! + subroutine evolve_3sstarplus() + + use blocks , only : block_data, list_data + use boundaries , only : boundary_fluxes + use equations , only : errors, ibp, cmax + use sources , only : update_sources + + implicit none + + type(block_data), pointer :: pdata + + logical :: test + integer :: i, l, nrej + real(kind=8) :: tm, dtm, fc +! +!------------------------------------------------------------------------------- +! +#ifdef PROFILE + call start_timer(imu) +#endif /* PROFILE */ + + test = .true. + nrej = 0 + +! at the entry point we assume the previous solution of conserved variables U(n) +! is stored in pdata%u(:,:,:,:,1) and the primitive variables are already +! updated from this array; +! +! pdata%uu(:,:,:,:,1) - the previous solution, U(n) +! pdata%uu(:,:,:,:,2) - the new 3rd order solution, U(1) +! pdata%uu(:,:,:,:,3) - the intermediate solution, U(2) +! pdata%uu(:,:,:,:,4) - the new 2nd order solution, U(3) +! + do while(test) + +!= preparatory step: U(1) = U(n), U(2) = 0, U(3) = U(n) +! + pdata => list_data + do while (associated(pdata)) + + pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) + pdata%uu(:,:,:,:,3) = pdata%uu(:,:,:,:,1) + pdata%uu(:,:,:,:,4) = pdata%uu(:,:,:,:,1) + + pdata%u => pdata%uu(:,:,:,:,2) + + pdata => pdata%next + end do + +!= 1st stage +! + if (fsal_update) then + pdata => list_data + do while (associated(pdata)) + + call update_increment(pdata) + call update_sources(pdata, time, 0.0d+00, pdata%du(:,:,:,:)) + + pdata => pdata%next + end do + + call boundary_fluxes() + end if + + pdata => list_data + do while (associated(pdata)) + + pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) & + + bt(1) * dt * pdata%du(:,:,:,:) + pdata%uu(:,:,:,:,4) = pdata%uu(:,:,:,:,4) & + + bh(1) * dt * pdata%du(:,:,:,:) + + pdata => pdata%next + end do + + call update_variables(tm, dtm) + +!= remaining stages +! + do i = 2, stages + + tm = time + c(i) * dt + dtm = (c(i) - c(i-1)) * dt + + pdata => list_data + do while (associated(pdata)) + + call update_increment(pdata) + call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) + + pdata => pdata%next + end do + + call boundary_fluxes() + + pdata => list_data + do while (associated(pdata)) + + pdata%uu(:,:,:,:,3) = pdata%uu(:,:,:,:,3) & + + dl(i) * pdata%uu(:,:,:,:,2) + pdata%uu(:,:,:,:,2) = gm(1,i) * pdata%uu(:,:,:,:,2) & + + gm(2,i) * pdata%uu(:,:,:,:,3) & + + gm(3,i) * pdata%uu(:,:,:,:,1) & + + bt(i) * dt * pdata%du(:,:,:,:) + pdata%uu(:,:,:,:,4) = pdata%uu(:,:,:,:,4) & + + bh(i) * dt * pdata%du(:,:,:,:) + + pdata => pdata%next + end do + + call update_variables(tm, dtm) + + end do ! i = 2, stages + +!= extra step: decay the magnetic divergence potential of both solutions +! + if (ibp > 0) then + adecay(:) = exp(aglm(:) * cmax * dt) + + pdata => list_data + do while (associated(pdata)) + + pdata%uu(ibp,:,:,:,2) = adecay(pdata%meta%level) & + * pdata%uu(ibp,:,:,:,2) + pdata%uu(ibp,:,:,:,4) = adecay(pdata%meta%level) & + * pdata%uu(ibp,:,:,:,4) + + pdata => pdata%next + end do + end if + +!= extra step: FSAL +! + if (fsal) then + + tm = time + dt + + pdata => list_data + do while (associated(pdata)) + + call update_increment(pdata) + call update_sources(pdata, tm, 0.0d+00, pdata%du(:,:,:,:)) + + pdata => pdata%next + end do + + call boundary_fluxes() + + pdata => list_data + do while (associated(pdata)) + + pdata%uu(:,:,:,:,4) = pdata%uu(:,:,:,:,4) & + + bh(stages+1) * dt * pdata%du(:,:,:,:) + + pdata => pdata%next + end do + + end if + +!= error estimation +! + call update_errors(2, 4) + +! calculate the tolerance and estimate the next time step due to the error +! + errtol = maxval(errors) + errs(1) = errtol + fc = product(errs(:)**betas(:)) + fc = 1.0d+00 + chi * atan((fc - 1.0d+00) / chi) + dte = dt * fc + + if (fc > fac .or. nrej >= mrej) then + test = .false. + + errs(3) = errs(2) + errs(2) = errs(1) + + if (fsal) fsal_update = .true. + else + dt = dte + + nrej = nrej + 1 ! rejection count in the current step + nrejections = nrejections + 1 + +! since the solution was rejected, we have to revert the primitive variables +! to the previous state +! + pdata => list_data + do while (associated(pdata)) + + pdata%u => pdata%uu(:,:,:,:,1) + + pdata => pdata%next + end do + + call update_variables(tm, dtm) + + if (fsal) fsal_update = .false. + end if + + niterations = niterations + 1 + + end do + +!= final step: U(n+1) = U(1) +! + dtm = dt + tm = time + dt + + pdata => list_data + do while (associated(pdata)) + + pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2) + + pdata%u => pdata%uu(:,:,:,:,1) + + pdata => pdata%next + end do + + call update_variables(tm, dtm) + +#ifdef PROFILE + call stop_timer(imu) +#endif /* PROFILE */ + +!------------------------------------------------------------------------------- +! + end subroutine evolve_3sstarplus +! +!=============================================================================== +! ! subroutine UPDATE_INCREMENT: ! --------------------------- ! diff --git a/sources/integrals.F90 b/sources/integrals.F90 index 64b7b9c..57b5f42 100644 --- a/sources/integrals.F90 +++ b/sources/integrals.F90 @@ -438,7 +438,7 @@ module integrals use equations , only : ien, imx, imy, imz use equations , only : magnetized, adiabatic_index, csnd use equations , only : errors - use evolution , only : step, time, dtn, dth, dte + use evolution , only : step, time, dt, dth, dte use forcing , only : einj, rinj, arms use helpers , only : flush_and_sync #ifdef MPI @@ -1268,7 +1268,7 @@ module integrals ! write down the integrals and statistics to appropriate files ! if (stored) then - write(funit,"(i9,13(1x,1es18.8e3))") step, time, dtn, inarr(1:10), arms + write(funit,"(i9,13(1x,1es18.8e3))") step, time, dt, inarr(1:10), arms write(sunit,"(i9,23(1x,1es18.8e3))") step, time & , avarr(1), mnarr(1), mxarr(1) & , avarr(2), mnarr(2), mxarr(2) & diff --git a/sources/io.F90 b/sources/io.F90 index 11a6b4f..c04b64b 100644 --- a/sources/io.F90 +++ b/sources/io.F90 @@ -213,7 +213,7 @@ module io public :: restart_snapshot_number, restart_from_snapshot public :: read_snapshot_parameter public :: read_restart_snapshot, write_restart_snapshot, write_snapshot - public :: next_tout, precise_snapshots + public :: update_dtp !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! @@ -916,31 +916,32 @@ module io ! !=============================================================================== ! -! function NEXT_TOUT: -! ------------------ -! -! Function returns the next data snapshot time. +! subroutine UPDATE_DTP: +! --------------------- ! +! Subroutine updates dtp from module EVOLUTION. ! !=============================================================================== ! - real(kind=8) function next_tout() + subroutine update_dtp() + + use evolution, only : time, dtp -! local variables are not implicit by default -! implicit none ! !------------------------------------------------------------------------------- ! - if (hsnap > 0.0d+00) then - next_tout = tsnap - else - next_tout = huge(hsnap) + if (precise_snapshots) then + if (tsnap > time) then + dtp = tsnap - time + else + dtp = hsnap + endif end if !------------------------------------------------------------------------------- ! - end function next_tout + end subroutine update_dtp ! !=============================================================================== !! @@ -1260,7 +1261,7 @@ module io use blocks , only : change_blocks_process use coordinates , only : nn => bcells, ncells, nghosts use coordinates , only : xmin, xmax, ymin, ymax, zmin, zmax - use evolution , only : step, time, dt, dtn, dte + use evolution , only : step, time, dt, dth, dte use evolution , only : niterations, nrejections, errs use forcing , only : nmodes, fcoefs, einj use hash , only : xxh64 @@ -1415,8 +1416,8 @@ module io read(svalue, fmt = *) time case('dt') read(svalue, fmt = *) dt - case('dtn') - read(svalue, fmt = *) dtn + case('dth') + read(svalue, fmt = *) dth case('dte') read(svalue, fmt = *) dte case('niterations') @@ -2248,7 +2249,7 @@ module io #endif /* NDIMS == 3 */ use coordinates , only : bdims => domain_base_dims use equations , only : eqsys, eos, nv - use evolution , only : step, time, dt, dtn, dte, cfl, glm_alpha, errs + use evolution , only : step, time, dt, dth, dte, cfl, glm_alpha, errs use evolution , only : atol, rtol, mrej, niterations, nrejections use forcing , only : nmodes, fcoefs, einj use iso_fortran_env, only : error_unit @@ -2396,7 +2397,7 @@ module io call write_attribute_xml(lun, "step" , step) call write_attribute_xml(lun, "time" , time) call write_attribute_xml(lun, "dt" , dt) - call write_attribute_xml(lun, "dtn" , dtn) + call write_attribute_xml(lun, "dth" , dth) call write_attribute_xml(lun, "dte" , dte) call write_attribute_xml(lun, "cfl" , cfl) call write_attribute_xml(lun, "glm_alpha", glm_alpha) @@ -4184,7 +4185,7 @@ module io use coordinates , only : xmin, xmax, ymin, ymax, zmin, zmax use coordinates , only : periodic use equations , only : eqsys, eos, adiabatic_index, csnd - use evolution , only : step, time, dt, dtn, dte, cfl, glm_alpha, errs + use evolution , only : step, time, dt, dth, dte, cfl, glm_alpha, errs use evolution , only : atol, rtol, mrej, niterations, nrejections use forcing , only : nmodes, einj, fcoefs use hdf5 , only : hid_t @@ -4287,7 +4288,7 @@ module io call write_attribute(gid, 'zmax', zmax) call write_attribute(gid, 'time', time) call write_attribute(gid, 'dt' , dt ) - call write_attribute(gid, 'dtn' , dtn ) + call write_attribute(gid, 'dth' , dth ) call write_attribute(gid, 'dte' , dte ) call write_attribute(gid, 'cfl' , cfl ) call write_attribute(gid, 'glm_alpha', glm_alpha) @@ -4385,7 +4386,7 @@ module io use blocks , only : get_mblocks, get_dblocks, get_nleafs use coordinates , only : ncells use coordinates , only : xmin, xmax, ymin, ymax, zmin, zmax - use evolution , only : step, time, dt, dtn, dte + use evolution , only : step, time, dt, dth, dte use evolution , only : niterations, nrejections, errs use forcing , only : nmodes, fcoefs use hdf5 , only : hid_t @@ -4460,7 +4461,7 @@ module io call read_attribute(gid, 'zmax', zmax) call read_attribute(gid, 'time', time) call read_attribute(gid, 'dt' , dt ) - call read_attribute(gid, 'dtn' , dtn ) + call read_attribute(gid, 'dth' , dth ) call read_attribute(gid, 'dte' , dte ) call read_attribute(gid, 'errs(1)', errs(1)) call read_attribute(gid, 'errs(2)', errs(2))