!!****************************************************************************** !! !! This file is part of the AMUN source code, a program to perform !! Newtonian or relativistic magnetohydrodynamical simulations on uniform or !! adaptive mesh. !! !! Copyright (C) 2008-2021 Grzegorz Kowal !! !! This program is free software: you can redistribute it and/or modify !! it under the terms of the GNU General Public License as published by !! the Free Software Foundation, either version 3 of the License, or !! (at your option) any later version. !! !! This program is distributed in the hope that it will be useful, !! but WITHOUT ANY WARRANTY; without even the implied warranty of !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !! GNU General Public License for more details. !! !! You should have received a copy of the GNU General Public License !! along with this program. If not, see . !! !!****************************************************************************** !! !! module: EVOLUTION !! !! This module provides an interface for temporal integration with !! the stability handling. New integration methods can be added by !! implementing more evolve_* subroutines. !! !!****************************************************************************** ! module evolution #ifdef PROFILE ! import external subroutines ! use timers, only : set_timer, start_timer, stop_timer #endif /* PROFILE */ ! module variables are not implicit by default ! implicit none #ifdef PROFILE ! timer indices ! integer , save :: imi, ima, imt, imu, imf, iui, imv #endif /* PROFILE */ ! interfaces for procedure pointers ! abstract interface subroutine evolve_iface() end subroutine subroutine update_errors_iface(n, m) integer, intent(in) :: n, m end subroutine end interface ! pointer to the temporal integration subroutine ! procedure(evolve_iface), pointer, save :: evolve => null() ! pointer to the error update subroutine ! procedure(update_errors_iface), pointer, save :: update_errors => null() ! 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 ! coefficients controlling the decay of scalar potential ѱ ! real(kind=8) , save :: glm_alpha = 5.0d-01 ! time variables ! integer , save :: step = 0 real(kind=8) , save :: time = 0.0d+00 real(kind=8) , save :: dt = 1.0d+00 real(kind=8) , save :: dth = 1.0d+00 real(kind=8) , save :: dte = 0.0d+00 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, ! the total number of integration iterations, the number of rejected iterations ! real(kind=8) , save :: atol = 1.0d-04 real(kind=8) , save :: rtol = 1.0d-04 real(kind=8) , save :: fac = 9.0d-01 real(kind=8) , save :: facmin = 1.0d-01 real(kind=8) , save :: facmax = 5.0d+00 real(kind=8) , save :: errtol = 1.0d+00 real(kind=8) , save :: chi = 1.0d+00 integer , save :: mrej = 5 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 ! GLM variables ! real(kind=8), dimension(:), allocatable, save :: adecay, aglm ! by default everything is private ! private ! declare public subroutines ! public :: initialize_evolution, finalize_evolution, print_evolution public :: initialize_time_step, new_time_step, advance ! declare public variables ! public :: step, time, dt, dth, dte, dtp, cfl, glm_alpha, registers public :: atol, rtol, mrej, niterations, nrejections, errs, errtol public :: error_control !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! contains ! !=============================================================================== !! !!*** PUBLIC SUBROUTINES ***************************************************** !! !=============================================================================== ! !=============================================================================== ! ! subroutine INITIALIZE_EVOLUTION: ! ------------------------------- ! ! Subroutine initializes module EVOLUTION by setting its parameters. ! ! Arguments: ! ! verbose - a logical flag turning the information printing; ! status - an integer flag for error return value; ! !=============================================================================== ! subroutine initialize_evolution(verbose, status) ! include external procedures ! use coordinates, only : toplev, adx, ady, adz use parameters , only : get_parameter ! local variables are not implicit by default ! implicit none ! subroutine arguments ! logical, intent(in) :: verbose integer, intent(out) :: status ! local variables ! character(len=255) :: integration = "rk2", error_norm = "l2" integer :: n ! local parameters ! character(len=*), parameter :: loc = 'EVOLUTION::initialize_evolution()' ! !------------------------------------------------------------------------------- ! #ifdef PROFILE ! set timer descriptions ! call set_timer('evolution:: initialization' , imi) call set_timer('evolution:: solution advance', ima) call set_timer('evolution:: new time step' , imt) call set_timer('evolution:: solution update' , imu) call set_timer('evolution:: flux update' , imf) call set_timer('evolution:: increment update', iui) call set_timer('evolution:: variable update' , imv) ! start accounting time for module initialization/finalization ! call start_timer(imi) #endif /* PROFILE */ ! reset the status flag ! status = 0 ! get the integration method and the value of the CFL coefficient ! call get_parameter("time_advance", integration) call get_parameter("stages" , stages ) call get_parameter("cfl" , cfl ) call get_parameter("glm_alpha" , glm_alpha ) ! get the absolute and relative tolerances, the maximum number of passes for ! the adaptive step, and limiting factors ! call get_parameter("absolute_tolerance", atol) call get_parameter("relative_tolerance", rtol) call get_parameter("maximum_rejections", mrej) call get_parameter("chi" , chi) call get_parameter("factor" , fac ) call get_parameter("minimum_factor" , facmin) call get_parameter("maximum_factor" , facmax) ! select the integration method, check the correctness of the integration ! parameters and adjust the CFL coefficient if necessary ! select case(trim(integration)) case ("euler", "EULER", "rk1", "RK1") evolve => evolve_euler order = 1 registers = 2 name_int = "1st order Euler" case ("rk2", "RK2") evolve => evolve_rk2 order = 2 registers = 2 name_int = "2nd order Runge-Kutta" 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)" case ("rk3.5", "RK3.5", "ssprk(5,3)", "SSPRK(5,3)") evolve => evolve_ssprk35 order = 3 registers = 3 cfl = 2.65062919143939d+00 * cfl name_int = "3rd order SSPRK(5,3)" case ("rk4.10", "RK4.10", "ssprk(10,4)", "SSPRK(10,4)") evolve => evolve_ssprk4_10 order = 4 registers = 3 cfl = 6.0d+00 * cfl name_int = "Optimal 4th order SSPRK(10,4)" case ("ssprk(m,2)", "SSPRK(m,2)") error_control = .true. evolve => evolve_ssprk2_m order = 2 registers = 3 stages = max(2, min(9, stages)) cfl = (stages - 1) * cfl write(name_int, "('2ⁿᵈ-order ',i0,'-step embedded " // & "SSPRK(',i0,',2) method')") stages, stages allocate(c(stages), dl(stages), bt(stages), bh(stages)) do n = 1, stages c(n) = n - 1.0d+00 end do c (:) = c(:) / (stages - 1) bt(:) = 1.0d+00 / (stages - 1) bt(stages) = 1.0d+00 / stages bh(:) = 1.0d+00 / stages bh(1) = (stages + 1.0d+00) / stages**2 bh(stages) = (stages - 1.0d+00) / stages**2 dl(1) = (stages - 1.0d+00) / stages dl(2) = 1.0d+00 - dl(1) betas(:) = [ 5.8d-01, -2.1d-01, 1.0d-01 ] call get_parameter("beta(1)", betas(1)) call get_parameter("beta(2)", betas(2)) call get_parameter("beta(3)", betas(3)) betas(:) = - betas(:) / 2.0d+00 case ("SSPRK3(2)4", "SSP3(2)4", "ssprk3(2)4", "ssp3(2)4") error_control = .true. evolve => evolve_ssp324 order = 3 registers = 3 stages = 4 cfl = 2.0d+00 * cfl name_int = "3ʳᵈ-order 4-step embedded SSP3(2)4 method" betas(:) = [ 5.5d-01, -2.7d-01, 5.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 ("SSPRK3(2)m", "SSP3(2)m", "ssprk3(2)m", "ssp3(2)m", "SSPRK(m,3)", "ssprk(m,3)") error_control = .true. evolve => evolve_ssprk32m order = 3 registers = 4 n = 2 do while(n**2 <= stages) n = n + 1 end do n = n - 1 stages = max(4, n**2) cfl = (n - 1) * n * cfl write(name_int, "('3ʳᵈ-order order ',i0,'-step embedded SSPRK3(2)n² method')") stages betas(:) = [ 5.5d-01, -2.7d-01, 5.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 if (stages == 4 .and. verbose) then write(*,*) write(*,*) "ADVICE:" write(*,*) "The method 'SSPRK3(2)m' is primarily designed for n²," // & " where n > 2. Since you intend to use 4 stages with" // & " this method, it is highly advised to use method" // & " 'SSPRK3(2)4', since it slightly better optimized and" // & " significantly used less memory." end if 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+1), 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, & 1.000000000000000000000000000000000000d+00 ] 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+1), 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, & 1.000000000000000000000000000000000000d+00 ] 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+1), 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, & 1.000000000000000000000000000000000000d+00 ] 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+1), 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, & 1.000000000000000000000000000000000000d+00 ] 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+1), 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, & 1.000000000000000000000000000000000000d+00 ] 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+1), 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, & 1.000000000000000000000000000000000000d+00 ] 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 write(*,*) write(*,"(1x,a)") "ERROR!" write(*,"(1x,a)") "The selected time advance method is not " // & "implemented: " // trim(integration) write(*,"(1x,a)") "Available methods: 'euler' 'rk2', 'rk3'," // & " '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', 'rk3(2)5', 'rk4(3)9', 'rk5(4)10',"// & "'rk3(2)5f', 'rk4(3)9f', 'rk5(4)10f'." end if status = 1 end select ! select error norm ! call get_parameter("error_norm", error_norm) select case(trim(error_norm)) case("max", "maximum", "infinity") update_errors => update_errors_max name_norm = "maximum norm" case default update_errors => update_errors_l2 name_norm = "L2-norm" end select ! reset recent error history ! errs = 1.0d+00 ! GLM coefficients for all refinement levels ! allocate(adecay(toplev), aglm(toplev)) #if NDIMS == 3 aglm(:) = - glm_alpha / min(adx(:), ady(:), adz(:)) #else /* NDIMS == 3 */ aglm(:) = - glm_alpha / min(adx(:), ady(:)) #endif /* NDIMS == 3 */ #ifdef PROFILE ! stop accounting time for module initialization/finalization ! call stop_timer(imi) #endif /* PROFILE */ !------------------------------------------------------------------------------- ! end subroutine initialize_evolution ! !=============================================================================== ! ! subroutine FINALIZE_EVOLUTION: ! ----------------------------- ! ! Subroutine releases memory used by the module. ! ! Arguments: ! ! status - an integer flag for error return value; ! !=============================================================================== ! subroutine finalize_evolution(verbose, status) ! include external procedures ! use helpers , only : print_section, print_parameter use iso_fortran_env, only : error_unit ! local variables are not implicit by default ! implicit none ! subroutine arguments ! logical, intent(in) :: verbose integer, intent(out) :: status ! local parameters ! character(len=*), parameter :: loc = 'EVOLUTION::finalize_evolution()' ! !------------------------------------------------------------------------------- ! #ifdef PROFILE ! start accounting time for module initialization/finalization ! call start_timer(imi) #endif /* PROFILE */ 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 ! status = 0 ! nullify pointers ! nullify(evolve) nullify(update_errors) ! print integration summary ! if (niterations > 0) then call print_section(verbose, "Integration summary") call print_parameter(verbose, "Number of iterations", niterations) call print_parameter(verbose, "Number of rejections", nrejections) end if #ifdef PROFILE ! stop accounting time for module initialization/finalization ! call stop_timer(imi) #endif /* PROFILE */ !------------------------------------------------------------------------------- ! end subroutine finalize_evolution ! !=============================================================================== ! ! subroutine PRINT_EVOLUTION: ! -------------------------- ! ! Subroutine prints module parameters. ! ! Arguments: ! ! verbose - a logical flag turning the information printing; ! !=============================================================================== ! subroutine print_evolution(verbose) ! import external procedures and variables ! use equations, only : magnetized use helpers , only : print_section, print_parameter ! local variables are not implicit by default ! implicit none ! subroutine arguments ! logical, intent(in) :: verbose ! !------------------------------------------------------------------------------- ! if (verbose) then call print_section(verbose, "Evolution") call print_parameter(verbose, "time advance", name_int) call print_parameter(verbose, "no. of registers", registers) call print_parameter(verbose, "CFL coefficient", cfl) if (magnetized) then call print_parameter(verbose, "GLM alpha coefficient", glm_alpha) end if if (error_control) then call print_parameter(verbose, "absolute tolerance", atol) call print_parameter(verbose, "relative tolerance", rtol) call print_parameter(verbose, "error norm" , name_norm) call print_parameter(verbose, "maximum rejections", mrej) call print_parameter(verbose, "chi" , chi) call print_parameter(verbose, "factor" , fac) call print_parameter(verbose, "minimum factor" , facmin) call print_parameter(verbose, "maximum factor" , facmax) end if end if !------------------------------------------------------------------------------- ! end subroutine print_evolution ! !=============================================================================== ! ! subroutine ADVANCE: ! ------------------ ! ! Subroutine advances the solution by one time step using the selected time ! integration method. ! ! Arguments: ! ! status - the subroutine call status: 0 for success, otherwise failure; ! !=============================================================================== ! subroutine advance(status) ! references ! use blocks , only : set_blocks_update use coordinates, only : toplev use forcing , only : update_forcing, forcing_enabled use mesh , only : update_mesh ! local variables are not implicit by default ! implicit none ! input variables ! integer, intent(out) :: status ! !------------------------------------------------------------------------------- ! #ifdef PROFILE ! start accounting time for solution advance ! call start_timer(ima) #endif /* PROFILE */ ! advance the solution using the selected method ! call evolve() ! add forcing contribution ! if (forcing_enabled) call update_forcing(dt) ! check if we need to perform the refinement step ! if (toplev > 1) then ! set all meta blocks to not be updated ! call set_blocks_update(.false.) ! check refinement and refine ! call update_mesh(status) if (status /= 0) go to 100 ! update primitive variables ! call update_variables(time + dt, 0.0d+00, status) if (status /= 0) go to 100 ! set all meta blocks to be updated ! call set_blocks_update(.true.) end if ! toplev > 1 #ifdef DEBUG ! check variables for NaNs ! call check_variables() #endif /* DEBUG */ ! error entry point ! 100 continue #ifdef PROFILE ! stop accounting time for solution advance ! call stop_timer(ima) #endif /* PROFILE */ !------------------------------------------------------------------------------- ! end subroutine advance ! !=============================================================================== ! ! subroutine INITIALIZE_TIME_STEP: ! ------------------------------- ! ! Subroutine estimates the initial time step. ! ! 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 initialize_time_step() 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 #endif /* NDIMS == 3 */ use equations , only : nf, maxspeed, cmax, cmax2 use iso_fortran_env, only : error_unit use mesh , only : work, nwork, work_in_use #ifdef MPI use mpitools , only : reduce_maximum, reduce_sum #endif /* MPI */ use sources , only : viscosity, resistivity, update_sources implicit none type(block_data), pointer :: pdata logical :: flag integer :: mlev, n, status real(kind=8) :: dx_min, fnorm, h0, h1 real(kind=8), dimension(3) :: d real(kind=8), dimension(:,:,:,:), pointer :: sc, df real(kind=8), parameter :: eps = tiny(cmax) character(len=*), parameter :: loc = 'EVOLUTION:initialize_time_step()' ! !------------------------------------------------------------------------------- ! #ifdef PROFILE call start_timer(imi) #endif /* PROFILE */ cmax = eps mlev = 1 pdata => list_data do while (associated(pdata)) mlev = max(mlev, pdata%meta%level) cmax = max(cmax, maxspeed(pdata%q(:,:,:,:))) 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 */ ! determine the time step due to the stability ! dth = dx_min / max(cmax & + 2.0d+00 * max(viscosity, resistivity) / dx_min, eps) dt = cfl * dth ! determine the time step due to the error ! if (error_control) then status = 0 if (work_in_use) then write(error_unit,"('[',a,']: ',a,3i4,a)") trim(loc), & "Workspace is already occupied!" status = 1 else n = 2 * nf * nc**NDIMS if (n > nwork) then write(error_unit,"('[',a,']: ',a,3i4,a)") trim(loc), & "Workspace has too be increased!" nwork = n deallocate(work) allocate(work(nwork), stat=status) if (status /= 0) then write(error_unit,"('[',a,']: ',a)") trim(loc), & "Cannot increase the workspace's size!" go to 100 end if end if n = n / 2 #if NDIMS == 3 sc(1:nf,1:nc,1:nc,1:nc) => work( 1: n) df(1:nf,1:nc,1:nc,1:nc) => work(n+1:2*n) #else /* NDIMS == 3 */ sc(1:nf,1:nc,1:nc,1: 1) => work( 1: n) df(1:nf,1:nc,1:nc,1: 1) => work(n+1:2*n) #endif /* NDIMS == 3 */ end if 100 continue #ifdef MPI call reduce_maximum(status) #endif /* MPI */ if (status > 0) then write(error_unit,"('[',a,']: ',a)") trim(loc), & "Could not use the workspace!" go to 200 end if work_in_use = .true. 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 h0 = 1.0d-06 end if n = 10 flag = .true. do while(flag) 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, status) if (status /= 0) then h0 = 2.5d-01 * h0 flag = n > 0 n = n - 1 else flag = .false. end if end do if (status /= 0) then write(error_unit,"('[', a, ']: ', a)") trim(loc), & "Could not estimate the initial time step due to " // & "the error. Setting it to the CFL time step." dte = dt else 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 work_in_use = .false. 200 continue end if #ifdef PROFILE 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 */ !------------------------------------------------------------------------------- ! end subroutine new_time_step ! !=============================================================================== !! !!*** PRIVATE SUBROUTINES **************************************************** !! !=============================================================================== ! !=============================================================================== ! ! subroutine EVOLVE_EULER: ! ----------------------- ! ! Subroutine advances the solution by one time step using the 1st order ! Euler integration method. ! ! References: ! ! [1] Press, W. H, Teukolsky, S. A., Vetterling, W. T., Flannery, B. P., ! "Numerical Recipes in Fortran", ! Cambridge University Press, Cambridge, 1992 ! !=============================================================================== ! subroutine evolve_euler() use blocks , only : block_data, list_data use boundaries, only : boundary_fluxes use equations , only : ibp, cmax use sources , only : update_sources implicit none type(block_data), pointer :: pdata integer :: status real(kind=8) :: t ! !------------------------------------------------------------------------------- ! #ifdef PROFILE call start_timer(imu) #endif /* PROFILE */ status = 1 do while(status /= 0) t = time + dt pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) pdata%u => pdata%uu(:,:,:,:,2) pdata => pdata%next end do pdata => list_data do while (associated(pdata)) call update_increment(pdata) call update_sources(pdata, t, dt, pdata%du(:,:,:,:)) pdata => pdata%next end do call boundary_fluxes() pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + dt * pdata%du(:,:,:,:) pdata => pdata%next end do if (ibp > 0) then adecay(:) = exp(aglm(:) * cmax * dt) pdata => list_data do while (associated(pdata)) pdata%u(ibp,:,:,:) = adecay(pdata%meta%level) * pdata%u(ibp,:,:,:) pdata => pdata%next end do end if call update_variables(t, dt, status) if (status /= 0) dt = 2.5d-01 * dt end do pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2) pdata%u => pdata%uu(:,:,:,:,1) pdata => pdata%next end do #ifdef PROFILE call stop_timer(imu) #endif /* PROFILE */ !------------------------------------------------------------------------------- ! end subroutine evolve_euler ! !=============================================================================== ! ! subroutine EVOLVE_RK2: ! --------------------- ! ! Subroutine advances the solution by one time step using the 2nd order ! Runge-Kutta time integration method. ! ! References: ! ! [1] Press, W. H, Teukolsky, S. A., Vetterling, W. T., Flannery, B. P., ! "Numerical Recipes in Fortran", ! Cambridge University Press, Cambridge, 1992 ! !=============================================================================== ! subroutine evolve_rk2() use blocks , only : block_data, list_data use boundaries, only : boundary_fluxes use equations , only : ibp, cmax use sources , only : update_sources implicit none type(block_data), pointer :: pdata integer :: status real(kind=8) :: t, ds ! !------------------------------------------------------------------------------- ! #ifdef PROFILE call start_timer(imu) #endif /* PROFILE */ status = 1 do while(status /= 0) ds = 5.0d-01 * dt != 1st step: U(1) = U(n) + dt * F[U(n)] ! t = time + dt pdata => list_data do while (associated(pdata)) call update_increment(pdata) call update_sources(pdata, t, 0.0d+00, pdata%du(:,:,:,:)) pdata => pdata%next end do call boundary_fluxes() pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) + dt * pdata%du(:,:,:,:) pdata%u => pdata%uu(:,:,:,:,2) pdata => pdata%next end do call update_variables(t, dt, status) if (status /= 0) go to 100 != 2nd step: U(n+1) = 1/2 U(n) + 1/2 U(1) + 1/2 dt * F[U(1)] ! pdata => list_data do while (associated(pdata)) call update_increment(pdata) call update_sources(pdata, t, ds, pdata%du(:,:,:,:)) pdata => pdata%next end do call boundary_fluxes() pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,2) = 5.0d-01 * (pdata%uu(:,:,:,:,1) & + pdata%uu(:,:,:,:,2)) & + ds * pdata%du(:,:,:,:) pdata => pdata%next end do if (ibp > 0) then adecay(:) = exp(aglm(:) * cmax * dt) pdata => list_data do while (associated(pdata)) pdata%u(ibp,:,:,:) = adecay(pdata%meta%level) * pdata%u(ibp,:,:,:) pdata => pdata%next end do end if call update_variables(t, ds, status) 100 continue if (status /= 0) dt = 2.5d-01 * dt end do pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2) pdata%u => pdata%uu(:,:,:,:,1) pdata => pdata%next end do #ifdef PROFILE call stop_timer(imu) #endif /* PROFILE */ !------------------------------------------------------------------------------- ! end subroutine evolve_rk2 ! !=============================================================================== ! ! subroutine EVOLVE_RK3: ! --------------------- ! ! Subroutine advances the solution by one time step using the 3rd order ! Runge-Kutta time integration method. ! ! References: ! ! [1] Press, W. H, Teukolsky, S. A., Vetterling, W. T., Flannery, B. P., ! "Numerical Recipes in Fortran", ! Cambridge University Press, Cambridge, 1992 ! !=============================================================================== ! subroutine evolve_rk3() use blocks , only : block_data, list_data use boundaries, only : boundary_fluxes use equations , only : ibp, cmax use sources , only : update_sources implicit none type(block_data), pointer :: pdata integer :: status real(kind=8) :: t, ds real(kind=8), parameter :: f21 = 3.0d+00 / 4.0d+00, f22 = 1.0d+00 / 4.0d+00 real(kind=8), parameter :: f31 = 1.0d+00 / 3.0d+00, f32 = 2.0d+00 / 3.0d+00 ! !------------------------------------------------------------------------------- ! #ifdef PROFILE call start_timer(imu) #endif /* PROFILE */ status = 1 do while(status /= 0) != 1st substep: U(1) = U(n) + dt F[U(n)] ! t = time + dt pdata => list_data do while (associated(pdata)) call update_increment(pdata) call update_sources(pdata, t, 0.0d+00, pdata%du(:,:,:,:)) pdata => pdata%next end do call boundary_fluxes() pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) + dt * pdata%du(:,:,:,:) pdata%u => pdata%uu(:,:,:,:,2) pdata => pdata%next end do call update_variables(t, dt, status) if (status /= 0) go to 100 != 2nd step: U(2) = 3/4 U(n) + 1/4 U(1) + 1/4 dt F[U(1)] ! ds = f22 * dt t = time + 0.5d+00 * dt pdata => list_data do while (associated(pdata)) call update_increment(pdata) call update_sources(pdata, t, ds, pdata%du(:,:,:,:)) pdata => pdata%next end do call boundary_fluxes() pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,2) = f21 * pdata%uu(:,:,:,:,1) & + f22 * pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:) pdata => pdata%next end do call update_variables(t, ds, status) if (status /= 0) go to 100 != 3rd step: U(n+1) = 1/3 U(n) + 2/3 U(2) + 2/3 dt F[U(2)] ! ds = f32 * dt t = time + dt pdata => list_data do while (associated(pdata)) call update_increment(pdata) call update_sources(pdata, t, ds, pdata%du(:,:,:,:)) pdata => pdata%next end do call boundary_fluxes() pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,2) = f31 * pdata%uu(:,:,:,:,1) & + f32 * pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:) pdata%u => pdata%uu(:,:,:,:,1) pdata => pdata%next end do if (ibp > 0) then adecay(:) = exp(aglm(:) * cmax * dt) pdata => list_data do while (associated(pdata)) pdata%u(ibp,:,:,:) = adecay(pdata%meta%level) * pdata%u(ibp,:,:,:) pdata => pdata%next end do end if call update_variables(t, dt, status) 100 continue if (status /= 0) dt = 2.5d-01 * dt end do pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2) pdata%u => pdata%uu(:,:,:,:,1) pdata => pdata%next end do #ifdef PROFILE call stop_timer(imu) #endif /* PROFILE */ !------------------------------------------------------------------------------- ! end subroutine evolve_rk3 ! !=============================================================================== ! ! subroutine EVOLVE_SSPRK34: ! ------------------------- ! ! Subroutine advances the solution by one time step using the 3rd order ! 4-stage Strong Stability Preserving Runge-Kutta time integration method. ! ! References: ! ! [1] Ruuth, S. J., ! "Global Optimization of Explicit Strong-Stability-Preserving ! Runge-Kutta methods", ! Mathematics of Computation, ! 2006, vol. 75, no. 253, pp. 183-207 ! !=============================================================================== ! subroutine evolve_ssprk34() use blocks , only : block_data, list_data use boundaries, only : boundary_fluxes use equations , only : ibp, cmax use sources , only : update_sources implicit none type(block_data), pointer :: pdata integer :: status real(kind=8) :: t, ds ! !------------------------------------------------------------------------------- ! #ifdef PROFILE call start_timer(imu) #endif /* PROFILE */ status = 1 do while(status /= 0) != 1st step: U(1) = U(n) + 1/2 dt F[U(n)] ! ds = 5.0d-01 * dt t = time + ds pdata => list_data do while (associated(pdata)) call update_increment(pdata) call update_sources(pdata, t, 0.0d+00, pdata%du(:,:,:,:)) pdata => pdata%next end do call boundary_fluxes() pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) + ds * pdata%du(:,:,:,:) pdata%u => pdata%uu(:,:,:,:,2) pdata => pdata%next end do call update_variables(t, ds, status) if (status /= 0) go to 100 != 2nd step: U(2) = U(2) + 1/2 dt F[U(1)] ! t = time + dt pdata => list_data do while (associated(pdata)) call update_increment(pdata) call update_sources(pdata, t, ds, pdata%du(:,:,:,:)) pdata => pdata%next end do call boundary_fluxes() pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:) pdata => pdata%next end do call update_variables(t, ds, status) if (status /= 0) go to 100 != 3rd step: U(3) = 2/3 U(n) + 1/3 (U(2) + 1/2 dt F[U(2)]) ! t = time + ds pdata => list_data do while (associated(pdata)) call update_increment(pdata) call update_sources(pdata, t, ds, pdata%du(:,:,:,:)) pdata => pdata%next end do call boundary_fluxes() pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,2) = (2.0d+00 * pdata%uu(:,:,:,:,1) & + pdata%uu(:,:,:,:,2) & + ds * pdata%du(:,:,:,:)) / 3.0d+00 pdata => pdata%next end do call update_variables(t, ds, status) if (status /= 0) go to 100 != the final step: U(n+1) = U(3) + 1/2 dt F[U(3)] ! t = time + dt pdata => list_data do while (associated(pdata)) call update_increment(pdata) call update_sources(pdata, t, ds, pdata%du(:,:,:,:)) pdata => pdata%next end do call boundary_fluxes() pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:) pdata => pdata%next end do if (ibp > 0) then adecay(:) = exp(aglm(:) * cmax * dt) pdata => list_data do while (associated(pdata)) pdata%u(ibp,:,:,:) = adecay(pdata%meta%level) * pdata%u(ibp,:,:,:) pdata => pdata%next end do end if call update_variables(t, dt, status) 100 continue if (status /= 0) dt = 2.5d-01 * dt end do pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2) pdata%u => pdata%uu(:,:,:,:,1) pdata => pdata%next end do #ifdef PROFILE call stop_timer(imu) #endif /* PROFILE */ !------------------------------------------------------------------------------- ! end subroutine evolve_ssprk34 ! !=============================================================================== ! ! subroutine EVOLVE_SSPRK35: ! ------------------------- ! ! Subroutine advances the solution by one time step using the 3rd order ! 5-stage Strong Stability Preserving Runge-Kutta time integration method. ! ! References: ! ! [1] Ruuth, S. J., ! "Global Optimization of Explicit Strong-Stability-Preserving ! Runge-Kutta methods", ! Mathematics of Computation, ! 2006, vol. 75, no. 253, pp. 183-207 ! !=============================================================================== ! subroutine evolve_ssprk35() use blocks , only : block_data, list_data use boundaries, only : boundary_fluxes use equations , only : ibp, cmax use sources , only : update_sources implicit none type(block_data), pointer :: pdata integer :: status real(kind=8) :: t, ds real(kind=8), parameter :: b1 = 3.77268915331368d-01 real(kind=8), parameter :: b3 = 2.42995220537396d-01 real(kind=8), parameter :: b4 = 2.38458932846290d-01 real(kind=8), parameter :: b5 = 2.87632146308408d-01 real(kind=8), parameter :: a31 = 3.55909775063327d-01 real(kind=8), parameter :: a33 = 6.44090224936673d-01 real(kind=8), parameter :: a41 = 3.67933791638137d-01 real(kind=8), parameter :: a44 = 6.32066208361863d-01 real(kind=8), parameter :: a53 = 2.37593836598569d-01 real(kind=8), parameter :: a55 = 7.62406163401431d-01 ! !------------------------------------------------------------------------------- ! #ifdef PROFILE call start_timer(imu) #endif /* PROFILE */ status = 1 do while(status /= 0) != 1st step: U(1) = U(n) + b1 dt F[U(n)] ! ds = b1 * dt t = time + ds pdata => list_data do while (associated(pdata)) call update_increment(pdata) call update_sources(pdata, t, 0.0d+00, pdata%du(:,:,:,:)) pdata => pdata%next end do call boundary_fluxes() pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) + ds * pdata%du(:,:,:,:) pdata%u => pdata%uu(:,:,:,:,2) pdata => pdata%next end do call update_variables(t, ds, status) if (status /= 0) go to 100 != 2nd step: U(2) = U(1) + b1 dt F[U(1)] ! t = time + 2.0d+00 * ds pdata => list_data do while (associated(pdata)) call update_increment(pdata) call update_sources(pdata, t, ds, pdata%du(:,:,:,:)) pdata => pdata%next end do call boundary_fluxes() pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:) pdata => pdata%next end do call update_variables(t, ds, status) if (status /= 0) go to 100 != 3rd step: U(3) = a31 U(n) + a33 U(2) + b3 dt F[U(2)] ! ds = b3 * dt t = time + (2.0d+00 * a33 * b1 + b3) * dt pdata => list_data do while (associated(pdata)) call update_increment(pdata) call update_sources(pdata, t, ds, pdata%du(:,:,:,:)) pdata => pdata%next end do call boundary_fluxes() pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,2) = a31 * pdata%uu(:,:,:,:,1) & + a33 * pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:) pdata => pdata%next end do call update_variables(t, ds, status) if (status /= 0) go to 100 != 4th step: U(4) = a41 U(n) + a44 U(3) + b4 dt F[U(3)] ! ds = b4 * dt t = time + ((2.0d+00 * b1 * a33 + b3) * a44 + b4) * dt pdata => list_data do while (associated(pdata)) call update_increment(pdata) call update_sources(pdata, t, ds, pdata%du(:,:,:,:)) pdata => pdata%next end do call boundary_fluxes() pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,3) = a41 * pdata%uu(:,:,:,:,1) & + a44 * pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:) pdata%u => pdata%uu(:,:,:,:,3) pdata => pdata%next end do call update_variables(t, ds, status) if (status /= 0) go to 100 != the final step: U(n+1) = a53 U(2) + a55 U(4) + b5 dt F[U(4)] ! ds = b5 * dt t = time + dt pdata => list_data do while (associated(pdata)) call update_increment(pdata) call update_sources(pdata, t, ds, pdata%du(:,:,:,:)) pdata => pdata%next end do call boundary_fluxes() pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,2) = a53 * pdata%uu(:,:,:,:,2) & + a55 * pdata%uu(:,:,:,:,3) + ds * pdata%du(:,:,:,:) pdata%u => pdata%uu(:,:,:,:,2) pdata => pdata%next end do if (ibp > 0) then adecay(:) = exp(aglm(:) * cmax * dt) pdata => list_data do while (associated(pdata)) pdata%u(ibp,:,:,:) = adecay(pdata%meta%level) * pdata%u(ibp,:,:,:) pdata => pdata%next end do end if call update_variables(t, dt, status) 100 continue if (status /= 0) dt = 2.5d-01 * dt end do pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2) pdata%u => pdata%uu(:,:,:,:,1) pdata => pdata%next end do #ifdef PROFILE call stop_timer(imu) #endif /* PROFILE */ !------------------------------------------------------------------------------- ! end subroutine evolve_ssprk35 ! !=============================================================================== ! ! subroutine EVOLVE_SSPRK4_10: ! --------------------------- ! ! Subroutine advances the solution by one time step using the 4rd order ! 10-stage Strong Stability Preserving Runge-Kutta time integration method. ! ! References: ! ! [1] Gottlieb, S., Ketcheson, D., Shu C. - W., ! "Strong stability preserving Runge-Kutta and multistep ! time discretizations", ! World Scientific Publishing, 2011 ! !=============================================================================== ! subroutine evolve_ssprk4_10() use blocks , only : block_data, list_data use boundaries, only : boundary_fluxes use equations , only : ibp, cmax use sources , only : update_sources implicit none type(block_data), pointer :: pdata integer :: i, status real(kind=8) :: t real(kind=8), dimension(9) :: ds real(kind=8), dimension(9), parameter :: c = & (/ 1.0d+00, 2.0d+00, 3.0d+00, 4.0d+00, 2.0d+00, & 3.0d+00, 4.0d+00, 5.0d+00, 6.0d+00 /) / 6.0d+00 ! !------------------------------------------------------------------------------- ! #ifdef PROFILE call start_timer(imu) #endif /* PROFILE */ ds(:) = c(:) * dt status = 1 do while(status /= 0) != 1st step: U(2) = U(1) ! pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) pdata%uu(:,:,:,:,3) = pdata%uu(:,:,:,:,1) pdata%u => pdata%uu(:,:,:,:,2) pdata => pdata%next end do != 2nd step: U(1) = [1 + dt/6 L] U(1), for i = 1, ..., 5 ! do i = 1, 5 t = time + ds(i) pdata => list_data do while (associated(pdata)) call update_increment(pdata) call update_sources(pdata, t, 0.0d+00, pdata%du(:,:,:,:)) pdata => pdata%next end do call boundary_fluxes() pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + ds(1) * pdata%du(:,:,:,:) pdata => pdata%next end do call update_variables(t, ds(1), status) if (status /= 0) go to 100 end do ! i = 1, 5 != 3rd step: U(2) = U(2)/25 + 9/25 U(1), U(1) = 15 U(2) - 5 U(1) ! pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,3) = (pdata%uu(:,:,:,:,3) & + 9.0d+00 * pdata%uu(:,:,:,:,2)) / 2.5d+01 pdata%uu(:,:,:,:,2) = 1.5d+01 * pdata%uu(:,:,:,:,3) & - 5.0d+00 * pdata%uu(:,:,:,:,2) pdata => pdata%next end do != 4th step: U(1) = [1 + dt/6 L] U(1), for i = 6, ..., 9 ! ! integrate the intermediate steps ! do i = 6, 9 t = time + ds(i) pdata => list_data do while (associated(pdata)) call update_increment(pdata) call update_sources(pdata, t, ds(1), pdata%du(:,:,:,:)) pdata => pdata%next end do call boundary_fluxes() pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + ds(1) * pdata%du(:,:,:,:) pdata => pdata%next end do call update_variables(t, ds(1), status) if (status /= 0) go to 100 end do ! i = 6, 9 != the final step: U(n+1) = U(2) + 3/5 U(1) + 1/10 dt F[U(1)] ! t = time + dt pdata => list_data do while (associated(pdata)) call update_increment(pdata) call update_sources(pdata, t, 1.0d-01 * dt, pdata%du(:,:,:,:)) pdata => pdata%next end do call boundary_fluxes() pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,3) & + 6.0d-01 * pdata%uu(:,:,:,:,2) & + 1.0d-01 * dt * pdata%du(:,:,:,:) pdata => pdata%next end do if (ibp > 0) then adecay(:) = exp(aglm(:) * cmax * dt) pdata => list_data do while (associated(pdata)) pdata%u(ibp,:,:,:) = adecay(pdata%meta%level) * pdata%u(ibp,:,:,:) pdata => pdata%next end do end if call update_variables(t, dt, status) 100 continue if (status /= 0) dt = 2.5d-01 * dt end do pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2) pdata%u => pdata%uu(:,:,:,:,1) pdata => pdata%next end do #ifdef PROFILE call stop_timer(imu) #endif /* PROFILE */ !------------------------------------------------------------------------------- ! end subroutine evolve_ssprk4_10 ! !=============================================================================== ! ! EMBEDDED METHODS ! !=============================================================================== ! ! subroutine EVOLVE_SSPRK2_M: ! -------------------------- ! ! Subroutine advances the solution by one time step using the 2nd order ! m-stage Strong Stability Preserving Runge-Kutta time integration method. ! Up to 9 stages are allowed, due to stability problems with more stages. ! ! References: ! ! [1] Conde, S., Fekete, I., Shadid, J. N., ! "Embedded error estimation and adaptive step-size control for ! optimal explicit strong stability preserving Runge–Kutta methods" ! 2018, arXiv:1806.08693v1 ! [2] Gottlieb, S. and Gottlieb, L.-A., J. ! "Strong Stability Preserving Properties of Runge-Kutta Time ! Discretization Methods for Linear Constant Coefficient Operators", ! Journal of Scientific Computing, ! 2003, vol. 18, no. 1, pp. 83-109 ! !=============================================================================== ! subroutine evolve_ssprk2_m() use blocks , only : block_data, list_data use boundaries , only : boundary_fluxes use coordinates , only : nb, ne use equations , only : errors, ibp, cmax use iso_fortran_env, only : error_unit use sources , only : update_sources implicit none type(block_data), pointer :: pdata logical :: test integer :: nrej, i, status real(kind=8) :: tm, dtm, dh, fc character(len=*), parameter :: loc = 'EVOLUTION:evolve_ssprk2_m()' !------------------------------------------------------------------------------- ! #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 2nd order solution, U(1) ! pdata%uu(:,:,:,:,3) - the new 1st order solution, U(2) ! 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%u => pdata%uu(:,:,:,:,2) pdata => pdata%next end do != 1st step: U(1) = U(1) + dt/(m-1) F[U(1)], for i = 1, ..., m-1 ! do i = 1, stages - 1 dtm = c(i) * dt tm = time + dtm 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(:,:,:,:,2) = pdata%uu(:,:,:,:,2) & + bt(i) * dt * pdata%du(:,:,:,:) pdata%uu(:,:,:,:,3) = pdata%uu(:,:,:,:,3) & + bh(i) * dt * pdata%du(:,:,:,:) pdata => pdata%next end do call update_variables(tm, dtm, status) if (status /= 0) go to 100 end do ! i = 1, stages - 1 != 2nd step: U(1) = (m-1)/m U(1) + 1/m U(2) + dt/m F[U(1)] ! i = stages dtm = c(i) * dt tm = time + dtm 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(:,:,:,:,2) = dl(1) * pdata%uu(:,:,:,:,2) & + dl(2) * pdata%uu(:,:,:,:,1) & + bt(i) * dt * pdata%du(:,:,:,:) pdata%uu(:,:,:,:,3) = pdata%uu(:,:,:,:,3) & + bh(i) * dt * pdata%du(:,:,:,:) pdata => pdata%next end do call update_variables(tm, dtm, status) if (status /= 0) go to 100 != 3rd 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,:,:,:,3) = adecay(pdata%meta%level) & * pdata%uu(ibp,:,:,:,3) pdata => pdata%next end do end if ! update the vector of errors ! call update_errors(2, 3) ! 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 100 continue if ((fc > fac .or. nrej >= mrej) .and. status == 0) then test = .false. errs(3) = errs(2) errs(2) = errs(1) else if (status == 0) then dt = dte nrej = nrej + 1 ! rejection count in the current step else dt = 2.5d-01 * dt end if 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, status) if (status /= 0) then write(error_unit,"('[', a, ']: ', a)") trim(loc), & "Possible bug: the revert to the previous state " //& "found it unphysical." end if end if niterations = niterations + 1 end do != final step: U(n+1) = U(1) - update the accepted solution ! pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2) pdata%u => pdata%uu(:,:,:,:,1) pdata => pdata%next end do #ifdef PROFILE call stop_timer(imu) #endif /* PROFILE */ !------------------------------------------------------------------------------- ! end subroutine evolve_ssprk2_m ! !=============================================================================== ! ! subroutine EVOLVE_SSP324: ! -------------------------- ! ! Subroutine advances the solution by one time step using the 3ʳᵈ-order ! 4-stage embedded Strong Stability Preserving Runge-Kutta time integration ! method SSP3(2)4 with the error control. ! ! 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.06836v2 ! [2] Conde, S., Fekete, I., Shadid, J. N., ! "Embedded error estimation and adaptive step-size control for ! optimal explicit strong stability preserving Runge–Kutta methods" ! 2018, arXiv:1806.08693v1 ! [3] Gottlieb, S., Ketcheson, D., Shu, C.-W., ! "Strong stability preserving Runge-Kutta and multistep time ! discretizations", ! 2011, World Scientific Publishing ! !=============================================================================== ! subroutine evolve_ssp324() use blocks , only : block_data, list_data use boundaries , only : boundary_fluxes use coordinates , only : nb, ne use equations , only : errors, ibp, cmax use iso_fortran_env, only : error_unit use sources , only : update_sources implicit none type(block_data), pointer :: pdata logical :: test integer :: nrej, i, status real(kind=8) :: tm, dtm, dh, fc real(kind=8), parameter :: onethird = 1.0d+00 / 3.0d+00, & twothird = 2.0d+00 / 3.0d+00 character(len=*), parameter :: loc = 'EVOLUTION:evolve_ssp324()' !------------------------------------------------------------------------------- ! #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 new 2nd order solution, U(2) ! do while(test) ! initiate the fractional time step ! dh = 5.0d-01 * dt != preparation step: U(1) = U(n) ! pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) pdata%u => pdata%uu(:,:,:,:,2) pdata => pdata%next end do != 1st to 3rd steps: U(1) = U(1) + ½ dt F[U(1)], for i = 1...3 ! do i = 1, 3 tm = time + (i - 1) * dh dtm = dh 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(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + dh * pdata%du(:,:,:,:) pdata => pdata%next end do call update_variables(tm, dtm, status) if (status /= 0) go to 100 end do ! i = 1, 3 != 4th step: U(1) = ⅔ U(n) + ⅓ U(1), U(2) = ⅓ U(n) + ⅔ U(1) ! pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,3) = onethird * pdata%uu(:,:,:,:,1) & + twothird * pdata%uu(:,:,:,:,2) pdata%uu(:,:,:,:,2) = twothird * pdata%uu(:,:,:,:,1) & + onethird * pdata%uu(:,:,:,:,2) pdata => pdata%next end do call update_variables(tm, dtm, status) if (status /= 0) go to 100 != 5th step: U(1) = U(1) + ½ dt F[U(1)] <- 3ʳᵈ-order candidate ! tm = time + dh dtm = dh 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(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + dh * pdata%du(:,:,:,:) pdata => pdata%next end do call update_variables(tm, dtm, status) if (status /= 0) go to 100 != 6th step: U(2) = ½ (U(1) + U(2)) <- 2ⁿᵈ-order approximation ! pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,3) = 5.0d-01 * (pdata%uu(:,:,:,:,3) & + pdata%uu(:,:,:,:,2)) pdata => pdata%next end do != 7th 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,:,:,:,3) = adecay(pdata%meta%level) & * pdata%uu(ibp,:,:,:,3) pdata => pdata%next end do end if ! update the vector of errors ! call update_errors(2, 3) ! 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 100 continue if ((fc > fac .or. nrej >= mrej) .and. status == 0) then test = .false. errs(3) = errs(2) errs(2) = errs(1) else if (status == 0) then dt = dte nrej = nrej + 1 ! rejection count in the current step else dt = 2.5d-01 * dt end if 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, status) if (status /= 0) then write(error_unit,"('[', a, ']: ', a)") trim(loc), & "Possible bug: the revert to the previous state " //& "found it unphysical." end if end if niterations = niterations + 1 end do != final step: U(n+1) = U(1) - update the accepted solution ! pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2) pdata%u => pdata%uu(:,:,:,:,1) pdata => pdata%next end do #ifdef PROFILE call stop_timer(imu) #endif /* PROFILE */ !------------------------------------------------------------------------------- ! end subroutine evolve_ssp324 ! !=============================================================================== ! ! subroutine EVOLVE_SSPRK32M: ! -------------------------- ! ! Subroutine advances the solution by one time step using the 3ʳᵈ-order ! m-stage (m=n²) embedded Strong Stability Preserving Runge-Kutta time ! integration method SSP3(2)m with the error control. ! ! References: ! ! [1] Conde, S., Fekete, I., Shadid, J. N., ! "Embedded error estimation and adaptive step-size control for ! optimal explicit strong stability preserving Runge–Kutta methods" ! 2018, arXiv:1806.08693v1 ! [2] Gottlieb, S., Ketcheson, D., Shu, C.-W., ! "Strong stability preserving Runge-Kutta and multistep time ! discretizations", ! 2011, World Scientific Publishing ! !=============================================================================== ! subroutine evolve_ssprk32m() use blocks , only : block_data, list_data use boundaries , only : boundary_fluxes use equations , only : errors, ibp, cmax use iso_fortran_env, only : error_unit use sources , only : update_sources implicit none type(block_data), pointer :: pdata logical , save :: first = .true. integer , save :: i1, i2, i3, i4, n, m real(kind=8), save :: f1, f2, f3, f4 logical :: test integer :: nrej, i, status real(kind=8) :: tm, dtm, dh, fc character(len=*), parameter :: loc = 'EVOLUTION:evolve_ssprk32m()' !------------------------------------------------------------------------------- ! #ifdef PROFILE call start_timer(imu) #endif /* PROFILE */ if (first) then n = 2 do while(n**2 < stages) n = n + 1 end do m = stages - n i1 = (n - 1) * (n - 2) / 2 i2 = i1 + 1 i3 = n * (n + 1) / 2 i4 = i3 + 1 fc = real(2 * n - 1, kind=8) f1 = real(n - 1, kind=8) / fc f2 = real(n , kind=8) / fc f4 = real(m , kind=8) / stages f3 = 1.0d+00 - f4 first = .false. end if 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) ! initiate the fractional time step ! dh = dt / m != preparation step: U(1) = U(n) ! pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) pdata%u => pdata%uu(:,:,:,:,2) pdata => pdata%next end do != 1st step: U(1) = U(1) + dt/r F[U(1)], for i = 1...(n-1)*(n-2)/2 ! do i = 1, i1 tm = time + (i - 1) * dh dtm = dh 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(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + dh * pdata%du(:,:,:,:) pdata => pdata%next end do call update_variables(tm, dtm, status) if (status /= 0) go to 100 end do ! i = 1, i1 != 2nd step: U(2) = U(1) ! pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,3) = pdata%uu(:,:,:,:,2) pdata => pdata%next end do != 3rd step: U(1) = U(1) + dt/r F[U(1)], for i = (n-1)*(n-2)/2+1...n*(n+1)/2 ! do i = i2, i3 tm = time + (i - 1) * dh dtm = dh 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(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + dh * pdata%du(:,:,:,:) pdata => pdata%next end do call update_variables(tm, dtm, status) if (status /= 0) go to 100 end do ! i = i2, i3 != 4th step: U(3) = (1 - r/s) * U(n) + r/s * U(1), ! U(1) = [(n - 1) * U(1) + n * U(2)] / (2*n - 1), ! U(3) = U(3) - r/s * U(1) ! pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,4) = f3 * pdata%uu(:,:,:,:,1) & + f4 * pdata%uu(:,:,:,:,2) pdata%uu(:,:,:,:,2) = f1 * pdata%uu(:,:,:,:,2) & + f2 * pdata%uu(:,:,:,:,3) pdata%uu(:,:,:,:,4) = pdata%uu(:,:,:,:,4) & - f4 * pdata%uu(:,:,:,:,2) pdata => pdata%next end do call update_variables(tm, dtm, status) if (status /= 0) go to 100 != 6th step: U(1) = U(1) + dt/r F[U(1)], for i = n*(n+1)/2+1...stages ! do i = i4, stages tm = time + (i - 1 - n) * dh dtm = dh 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(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + dh * pdata%du(:,:,:,:) pdata => pdata%next end do call update_variables(tm, dtm, status) if (status /= 0) go to 100 end do ! i = i4, stages != 7th step: U(3) = U(3) + r/s * U(1) ! pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,4) = pdata%uu(:,:,:,:,4) + f4 * pdata%uu(:,:,:,:,2) pdata => pdata%next end do != 8th 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 ! update the vector of errors ! 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 100 continue if ((fc > fac .or. nrej >= mrej) .and. status == 0) then test = .false. errs(3) = errs(2) errs(2) = errs(1) else if (status == 0) then dt = dte nrej = nrej + 1 ! rejection count in the current step else dt = 2.5d-01 * dt end if 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, status) if (status /= 0) then write(error_unit,"('[', a, ']: ', a)") trim(loc), & "Possible bug: the revert to the previous state " //& "found it unphysical." end if end if niterations = niterations + 1 end do != final step: U(n+1) = U(1) - update the accepted solution ! pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2) pdata%u => pdata%uu(:,:,:,:,1) pdata => pdata%next end do #ifdef PROFILE call stop_timer(imu) #endif /* PROFILE */ !------------------------------------------------------------------------------- ! end subroutine evolve_ssprk32m ! !=============================================================================== ! ! 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 iso_fortran_env, only : error_unit use sources , only : update_sources implicit none type(block_data), pointer :: pdata logical :: test integer :: i, l, nrej, status real(kind=8) :: tm, dtm, fc character(len=*), parameter :: loc = 'EVOLUTION:evolve_3sstarplus()' !------------------------------------------------------------------------------- ! #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 dtm = c(2) * dt call update_variables(time, dtm, status) if (status /= 0) go to 100 != remaining stages ! do i = 2, stages dtm = c(i) * dt tm = time + dtm 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 dtm = c(i+1) * dt call update_variables(tm, dtm, status) if (status /= 0) go to 100 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 100 continue if ((fc > fac .or. nrej >= mrej) .and. status == 0) then test = .false. errs(3) = errs(2) errs(2) = errs(1) if (fsal) fsal_update = .true. else if (status == 0) then dt = dte nrej = nrej + 1 ! rejection count in the current step else dt = 2.5d-01 * dt end if 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, status) if (status /= 0) then write(error_unit,"('[', a, ']: ', a)") trim(loc), & "Possible bug: the revert to the previous state " //& "found it unphysical." end if if (fsal) fsal_update = .false. end if niterations = niterations + 1 end do != final step: U(n+1) = U(1) ! pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2) pdata%u => pdata%uu(:,:,:,:,1) pdata => pdata%next end do tm = time + dt call update_variables(tm, dt, status) if (status /= 0) then write(error_unit,"('[', a, ']: ', a)") trim(loc), & "Possible bug: the new state has been accepted, " // & "nevertheless the variable update found it unphysical." end if #ifdef PROFILE call stop_timer(imu) #endif /* PROFILE */ !------------------------------------------------------------------------------- ! end subroutine evolve_3sstarplus ! !=============================================================================== ! ! subroutine UPDATE_INCREMENT: ! --------------------------- ! ! Subroutine calculates the conservative variable increment from ! directional fluxes. ! ! Arguments: ! ! pdata - the point to data block storing the directional fluxes; ! !=============================================================================== ! subroutine update_increment(pdata) ! include external variables ! use blocks , only : block_data use coordinates, only : nbl, ne, neu, nn => bcells use coordinates, only : adx, ady, adxi, adyi #if NDIMS == 3 use coordinates, only : adz, adzi #endif /* NDIMS == 3 */ use equations , only : nf, ns use equations , only : idn, isl, isu use schemes , only : update_flux ! local variables are not implicit by default ! implicit none ! subroutine arguments ! type(block_data), pointer, intent(inout) :: pdata ! local variables ! integer :: i , j , k = 1, p integer :: im1, ip1 integer :: jm1, jp1 #if NDIMS == 3 integer :: km1, kp1 #endif /* NDIMS == 3 */ real(kind=8) :: df real(kind=8), dimension(NDIMS) :: dh, dhi #if NDIMS == 3 real(kind=8), dimension(nf,nn,nn,nn,3) :: f #else /* NDIMS == 3 */ real(kind=8), dimension(nf,nn,nn, 1,2) :: f #endif /* NDIMS == 3 */ ! !------------------------------------------------------------------------------- ! #ifdef PROFILE ! start accounting time for the increment update ! call start_timer(iui) #endif /* PROFILE */ ! prepare the coordinate intervals ! dh(1) = adx(pdata%meta%level) dh(2) = ady(pdata%meta%level) #if NDIMS == 3 dh(3) = adz(pdata%meta%level) #endif /* NDIMS == 3 */ dhi(1) = adxi(pdata%meta%level) dhi(2) = adyi(pdata%meta%level) #if NDIMS == 3 dhi(3) = adzi(pdata%meta%level) #endif /* NDIMS == 3 */ ! calculate the interface fluxes ! call update_flux(dh(1:NDIMS), pdata%q(:,:,:,:), f(:,:,:,:,1:NDIMS)) ! store block interface fluxes ! pdata%fx(:,1,:,:) = f(:,nbl,:,:,1) pdata%fx(:,2,:,:) = f(:,ne ,:,:,1) pdata%fy(:,:,1,:) = f(:,:,nbl,:,2) pdata%fy(:,:,2,:) = f(:,:,ne ,:,2) #if NDIMS == 3 pdata%fz(:,:,:,1) = f(:,:,:,nbl,3) pdata%fz(:,:,:,2) = f(:,:,:,ne ,3) #endif /* NDIMS == 3 */ ! calculate the variable update from the directional fluxes ! #if NDIMS == 3 do k = nbl, neu km1 = k - 1 #endif /* NDIMS == 3 */ do j = nbl, neu jm1 = j - 1 do i = nbl, neu im1 = i - 1 #if NDIMS == 3 pdata%du(1:nf,i,j,k) = - dhi(1) * (f(:,i,j,k,1) - f(:,im1,j,k,1)) & - dhi(2) * (f(:,i,j,k,2) - f(:,i,jm1,k,2)) & - dhi(3) * (f(:,i,j,k,3) - f(:,i,j,km1,3)) #else /* NDIMS == 3 */ pdata%du(1:nf,i,j,k) = - dhi(1) * (f(:,i,j,k,1) - f(:,im1,j,k,1)) & - dhi(2) * (f(:,i,j,k,2) - f(:,i,jm1,k,2)) #endif /* NDIMS == 3 */ end do ! i = nbl, neu end do ! j = nbl, neu #if NDIMS == 3 end do ! k = nbl, neu #endif /* NDIMS == 3 */ ! update passive scalars ! if (ns > 0) then ! reset passive scalar increments ! pdata%du(isl:isu,:,:,:) = 0.0d+00 #if NDIMS == 3 do k = nbl - 1, neu kp1 = k + 1 #endif /* NDIMS == 3 */ do j = nbl - 1, neu jp1 = j + 1 do i = nbl - 1, neu ip1 = i + 1 ! X-face ! if (f(idn,i,j,k,1) >= 0.0d+00) then do p = isl, isu df = dhi(1) * f(idn,i,j,k,1) * pdata%q(p,i ,j,k) pdata%du(p,i ,j,k) = pdata%du(p,i ,j,k) - df pdata%du(p,ip1,j,k) = pdata%du(p,ip1,j,k) + df end do else do p = isl, isu df = dhi(1) * f(idn,i,j,k,1) * pdata%q(p,ip1,j,k) pdata%du(p,i ,j,k) = pdata%du(p,i ,j,k) - df pdata%du(p,ip1,j,k) = pdata%du(p,ip1,j,k) + df end do end if ! Y-face ! if (f(idn,i,j,k,2) >= 0.0d+00) then do p = isl, isu df = dhi(2) * f(idn,i,j,k,2) * pdata%q(p,i,j ,k) pdata%du(p,i,j ,k) = pdata%du(p,i,j ,k) - df pdata%du(p,i,jp1,k) = pdata%du(p,i,jp1,k) + df end do else do p = isl, isu df = dhi(2) * f(idn,i,j,k,2) * pdata%q(p,i,jp1,k) pdata%du(p,i,j ,k) = pdata%du(p,i,j ,k) - df pdata%du(p,i,jp1,k) = pdata%du(p,i,jp1,k) + df end do end if #if NDIMS == 3 ! Z-face ! if (f(idn,i,j,k,3) >= 0.0d+00) then do p = isl, isu df = dhi(3) * f(idn,i,j,k,3) * pdata%q(p,i,j,k ) pdata%du(p,i,j,k ) = pdata%du(p,i,j,k ) - df pdata%du(p,i,j,kp1) = pdata%du(p,i,j,kp1) + df end do else do p = isl, isu df = dhi(3) * f(idn,i,j,k,3) * pdata%q(p,i,j,kp1) pdata%du(p,i,j,k ) = pdata%du(p,i,j,k ) - df pdata%du(p,i,j,kp1) = pdata%du(p,i,j,kp1) + df end do end if #endif /* NDIMS == 3 */ end do ! i = nbl, neu end do ! j = nbl, neu #if NDIMS == 3 end do ! k = nbl, neu #endif /* NDIMS == 3 */ end if #ifdef PROFILE ! stop accounting time for the increment update ! call stop_timer(iui) #endif /* PROFILE */ !------------------------------------------------------------------------------- ! end subroutine update_increment ! !=============================================================================== ! ! subroutine UPDATE_VARIABLES: ! --------------------------- ! ! Subroutine iterates over all data blocks and converts the conservative ! variables to their primitive representation. ! ! Arguments: ! ! tm - time at the moment of update; ! dtm - time step since the last update; ! status - status flag indicating if the update was successful; ! !=============================================================================== ! subroutine update_variables(tm, dtm, status) use blocks , only : block_meta, list_meta use blocks , only : block_data, list_data use blocks , only : set_neighbors_update use boundaries, only : boundary_variables use equations , only : update_primitive_variables use equations , only : fix_unphysical_cells, correct_unphysical_states #ifdef MPI use mpitools , only : reduce_maximum #endif /* MPI */ use shapes , only : update_shapes implicit none real(kind=8), intent(in) :: tm, dtm integer , intent(out) :: status type(block_meta), pointer :: pmeta type(block_data), pointer :: pdata !------------------------------------------------------------------------------- ! #ifdef PROFILE call start_timer(imv) #endif /* PROFILE */ status = 0 pdata => list_data do while (associated(pdata)) pmeta => pdata%meta if (pmeta%update) then call update_primitive_variables(pdata%u, pdata%q, status) if (status /= 0) go to 100 end if pdata => pdata%next end do 100 continue #ifdef MPI call reduce_maximum(status) #endif /* MPI */ if (status /= 0) go to 200 call boundary_variables(tm, dtm) pdata => list_data do while (associated(pdata)) pmeta => pdata%meta if (pmeta%update) call update_shapes(pdata, tm, dtm) pdata => pdata%next end do if (fix_unphysical_cells) then ! if an unphysical cell appeared in a block while updating its primitive ! variables it could be propagated to its neighbors through boundary update; ! mark all neighbors of such a block to be verified and corrected for ! unphysical cells too ! pmeta => list_meta do while (associated(pmeta)) if (pmeta%update) call set_neighbors_update(pmeta) pmeta => pmeta%next end do ! verify and correct, if necessary, unphysical cells in recently updated blocks ! pdata => list_data do while (associated(pdata)) pmeta => pdata%meta if (pmeta%update) & call correct_unphysical_states(step, pmeta%id, pdata%q, pdata%u) pdata => pdata%next end do end if 200 continue #ifdef PROFILE call stop_timer(imv) #endif /* PROFILE */ !------------------------------------------------------------------------------- ! end subroutine update_variables ! !=============================================================================== ! ! subroutine UPDATE_ERRORS_L2: ! --------------------------- ! ! Subroutine updates the errors with respect to the absolute and relative ! tolerances using the L2 norm. The errors are calculated per variable. ! ! Arguments: ! ! n - the index of the higher order solution ! m - the index of the lower order solution ! !=============================================================================== ! subroutine update_errors_l2(n, m) use blocks , only : block_data, list_data, get_nleafs use coordinates, only : ncells, nb, ne use equations , only : nf, errors #ifdef MPI use mpitools , only : reduce_sum #endif /* MPI */ implicit none integer, intent(in) :: n, m integer :: l real(kind=8) :: fnorm type(block_data), pointer :: pdata !------------------------------------------------------------------------------- ! errors(:) = 0.0d+00 fnorm = 1.0d+00 / (get_nleafs() * ncells**NDIMS) pdata => list_data do while (associated(pdata)) do l = 1, nf #if NDIMS == 3 errors(l) = errors(l) + & sum((abs(pdata%uu(l,nb:ne,nb:ne,nb:ne,n) & - pdata%uu(l,nb:ne,nb:ne,nb:ne,m)) / & (atol + rtol * & max(abs(pdata%uu(l,nb:ne,nb:ne,nb:ne,n)), & abs(pdata%uu(l,nb:ne,nb:ne,nb:ne,m)))))**2) #else /* NDIMS == 3 */ errors(l) = errors(l) + & sum((abs(pdata%uu(l,nb:ne,nb:ne, : ,n) & - pdata%uu(l,nb:ne,nb:ne, : ,m)) / & (atol + rtol * & max(abs(pdata%uu(l,nb:ne,nb:ne, : ,n)), & abs(pdata%uu(l,nb:ne,nb:ne, : ,m)))))**2) #endif /* NDIMS == 3 */ end do pdata => pdata%next end do #ifdef MPI call reduce_sum(errors) #endif /* MPI */ errors(:) = sqrt(fnorm * errors(:)) !------------------------------------------------------------------------------- ! end subroutine update_errors_l2 ! !=============================================================================== ! ! subroutine UPDATE_ERRORS_MAX: ! ---------------------------- ! ! Subroutine updates the errors with respect to the absolute and relative ! tolerances using the maximum norm. The errors are calculated per variable. ! ! Arguments: ! ! n - the index of the higher order solution ! m - the index of the lower order solution ! !=============================================================================== ! subroutine update_errors_max(n, m) use blocks , only : block_data, list_data use coordinates, only : nb, ne use equations , only : nf, errors #ifdef MPI use mpitools , only : reduce_maximum #endif /* MPI */ implicit none integer, intent(in) :: n, m integer :: l type(block_data), pointer :: pdata !------------------------------------------------------------------------------- ! errors(:) = 0.0d+00 pdata => list_data do while (associated(pdata)) do l = 1, nf #if NDIMS == 3 errors(l) = max(errors(l), & maxval(abs(pdata%uu(l,nb:ne,nb:ne,nb:ne,n) & - pdata%uu(l,nb:ne,nb:ne,nb:ne,m)) / & (atol + rtol * & max(abs(pdata%uu(l,nb:ne,nb:ne,nb:ne,n)), & abs(pdata%uu(l,nb:ne,nb:ne,nb:ne,m)))))) #else /* NDIMS == 3 */ errors(l) = max(errors(l), & maxval(abs(pdata%uu(l,nb:ne,nb:ne, : ,n) & - pdata%uu(l,nb:ne,nb:ne, : ,m)) / & (atol + rtol * & max(abs(pdata%uu(l,nb:ne,nb:ne, : ,n)), & abs(pdata%uu(l,nb:ne,nb:ne, : ,m)))))) #endif /* NDIMS == 3 */ end do pdata => pdata%next end do #ifdef MPI call reduce_maximum(errors) #endif /* MPI */ !------------------------------------------------------------------------------- ! end subroutine update_errors_max #ifdef DEBUG ! !=============================================================================== ! ! subroutine CHECK_VARIABLES: ! -------------------------- ! ! Subroutine iterates over all data blocks and converts the conservative ! variables to their primitive representation. ! ! !=============================================================================== ! subroutine check_variables() ! include external procedures ! use coordinates , only : nn => bcells use equations , only : nv, pvars, cvars use ieee_arithmetic, only : ieee_is_nan ! include external variables ! use blocks , only : block_meta, list_meta use blocks , only : block_data, list_data ! local variables are not implicit by default ! implicit none ! local variables ! integer :: i, j, k = 1, p ! local pointers ! type(block_meta), pointer :: pmeta type(block_data), pointer :: pdata ! !------------------------------------------------------------------------------- ! ! associate the pointer with the first block on the data block list ! pdata => list_data ! iterate over all data blocks ! do while (associated(pdata)) ! associate pmeta with the corresponding meta block ! pmeta => pdata%meta ! check if there are NaNs in primitive variables ! #if NDIMS == 3 do k = 1, nn #endif /* NDIMS == 3 */ do j = 1, nn do i = 1, nn do p = 1, nv if (ieee_is_nan(pdata%u(p,i,j,k))) then print *, 'U NaN:', cvars(p), pdata%meta%id, i, j, k end if if (ieee_is_nan(pdata%q(p,i,j,k))) then print *, 'Q NaN:', pvars(p), pdata%meta%id, i, j, k end if end do ! p = 1, nv end do ! i = 1, nn end do ! j = 1, nn #if NDIMS == 3 end do ! k = 1, nn #endif /* NDIMS == 3 */ ! assign pointer to the next block ! pdata => pdata%next end do !------------------------------------------------------------------------------- ! end subroutine check_variables #endif /* DEBUG */ !=============================================================================== ! end module evolution