gkowal/amun-code
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
amun-code/sources/evolution.F90
Grzegorz Kowal df9b301734 EVOLUTION: Restore the initial state in initialize_time_step(). 
The estimation of the initial time step requires the advance of the
system by one step. This, however, alters the primitive variables of
data blocks. Restore the primitive variables to the initial state, once
the initial time step has been estimated, in order to solve this
inconsistency.

Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
2021-12-15 09:11:49 -03:00

4313 lines
136 KiB
Fortran
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

!!******************************************************************************
!!
!! 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 <grzegorz@amuncode.org>
!!
!! 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 <http://www.gnu.org/licenses/>.
!!
!!******************************************************************************
!!
!! 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
implicit none
! 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)
use coordinates, only : toplev, adx, ady
#if NDIMS == 3
use coordinates, only : adz
#endif /* NDIMS == 3 */
use forcing , only : initialize_forcing
use gravity , only : initialize_gravity
use helpers , only : print_message
use mpitools , only : check_status
use parameters , only : get_parameter
use schemes , only : initialize_schemes
use shapes , only : initialize_shapes
use sources , only : initialize_sources
implicit none
logical, intent(in) :: verbose
integer, intent(out) :: status
character(len=255) :: integration = "rk2", error_norm = "l2"
integer :: n
character(len=*), parameter :: loc = 'EVOLUTION::initialize_evolution()'
!-------------------------------------------------------------------------------
!
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")
evolve => evolve_euler
order = 1
registers = 2
name_int = "1st order Euler"
case ("ssprk(2,2)", "SSPRK(2,2)", "rk2", "RK2")
evolve => evolve_ssprk2
order = 2
registers = 2
name_int = "2nd order SSPRK(2,2)"
case ("ssprk(3,3)", "SSPRK(3,3)", "rk3", "RK3")
evolve => evolve_ssprk3
order = 3
registers = 2
name_int = "3ʳᵈ order SSPRK(3,3)"
case ("ssprk(4,3)", "SSPRK(4,3)", "rk3.4", "RK3.4")
evolve => evolve_ssprk34
order = 3
registers = 2
cfl = 2.0d+00 * cfl
name_int = "3ʳᵈ order 4-step SSPRK(4,3)"
case ("ssprk(5,3)", "SSPRK(5,3)", "rk3.5", "RK3.5")
evolve => evolve_ssprk35
order = 3
registers = 3
cfl = 2.65062919143939d+00 * cfl
name_int = "3ʳᵈ order 5-step SSPRK(5,3)"
case ("ssprk(5,4)", "SSPRK(5,4)", "rk4.5", "RK4.5")
evolve => evolve_ssprk45
order = 4
registers = 3
cfl = 1.50818d+00 * cfl
name_int = "4th order 5-step SSPRK(5,4)"
case ("ssprk(10,4)", "SSPRK(10,4)", "rk4.10", "RK4.10")
evolve => evolve_ssprk410
order = 4
registers = 3
cfl = 6.0d+00 * cfl
name_int = "4th order 10-step SSPRK(10,4)"
case ("ssprk(m,2)", "SSPRK(m,2)", "ssprk2(1)m", "SSPRK2(1)M")
error_control = .true.
evolve => evolve_ssprk21m
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", "ssprk(m,3)", "SSPRK3(2)M", "SSP3(2)M", "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 */
! initialize other modules
!
call initialize_schemes(verbose, status)
if (status /=0) &
call print_message(loc, "Could not initialize module SCHEMES!")
if (check_status(status /= 0)) return
call initialize_sources(verbose, status)
if (status /=0) &
call print_message(loc, "Could not initialize module SOURCES!")
if (check_status(status /= 0)) return
call initialize_forcing(verbose, status)
if (status /=0) &
call print_message(loc, "Could not initialize module FORCING!")
if (check_status(status /= 0)) return
call initialize_gravity(verbose, status)
if (status /=0) &
call print_message(loc, "Could not initialize module GRAVITY!")
if (check_status(status /= 0)) return
call initialize_shapes(status)
if (status /=0) &
call print_message(loc, "Could not initialize module SHAPES!")
if (check_status(status /= 0)) return
!-------------------------------------------------------------------------------
!
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)
use forcing, only : finalize_forcing
use gravity, only : finalize_gravity
use helpers, only : print_section, print_parameter, print_message
use schemes, only : finalize_schemes
use shapes , only : finalize_shapes
use sources, only : finalize_sources
implicit none
logical, intent(in) :: verbose
integer, intent(out) :: status
character(len=*), parameter :: loc = 'EVOLUTION::finalize_evolution()'
!-------------------------------------------------------------------------------
!
status = 0
call finalize_shapes(status)
if (status /=0) &
call print_message(loc, "Could not finalize module SHAPES!")
call finalize_gravity(status)
if (status /=0) &
call print_message(loc, "Could not finalize module GRAVITY!")
call finalize_forcing(status)
if (status /=0) &
call print_message(loc, "Could not finalize module FORCING!")
call finalize_sources(status)
if (status /=0) &
call print_message(loc, "Could not finalize module SOURCES!")
call finalize_schemes(status)
if (status /=0) &
call print_message(loc, "Could not finalize module SCHEMES!")
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)
! 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
!-------------------------------------------------------------------------------
!
end subroutine finalize_evolution
!
!===============================================================================
!
! subroutine PRINT_EVOLUTION:
! --------------------------
!
! Subroutine prints module parameters.
!
! Arguments:
!
! verbose - a logical flag turning the information printing;
!
!===============================================================================
!
subroutine print_evolution(verbose)
use equations, only : magnetized
use forcing , only : print_forcing
use helpers , only : print_section, print_parameter
use schemes , only : print_schemes
use shapes , only : print_shapes
use sources , only : print_sources
implicit none
logical, intent(in) :: verbose
!-------------------------------------------------------------------------------
!
call print_sources(verbose)
call print_forcing(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
call print_shapes(verbose)
call print_schemes(verbose)
!-------------------------------------------------------------------------------
!
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
!
!-------------------------------------------------------------------------------
!
! 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
!-------------------------------------------------------------------------------
!
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, data_blocks, get_nleafs, get_dblocks
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 helpers , only : print_message
#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 :: l, m, n, status
real(kind=8) :: cm, dx_min, fnorm, h0, h1
real(kind=8), dimension(3) :: d
real(kind=8), dimension(:,:), allocatable :: s
#if NDIMS == 3
real(kind=8), dimension(nf,nc,nc,nc) :: sc, df
#else /* NDIMS == 3 */
real(kind=8), dimension(nf,nc,nc, 1) :: sc, df
#endif /* NDIMS == 3 */
real(kind=8), parameter :: eps = tiny(cmax)
character(len=*), parameter :: loc = 'EVOLUTION:initialize_time_step()'
!
!-------------------------------------------------------------------------------
!
cmax = eps
m = 1
n = get_dblocks()
!$omp parallel do default(shared) private(pdata,cm)
do l = 1, n
pdata => data_blocks(l)%ptr
cm = maxspeed(pdata%q)
!$omp atomic update
cmax = max(cmax, cm)
!$omp atomic update
m = max(m, pdata%meta%level)
end do
!$omp end parallel do
#ifdef MPI
call reduce_maximum(cmax)
call reduce_maximum(m)
#endif /* MPI */
cmax2 = cmax * cmax
#if NDIMS == 2
dx_min = min(adx(m), ady(m))
#endif /* NDIMS == 2 */
#if NDIMS == 3
dx_min = min(adx(m), ady(m), adz(m))
#endif /* NDIMS == 3 */
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
allocate(s(n,3))
s(:,:) = 0.0d+00
fnorm = 1.0d+00 / (get_nleafs() * nc**NDIMS)
d(:) = 0.0d+00
!$omp parallel do default(shared) private(pdata,sc)
do l = 1, n
pdata => data_blocks(l)%ptr
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))
s(l,1) = sum((pdata%uu(:,nb:ne,nb:ne,nb:ne,1) / sc(:,:,:,:))**2)
s(l,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))
s(l,1) = sum((pdata%uu(:,nb:ne,nb:ne, : ,1) / sc(:,:,:,:))**2)
s(l,2) = sum((pdata%du(:,nb:ne,nb:ne, : ) / sc(:,:,:,:))**2)
#endif /* NDIMS == 3 */
end do
!$omp end parallel do
d(1) = sum(s(:,1))
d(2) = sum(s(:,2))
#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
m = 10
flag = .true.
do while(flag)
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) + h0 * pdata%du(:,:,:,:)
pdata%u => pdata%uu(:,:,:,:,2)
end do
!$omp end parallel do
call update_variables(time + h0, h0, status)
if (status /= 0) then
h0 = 2.5d-01 * h0
flag = m > 0
m = m - 1
else
flag = .false.
end if
end do
if (status /= 0) then
call print_message(loc, "Could not estimate the initial time step " // &
"due to the error. Setting it to the CFL time step.")
dte = dt
else
!$omp parallel do default(shared) private(pdata,df,sc)
do l = 1, n
pdata => data_blocks(l)%ptr
#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 */
s(l,3) = sum((df(:,:,:,:) / sc(:,:,:,:))**2)
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1)
pdata%u => pdata%uu(:,:,:,:,1)
end do
!$omp end parallel do
d(3) = sum(s(:,3))
deallocate(s)
call update_variables(time, 0.0d+00, status)
#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
end if
!-------------------------------------------------------------------------------
!
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, data_blocks, get_dblocks
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 :: l, n, m
real(kind=8) :: cm, dx_min
real(kind=8), parameter :: eps = tiny(cmax)
!-------------------------------------------------------------------------------
!
cmax = eps
m = 1
n = get_dblocks()
!$omp parallel do default(shared) private(pdata,cm)
do l = 1, n
pdata => data_blocks(l)%ptr
cm = maxspeed(pdata%q)
!$omp atomic update
cmax = max(cmax, cm)
!$omp atomic update
m = max(m, pdata%meta%level)
end do
!$omp end parallel do
#ifdef MPI
call reduce_maximum(cmax)
call reduce_maximum(m)
#endif /* MPI */
cmax2 = cmax * cmax
#if NDIMS == 2
dx_min = min(adx(m), ady(m))
#endif /* NDIMS == 2 */
#if NDIMS == 3
dx_min = min(adx(m), ady(m), adz(m))
#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)
!-------------------------------------------------------------------------------
!
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, data_blocks, get_dblocks
use boundaries, only : boundary_fluxes
use equations , only : ibp, cmax
use sources , only : update_sources
implicit none
type(block_data), pointer :: pdata
integer :: l, n, status
real(kind=8) :: t
!-------------------------------------------------------------------------------
!
status = 1
n = get_dblocks()
do while(status /= 0)
t = time + dt
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1)
pdata%u => pdata%uu(:,:,:,:,2)
end do
!$omp end parallel do
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
call update_increment(pdata)
call update_sources(pdata, t, dt, pdata%du(:,:,:,:))
end do
!$omp end parallel do
call boundary_fluxes()
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + dt * pdata%du(:,:,:,:)
end do
!$omp end parallel do
if (ibp > 0) then
adecay(:) = exp(aglm(:) * cmax * dt)
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%u(ibp,:,:,:) = adecay(pdata%meta%level) * pdata%u(ibp,:,:,:)
end do
!$omp end parallel do
end if
call update_variables(t, dt, status)
if (status /= 0) dt = 2.5d-01 * dt
end do
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2)
pdata%u => pdata%uu(:,:,:,:,1)
end do
!$omp end parallel do
!-------------------------------------------------------------------------------
!
end subroutine evolve_euler
!
!===============================================================================
!
! subroutine EVOLVE_SSPRK2:
! ------------------------
!
! Subroutine advances the solution by one time step using the 2nd order
! Strong Stability Preserving 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_ssprk2()
use blocks , only : block_data, data_blocks, get_dblocks
use boundaries, only : boundary_fluxes
use equations , only : ibp, cmax
use sources , only : update_sources
implicit none
type(block_data), pointer :: pdata
integer :: l, n, status
real(kind=8) :: t, ds
!-------------------------------------------------------------------------------
!
status = 1
n = get_dblocks()
do while(status /= 0)
ds = 5.0d-01 * dt
!= 1st step: U(1) = U(n) + dt * F[U(n)]
!
t = time + dt
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
call update_increment(pdata)
call update_sources(pdata, t, 0.0d+00, pdata%du(:,:,:,:))
end do
!$omp end parallel do
call boundary_fluxes()
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) + dt * pdata%du(:,:,:,:)
pdata%u => pdata%uu(:,:,:,:,2)
end do
!$omp end parallel 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)]
!
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
call update_increment(pdata)
call update_sources(pdata, t, ds, pdata%du(:,:,:,:))
end do
!$omp end parallel do
call boundary_fluxes()
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,2) = 5.0d-01 * (pdata%uu(:,:,:,:,1) &
+ pdata%uu(:,:,:,:,2)) &
+ ds * pdata%du(:,:,:,:)
end do
!$omp end parallel do
if (ibp > 0) then
adecay(:) = exp(aglm(:) * cmax * dt)
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%u(ibp,:,:,:) = adecay(pdata%meta%level) * pdata%u(ibp,:,:,:)
end do
!$omp end parallel do
end if
call update_variables(t, ds, status)
100 continue
if (status /= 0) dt = 2.5d-01 * dt
end do
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2)
pdata%u => pdata%uu(:,:,:,:,1)
end do
!$omp end parallel do
!-------------------------------------------------------------------------------
!
end subroutine evolve_ssprk2
!
!===============================================================================
!
! subroutine EVOLVE_SSPRK3:
! ------------------------
!
! Subroutine advances the solution by one time step using the 3rd order
! Strong Stability Preserving 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_ssprk3()
use blocks , only : block_data, data_blocks, get_dblocks
use boundaries, only : boundary_fluxes
use equations , only : ibp, cmax
use sources , only : update_sources
implicit none
type(block_data), pointer :: pdata
integer :: l, n, 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
!-------------------------------------------------------------------------------
!
status = 1
n = get_dblocks()
do while(status /= 0)
!= 1st substep: U(1) = U(n) + dt F[U(n)]
!
t = time + dt
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
call update_increment(pdata)
call update_sources(pdata, t, 0.0d+00, pdata%du(:,:,:,:))
end do
!$omp end parallel do
call boundary_fluxes()
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) + dt * pdata%du(:,:,:,:)
pdata%u => pdata%uu(:,:,:,:,2)
end do
!$omp end parallel 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
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
call update_increment(pdata)
call update_sources(pdata, t, ds, pdata%du(:,:,:,:))
end do
!$omp end parallel do
call boundary_fluxes()
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,2) = f21 * pdata%uu(:,:,:,:,1) &
+ f22 * pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:)
end do
!$omp end parallel 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
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
call update_increment(pdata)
call update_sources(pdata, t, ds, pdata%du(:,:,:,:))
end do
!$omp end parallel do
call boundary_fluxes()
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,2) = f31 * pdata%uu(:,:,:,:,1) &
+ f32 * pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:)
pdata%u => pdata%uu(:,:,:,:,1)
end do
!$omp end parallel do
if (ibp > 0) then
adecay(:) = exp(aglm(:) * cmax * dt)
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%u(ibp,:,:,:) = adecay(pdata%meta%level) * pdata%u(ibp,:,:,:)
end do
!$omp end parallel do
end if
call update_variables(t, dt, status)
100 continue
if (status /= 0) dt = 2.5d-01 * dt
end do
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2)
pdata%u => pdata%uu(:,:,:,:,1)
end do
!$omp end parallel do
!-------------------------------------------------------------------------------
!
end subroutine evolve_ssprk3
!
!===============================================================================
!
! 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, data_blocks, get_dblocks
use boundaries, only : boundary_fluxes
use equations , only : ibp, cmax
use sources , only : update_sources
implicit none
type(block_data), pointer :: pdata
integer :: l, n, status
real(kind=8) :: t, ds
real(kind=8), parameter :: b1 = 2.0d+00 / 3.0d+00, b2 = 1.0d+00 / 3.0d+00, &
b3 = 1.0d+00 / 6.0d+00
!-------------------------------------------------------------------------------
!
status = 1
n = get_dblocks()
do while(status /= 0)
!= 1st step: U(1) = U(n) + 1/2 dt F[U(n)]
!
ds = 5.0d-01 * dt
t = time + ds
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
call update_increment(pdata)
call update_sources(pdata, t, 0.0d+00, pdata%du(:,:,:,:))
end do
!$omp end parallel do
call boundary_fluxes()
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) + ds * pdata%du(:,:,:,:)
pdata%u => pdata%uu(:,:,:,:,2)
end do
!$omp end parallel 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
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
call update_increment(pdata)
call update_sources(pdata, t, ds, pdata%du(:,:,:,:))
end do
!$omp end parallel do
call boundary_fluxes()
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:)
end do
!$omp end parallel 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/6 dt F[U(2)])
!
t = time + ds
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
call update_increment(pdata)
call update_sources(pdata, t, ds, pdata%du(:,:,:,:))
end do
!$omp end parallel do
call boundary_fluxes()
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,2) = b1 * pdata%uu(:,:,:,:,1) &
+ b2 * pdata%uu(:,:,:,:,2) &
+ b3 * dt * pdata%du(:,:,:,:)
end do
!$omp end parallel 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
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
call update_increment(pdata)
call update_sources(pdata, t, ds, pdata%du(:,:,:,:))
end do
!$omp end parallel do
call boundary_fluxes()
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:)
end do
!$omp end parallel do
if (ibp > 0) then
adecay(:) = exp(aglm(:) * cmax * dt)
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%u(ibp,:,:,:) = adecay(pdata%meta%level) * pdata%u(ibp,:,:,:)
end do
!$omp end parallel do
end if
call update_variables(t, dt, status)
100 continue
if (status /= 0) dt = 2.5d-01 * dt
end do
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2)
pdata%u => pdata%uu(:,:,:,:,1)
end do
!$omp end parallel do
!-------------------------------------------------------------------------------
!
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, data_blocks, get_dblocks
use boundaries, only : boundary_fluxes
use equations , only : ibp, cmax
use sources , only : update_sources
implicit none
type(block_data), pointer :: pdata
integer :: l, n, 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
!
!-------------------------------------------------------------------------------
!
status = 1
n = get_dblocks()
do while(status /= 0)
!= 1st step: U(1) = U(n) + b1 dt F[U(n)]
!
ds = b1 * dt
t = time + ds
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
call update_increment(pdata)
call update_sources(pdata, t, 0.0d+00, pdata%du(:,:,:,:))
end do
!$omp end parallel do
call boundary_fluxes()
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) + ds * pdata%du(:,:,:,:)
pdata%u => pdata%uu(:,:,:,:,2)
end do
!$omp end parallel 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
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
call update_increment(pdata)
call update_sources(pdata, t, ds, pdata%du(:,:,:,:))
end do
!$omp end parallel do
call boundary_fluxes()
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:)
end do
!$omp end parallel 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
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
call update_increment(pdata)
call update_sources(pdata, t, ds, pdata%du(:,:,:,:))
end do
!$omp end parallel do
call boundary_fluxes()
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,2) = a31 * pdata%uu(:,:,:,:,1) &
+ a33 * pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:)
end do
!$omp end parallel 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
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
call update_increment(pdata)
call update_sources(pdata, t, ds, pdata%du(:,:,:,:))
end do
!$omp end parallel do
call boundary_fluxes()
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,3) = a41 * pdata%uu(:,:,:,:,1) &
+ a44 * pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:)
pdata%u => pdata%uu(:,:,:,:,3)
end do
!$omp end parallel 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
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
call update_increment(pdata)
call update_sources(pdata, t, ds, pdata%du(:,:,:,:))
end do
!$omp end parallel do
call boundary_fluxes()
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,2) = a53 * pdata%uu(:,:,:,:,2) &
+ a55 * pdata%uu(:,:,:,:,3) + ds * pdata%du(:,:,:,:)
pdata%u => pdata%uu(:,:,:,:,2)
end do
!$omp end parallel do
if (ibp > 0) then
adecay(:) = exp(aglm(:) * cmax * dt)
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%u(ibp,:,:,:) = adecay(pdata%meta%level) * pdata%u(ibp,:,:,:)
end do
!$omp end parallel do
end if
call update_variables(t, dt, status)
100 continue
if (status /= 0) dt = 2.5d-01 * dt
end do
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2)
pdata%u => pdata%uu(:,:,:,:,1)
end do
!$omp end parallel do
!-------------------------------------------------------------------------------
!
end subroutine evolve_ssprk35
!
!===============================================================================
!
! subroutine EVOLVE_SSPRK45:
! -------------------------
!
! Subroutine advances the solution by one time step using the 4th 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_ssprk45()
use blocks , only : block_data, data_blocks, get_dblocks
use boundaries, only : boundary_fluxes
use equations , only : ibp, cmax
use sources , only : update_sources
implicit none
type(block_data), pointer :: pdata
integer :: l, n, status
real(kind=8) :: t, ds
real(kind=8), parameter :: b11 = 3.91752226571890d-01, &
b22 = 3.68410593050371d-01, &
b33 = 2.51891774271694d-01, &
b44 = 5.44974750228521d-01, &
b54 = 6.36924686662900d-02, &
b55 = 2.26007483236906d-01
real(kind=8), parameter :: a21 = 4.44370493651235d-01, &
a22 = 5.55629506348765d-01
real(kind=8), parameter :: a31 = 6.20101851488403d-01, &
a33 = 3.79898148511597d-01
real(kind=8), parameter :: a41 = 1.78079954393132d-01, &
a44 = 8.21920045606868d-01
real(kind=8), parameter :: a53 = 5.17231671970585d-01, &
a54 = 9.60597105261470d-02, &
a55 = 3.86708617503269d-01
real(kind=8), parameter :: c2 = 3.91752226571890d-01, &
c3 = 5.86079689311541d-01, &
c4 = 4.74542363121399d-01, &
c5 = 9.35010630967652d-01
!-------------------------------------------------------------------------------
!
status = 1
n = get_dblocks()
do while(status /= 0)
!= 1st step: U(1) = U(n) + b11 dt F[U(n)]
!
ds = b11 * dt
t = time
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
call update_increment(pdata)
call update_sources(pdata, t, 0.0d+00, pdata%du(:,:,:,:))
end do
!$omp end parallel do
call boundary_fluxes()
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) + ds * pdata%du(:,:,:,:)
pdata%u => pdata%uu(:,:,:,:,2)
end do
!$omp end parallel do
call update_variables(t, ds, status)
if (status /= 0) go to 100
!= 2nd step: U(2) = a21 U(n) + a22 U(1) + b22 dt F[U(1)]
!
ds = b22 * dt
t = time + c2 * dt
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
call update_increment(pdata)
call update_sources(pdata, t, ds, pdata%du(:,:,:,:))
end do
!$omp end parallel do
call boundary_fluxes()
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,2) = a21 * pdata%uu(:,:,:,:,1) &
+ a22 * pdata%uu(:,:,:,:,2) &
+ ds * pdata%du(:,:,:,:)
end do
!$omp end parallel do
call update_variables(t, ds, status)
if (status /= 0) go to 100
!= 3rd step: U(3) = a31 U(n) + a33 U(2) + b33 dt F[U(2)]
!
ds = b33 * dt
t = time + c3 * dt
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
call update_increment(pdata)
call update_sources(pdata, t, ds, pdata%du(:,:,:,:))
end do
!$omp end parallel do
call boundary_fluxes()
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,3) = a53 * pdata%uu(:,:,:,:,2)
pdata%uu(:,:,:,:,2) = a31 * pdata%uu(:,:,:,:,1) &
+ a33 * pdata%uu(:,:,:,:,2) &
+ ds * pdata%du(:,:,:,:)
end do
!$omp end parallel do
call update_variables(t, ds, status)
if (status /= 0) go to 100
!= 4th step: U(4) = a41 U(n) + a44 U(3) + b44 dt F[U(3)]
!
ds = b44 * dt
t = time + c4 * dt
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
call update_increment(pdata)
call update_sources(pdata, t, ds, pdata%du(:,:,:,:))
end do
!$omp end parallel do
call boundary_fluxes()
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,3) = pdata%uu(:,:,:,:,3) &
+ a54 * pdata%uu(:,:,:,:,2) &
+ b54 * dt * pdata%du(:,:,:,:)
pdata%uu(:,:,:,:,2) = a41 * pdata%uu(:,:,:,:,1) &
+ a44 * pdata%uu(:,:,:,:,2) &
+ ds * pdata%du(:,:,:,:)
end do
!$omp end parallel do
call update_variables(t, ds, status)
if (status /= 0) go to 100
!= the final step: U(n+1) = a54 U(3) + a55 U(4) + b54 dt F[U(3)] + b55 dt F[U(4)]
!
ds = b55 * dt
t = time + c5 * dt
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
call update_increment(pdata)
call update_sources(pdata, t, ds, pdata%du(:,:,:,:))
end do
!$omp end parallel do
call boundary_fluxes()
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,3) &
+ a55 * pdata%uu(:,:,:,:,2) &
+ ds * pdata%du(:,:,:,:)
pdata%u => pdata%uu(:,:,:,:,2)
end do
!$omp end parallel do
if (ibp > 0) then
adecay(:) = exp(aglm(:) * cmax * dt)
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%u(ibp,:,:,:) = adecay(pdata%meta%level) * pdata%u(ibp,:,:,:)
end do
!$omp end parallel do
end if
call update_variables(t, dt, status)
100 continue
if (status /= 0) dt = 2.5d-01 * dt
end do
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2)
pdata%u => pdata%uu(:,:,:,:,1)
end do
!$omp end parallel do
!-------------------------------------------------------------------------------
!
end subroutine evolve_ssprk45
!
!===============================================================================
!
! subroutine EVOLVE_SSPRK410:
! --------------------------
!
! Subroutine advances the solution by one time step using the 4th 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_ssprk410()
use blocks , only : block_data, data_blocks, get_dblocks
use boundaries, only : boundary_fluxes
use equations , only : ibp, cmax
use sources , only : update_sources
implicit none
type(block_data), pointer :: pdata
integer :: i, l, n, status
real(kind=8) :: t
real(kind=8), dimension(9) :: ds
real(kind=8), dimension(9), parameter :: cc = &
[ 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
!-------------------------------------------------------------------------------
!
ds(:) = cc(:) * dt
status = 1
n = get_dblocks()
do while(status /= 0)
!= 1st step: U(2) = U(1)
!
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1)
pdata%uu(:,:,:,:,3) = pdata%uu(:,:,:,:,1)
pdata%u => pdata%uu(:,:,:,:,2)
end do
!$omp end parallel do
!= 2nd step: U(1) = [1 + dt/6 L] U(1), for i = 1, ..., 5
!
do i = 1, 5
t = time + ds(i)
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
call update_increment(pdata)
call update_sources(pdata, t, 0.0d+00, pdata%du(:,:,:,:))
end do
!$omp end parallel do
call boundary_fluxes()
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + ds(1) * pdata%du(:,:,:,:)
end do
!$omp end parallel 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)
!
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
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)
end do
!$omp end parallel 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)
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
call update_increment(pdata)
call update_sources(pdata, t, ds(1), pdata%du(:,:,:,:))
end do
!$omp end parallel do
call boundary_fluxes()
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + ds(1) * pdata%du(:,:,:,:)
end do
!$omp end parallel 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
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
call update_increment(pdata)
call update_sources(pdata, t, 1.0d-01 * dt, pdata%du(:,:,:,:))
end do
!$omp end parallel do
call boundary_fluxes()
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,3) &
+ 6.0d-01 * pdata%uu(:,:,:,:,2) &
+ 1.0d-01 * dt * pdata%du(:,:,:,:)
end do
!$omp end parallel do
if (ibp > 0) then
adecay(:) = exp(aglm(:) * cmax * dt)
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%u(ibp,:,:,:) = adecay(pdata%meta%level) * pdata%u(ibp,:,:,:)
end do
!$omp end parallel do
end if
call update_variables(t, dt, status)
100 continue
if (status /= 0) dt = 2.5d-01 * dt
end do
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2)
pdata%u => pdata%uu(:,:,:,:,1)
end do
!$omp end parallel do
!-------------------------------------------------------------------------------
!
end subroutine evolve_ssprk410
!
!===============================================================================
!
! EMBEDDED METHODS
!
!===============================================================================
!
! subroutine EVOLVE_SSPRK21M:
! --------------------------
!
! 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_ssprk21m()
use blocks , only : block_data, data_blocks, get_dblocks
use boundaries, only : boundary_fluxes
use equations , only : errors, ibp, cmax
use helpers , only : print_message
use sources , only : update_sources
implicit none
type(block_data), pointer :: pdata
logical :: test, cond
integer :: i, l, n, nrej, status
real(kind=8) :: tm, dtm, fc
character(len=*), parameter :: loc = 'EVOLUTION:evolve_ssprk21m()'
!-------------------------------------------------------------------------------
!
test = .true.
cond = .true.
nrej = 0
n = get_dblocks()
! 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)
!
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1)
pdata%uu(:,:,:,:,3) = pdata%uu(:,:,:,:,1)
pdata%u => pdata%uu(:,:,:,:,2)
end do
!$omp end parallel 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
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
call update_increment(pdata)
call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:))
end do
!$omp end parallel do
call boundary_fluxes()
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) &
+ bt(i) * dt * pdata%du(:,:,:,:)
pdata%uu(:,:,:,:,3) = pdata%uu(:,:,:,:,3) &
+ bh(i) * dt * pdata%du(:,:,:,:)
end do
!$omp end parallel 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
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
call update_increment(pdata)
call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:))
end do
!$omp end parallel do
call boundary_fluxes()
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
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(:,:,:,:)
end do
!$omp end parallel 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)
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(ibp,:,:,:,2) = adecay(pdata%meta%level) &
* pdata%uu(ibp,:,:,:,2)
pdata%uu(ibp,:,:,:,3) = adecay(pdata%meta%level) &
* pdata%uu(ibp,:,:,:,3)
end do
!$omp end parallel 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
cond = fc > fac .or. nrej >= mrej
100 continue
cond = cond .and. status == 0
if (cond) 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
!
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%u => pdata%uu(:,:,:,:,1)
end do
!$omp end parallel do
call update_variables(tm, dtm, status)
if (status /= 0) &
call print_message(loc, "Possible bug: the revert to " // &
"the previous state found it unphysical.")
end if
niterations = niterations + 1
end do
!= final step: U(n+1) = U(1) - update the accepted solution
!
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2)
pdata%u => pdata%uu(:,:,:,:,1)
end do
!$omp end parallel do
!-------------------------------------------------------------------------------
!
end subroutine evolve_ssprk21m
!
!===============================================================================
!
! 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, data_blocks, get_dblocks
use boundaries, only : boundary_fluxes
use equations , only : errors, ibp, cmax
use helpers , only : print_message
use sources , only : update_sources
implicit none
type(block_data), pointer :: pdata
logical :: test, cond
integer :: i, l, n, nrej, 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()'
!-------------------------------------------------------------------------------
!
test = .true.
cond = .true.
nrej = 0
n = get_dblocks()
! 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)
!
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1)
pdata%u => pdata%uu(:,:,:,:,2)
end do
!$omp end parallel 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
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
call update_increment(pdata)
call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:))
end do
!$omp end parallel do
call boundary_fluxes()
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + dh * pdata%du(:,:,:,:)
end do
!$omp end parallel 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)
!
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,3) = onethird * pdata%uu(:,:,:,:,1) &
+ twothird * pdata%uu(:,:,:,:,2)
pdata%uu(:,:,:,:,2) = twothird * pdata%uu(:,:,:,:,1) &
+ onethird * pdata%uu(:,:,:,:,2)
end do
!$omp end parallel 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
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
call update_increment(pdata)
call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:))
end do
!$omp end parallel do
call boundary_fluxes()
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + dh * pdata%du(:,:,:,:)
end do
!$omp end parallel do
call update_variables(tm, dtm, status)
if (status /= 0) go to 100
!= 6th step: U(2) = ½ (U(1) + U(2)) <- 2ⁿᵈ-order approximation
!
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,3) = 5.0d-01 * (pdata%uu(:,:,:,:,3) &
+ pdata%uu(:,:,:,:,2))
end do
!$omp end parallel do
!= 7th step: decay the magnetic divergence potential of both solutions
!
if (ibp > 0) then
adecay(:) = exp(aglm(:) * cmax * dt)
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(ibp,:,:,:,2) = adecay(pdata%meta%level) &
* pdata%uu(ibp,:,:,:,2)
pdata%uu(ibp,:,:,:,3) = adecay(pdata%meta%level) &
* pdata%uu(ibp,:,:,:,3)
end do
!$omp end parallel 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
cond = fc > fac .or. nrej >= mrej
100 continue
cond = cond .and. status == 0
if (cond) 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
!
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%u => pdata%uu(:,:,:,:,1)
end do
!$omp end parallel do
call update_variables(tm, dtm, status)
if (status /= 0) &
call print_message(loc, "Possible bug: the revert to " // &
"the previous state found it unphysical.")
end if
niterations = niterations + 1
end do
!= final step: U(n+1) = U(1) - update the accepted solution
!
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2)
pdata%u => pdata%uu(:,:,:,:,1)
end do
!$omp end parallel do
!-------------------------------------------------------------------------------
!
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, data_blocks, get_dblocks
use boundaries, only : boundary_fluxes
use equations , only : errors, ibp, cmax
use helpers , only : print_message
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, cond
integer :: i, l, nrej, status
real(kind=8) :: tm, dtm, dh, fc
character(len=*), parameter :: loc = 'EVOLUTION:evolve_ssprk32m()'
!-------------------------------------------------------------------------------
!
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.
cond = .true.
nrej = 0
n = get_dblocks()
! 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)
!
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1)
pdata%u => pdata%uu(:,:,:,:,2)
end do
!$omp end parallel 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
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
call update_increment(pdata)
call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:))
end do
!$omp end parallel do
call boundary_fluxes()
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + dh * pdata%du(:,:,:,:)
end do
!$omp end parallel do
call update_variables(tm, dtm, status)
if (status /= 0) go to 100
end do ! i = 1, i1
!= 2nd step: U(2) = U(1)
!
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,3) = pdata%uu(:,:,:,:,2)
end do
!$omp end parallel 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
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
call update_increment(pdata)
call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:))
end do
!$omp end parallel do
call boundary_fluxes()
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + dh * pdata%du(:,:,:,:)
end do
!$omp end parallel 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)
!
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
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)
end do
!$omp end parallel 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
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
call update_increment(pdata)
call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:))
end do
!$omp end parallel do
call boundary_fluxes()
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + dh * pdata%du(:,:,:,:)
end do
!$omp end parallel 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)
!
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,4) = pdata%uu(:,:,:,:,4) + f4 * pdata%uu(:,:,:,:,2)
end do
!$omp end parallel do
!= 8th step: decay the magnetic divergence potential of both solutions
!
if (ibp > 0) then
adecay(:) = exp(aglm(:) * cmax * dt)
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(ibp,:,:,:,2) = adecay(pdata%meta%level) &
* pdata%uu(ibp,:,:,:,2)
pdata%uu(ibp,:,:,:,4) = adecay(pdata%meta%level) &
* pdata%uu(ibp,:,:,:,4)
end do
!$omp end parallel 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
cond = fc > fac .or. nrej >= mrej
100 continue
cond = cond .and. status == 0
if (cond) 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
!
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%u => pdata%uu(:,:,:,:,1)
end do
!$omp end parallel do
call update_variables(tm, dtm, status)
if (status /= 0) &
call print_message(loc, "Possible bug: the revert to " // &
"the previous state found it unphysical.")
end if
niterations = niterations + 1
end do
!= final step: U(n+1) = U(1) - update the accepted solution
!
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2)
pdata%u => pdata%uu(:,:,:,:,1)
end do
!$omp end parallel do
!-------------------------------------------------------------------------------
!
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, data_blocks, get_dblocks
use boundaries, only : boundary_fluxes
use equations , only : errors, ibp, cmax
use helpers , only : print_message
use sources , only : update_sources
implicit none
type(block_data), pointer :: pdata
logical :: test, cond
integer :: i, l, n, nrej, status
real(kind=8) :: tm, dtm, fc
character(len=*), parameter :: loc = 'EVOLUTION:evolve_3sstarplus()'
!-------------------------------------------------------------------------------
!
test = .true.
cond = .true.
nrej = 0
n = get_dblocks()
! 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)
!
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1)
pdata%uu(:,:,:,:,3) = pdata%uu(:,:,:,:,1)
pdata%uu(:,:,:,:,4) = pdata%uu(:,:,:,:,1)
pdata%u => pdata%uu(:,:,:,:,2)
end do
!$omp end parallel do
!= 1st stage
!
if (fsal_update) then
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
call update_increment(pdata)
call update_sources(pdata, time, 0.0d+00, pdata%du(:,:,:,:))
end do
!$omp end parallel do
call boundary_fluxes()
end if
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) &
+ bt(1) * dt * pdata%du(:,:,:,:)
pdata%uu(:,:,:,:,4) = pdata%uu(:,:,:,:,4) &
+ bh(1) * dt * pdata%du(:,:,:,:)
end do
!$omp end parallel 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
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
call update_increment(pdata)
call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:))
end do
!$omp end parallel do
call boundary_fluxes()
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
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(:,:,:,:)
end do
!$omp end parallel 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)
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(ibp,:,:,:,2) = adecay(pdata%meta%level) &
* pdata%uu(ibp,:,:,:,2)
pdata%uu(ibp,:,:,:,4) = adecay(pdata%meta%level) &
* pdata%uu(ibp,:,:,:,4)
end do
!$omp end parallel do
end if
!= extra step: FSAL
!
if (fsal) then
tm = time + dt
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
call update_increment(pdata)
call update_sources(pdata, tm, 0.0d+00, pdata%du(:,:,:,:))
end do
!$omp end parallel do
call boundary_fluxes()
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,4) = pdata%uu(:,:,:,:,4) &
+ bh(stages+1) * dt * pdata%du(:,:,:,:)
end do
!$omp end parallel 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
cond = fc > fac .or. nrej >= mrej
100 continue
cond = cond .and. status == 0
if (cond) 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
!
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%u => pdata%uu(:,:,:,:,1)
end do
!$omp end parallel do
call update_variables(tm, dtm, status)
if (status /= 0) &
call print_message(loc, "Possible bug: the revert to " // &
"the previous state found it unphysical.")
if (fsal) fsal_update = .false.
end if
niterations = niterations + 1
end do
!= final step: U(n+1) = U(1)
!
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2)
pdata%u => pdata%uu(:,:,:,:,1)
end do
!$omp end parallel do
tm = time + dt
call update_variables(tm, dt, status)
if (status /= 0) &
call print_message(loc, "Possible bug: the new state has been " // &
"accepted, nevertheless the variable " // &
"update found it unphysical.")
!-------------------------------------------------------------------------------
!
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)
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 helpers , only : print_message
use schemes , only : update_flux
use workspace , only : resize_workspace, work, work_in_use
implicit none
type(block_data), pointer, intent(inout) :: pdata
logical, save :: first = .true.
integer :: i , j , k = 1, p
integer :: im1, ip1
integer :: jm1, jp1
#if NDIMS == 3
integer :: km1, kp1
#endif /* NDIMS == 3 */
integer :: status
real(kind=8) :: df
real(kind=8), dimension(NDIMS) :: dh, dhi
integer :: nt = 0
!$ integer :: omp_get_thread_num
real(kind=8), dimension(:,:,:,:,:) , pointer, save :: f
real(kind=8), dimension(:,:,:,:,:,:), pointer, save :: s
!$omp threadprivate(first,nt,f,s)
character(len=*), parameter :: loc = 'EVOLUTION:update_increment()'
!-------------------------------------------------------------------------------
!
!$ nt = omp_get_thread_num()
if (first) then
i = NDIMS * nf * nn**NDIMS
j = 3 * i
call resize_workspace(j, status)
if (status /= 0) then
call print_message(loc, "Could not resize the workspace!")
go to 100
end if
#if NDIMS == 3
f(1:nf,1:nn,1:nn,1:nn,1:3) => work( 1:i,nt)
s(1:nf,1:nn,1:nn,1:nn,1:2,1:3) => work(i+1:j,nt)
#else /* NDIMS == 3 */
f(1:nf,1:nn,1:nn,1: 1,1:2) => work( 1:i,nt)
s(1:nf,1:nn,1:nn,1: 1,1:2,1:2) => work(i+1:j,nt)
#endif /* NDIMS == 3 */
first = .false.
end if
if (work_in_use(nt)) &
call print_message(loc, "Workspace is being used right now! " // &
"Corruptions can occur!")
work_in_use(nt) = .true.
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(:,:,:,:), &
s(:,:,:,:,:,:), f(:,:,:,:,:))
! store the 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
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
work_in_use(nt) = .false.
100 continue
!-------------------------------------------------------------------------------
!
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, data_blocks
use blocks , only : set_neighbors_update, get_dblocks
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
integer :: l, n, s
!-------------------------------------------------------------------------------
!
status = 0
n = get_dblocks()
!$omp parallel do default(shared) private(pdata,s)
do l = 1, n
pdata => data_blocks(l)%ptr
if (pdata%meta%update) then
call update_primitive_variables(pdata%u, pdata%q, s)
!$omp critical
if (s /= 0) status = 1
!$omp end critical
end if
end do
!$omp end parallel do
#ifdef MPI
call reduce_maximum(status)
#endif /* MPI */
if (status /= 0) go to 100
call boundary_variables(tm, dtm, status)
if (status /= 0) go to 100
!$omp parallel do default(shared) private(pdata)
do l = 1, n
pdata => data_blocks(l)%ptr
if (pdata%meta%update) &
call update_shapes(pdata, tm, dtm)
end do
!$omp end parallel 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
100 continue
!-------------------------------------------------------------------------------
!
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:
!
! nh - the index of the higher order solution
! nl - the index of the lower order solution
!
!===============================================================================
!
subroutine update_errors_l2(nh, nl)
use blocks , only : block_data, data_blocks
use blocks , only : get_nleafs, get_dblocks
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) :: nh, nl
integer :: l, n, p
real(kind=8) :: fnorm, err
type(block_data), pointer :: pdata
!-------------------------------------------------------------------------------
!
errors(:) = 0.0d+00
fnorm = 1.0d+00 / (get_nleafs() * ncells**NDIMS)
n = get_dblocks()
!$omp parallel do default(shared) private(pdata,p,err)
do l = 1, n
pdata => data_blocks(l)%ptr
do p = 1, nf
#if NDIMS == 3
err = sum((abs(pdata%uu(p,nb:ne,nb:ne,nb:ne,nh) &
- pdata%uu(p,nb:ne,nb:ne,nb:ne,nl)) / (atol + rtol * &
max(abs(pdata%uu(p,nb:ne,nb:ne,nb:ne,nh)), &
abs(pdata%uu(p,nb:ne,nb:ne,nb:ne,nl)))))**2)
#else /* NDIMS == 3 */
err = sum((abs(pdata%uu(p,nb:ne,nb:ne, : ,nh) &
- pdata%uu(p,nb:ne,nb:ne, : ,nl)) / (atol + rtol * &
max(abs(pdata%uu(p,nb:ne,nb:ne, : ,nh)), &
abs(pdata%uu(p,nb:ne,nb:ne, : ,nl)))))**2)
#endif /* NDIMS == 3 */
!$omp atomic update
errors(p) = errors(p) + err
end do
end do
!$omp end parallel 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:
!
! nh - the index of the higher order solution
! nl - the index of the lower order solution
!
!===============================================================================
!
subroutine update_errors_max(nh, nl)
use blocks , only : block_data, data_blocks, get_dblocks
use coordinates, only : nb, ne
use equations , only : nf, errors
#ifdef MPI
use mpitools , only : reduce_maximum
#endif /* MPI */
implicit none
integer, intent(in) :: nh, nl
integer :: l, n, p
real(kind=8) :: err
type(block_data), pointer :: pdata
!-------------------------------------------------------------------------------
!
errors(:) = 0.0d+00
n = get_dblocks()
!$omp parallel do default(shared) private(pdata,p,err)
do l = 1, n
pdata => data_blocks(l)%ptr
do p = 1, nf
#if NDIMS == 3
err = maxval(abs(pdata%uu(p,nb:ne,nb:ne,nb:ne,nh) &
- pdata%uu(p,nb:ne,nb:ne,nb:ne,nl)) / (atol + rtol * &
max(abs(pdata%uu(p,nb:ne,nb:ne,nb:ne,nh)), &
abs(pdata%uu(p,nb:ne,nb:ne,nb:ne,nl)))))
#else /* NDIMS == 3 */
err = maxval(abs(pdata%uu(p,nb:ne,nb:ne, : ,nh) &
- pdata%uu(p,nb:ne,nb:ne, : ,nl)) / (atol + rtol * &
max(abs(pdata%uu(p,nb:ne,nb:ne, : ,nh)), &
abs(pdata%uu(p,nb:ne,nb:ne, : ,nl)))))
#endif /* NDIMS == 3 */
!$omp atomic update
errors(p) = max(errors(p), err)
end do
end do
!$omp end parallel 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
Reference in New Issue View Git Blame Copy Permalink