gkowal/amun-code
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

EVOLUTION: Add limiting factors and correct exponents in error control.

Browse Source 
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
Grzegorz Kowal 2020-09-15 22:59:20 -03:00
parent 2450ec2a34
commit 36b10a40d7
1 changed files with 41 additions and 9 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

50
sources/evolution.F90
View File

@ -83,6 +83,9 @@ module evolution
! !
real(kind=8) , save :: atol = 1.0d-04 real(kind=8) , save :: atol = 1.0d-04
real(kind=8) , save :: rtol = 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
integer , save :: mrej = 5 integer , save :: mrej = 5
integer , save :: niterations = 0 integer , save :: niterations = 0
integer , save :: nrejections = 0 integer , save :: nrejections = 0
@ -190,11 +193,14 @@ module evolution
call get_parameter("glm_alpha" , glm_alpha ) call get_parameter("glm_alpha" , glm_alpha )
! get the absolute and relative tolerances, the maximum number of passes for ! get the absolute and relative tolerances, the maximum number of passes for
! the adaptive step ! the adaptive step, and limiting factors
! !
call get_parameter("absolute_tolerance", atol) call get_parameter("absolute_tolerance", atol)
call get_parameter("relative_tolerance", rtol) call get_parameter("relative_tolerance", rtol)
call get_parameter("maximum_rejections", mrej) call get_parameter("maximum_rejections", mrej)
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 ! select the integration method, check the correctness of the integration
! parameters and adjust the CFL coefficient if necessary ! parameters and adjust the CFL coefficient if necessary
@ -426,6 +432,9 @@ module evolution
call print_parameter(verbose, "absolute tolerance", atol) call print_parameter(verbose, "absolute tolerance", atol)
call print_parameter(verbose, "relative tolerance", rtol) call print_parameter(verbose, "relative tolerance", rtol)
call print_parameter(verbose, "maximum rejections", mrej) call print_parameter(verbose, "maximum rejections", mrej)
call print_parameter(verbose, "factor" , fac)
call print_parameter(verbose, "minimum factor" , facmin)
call print_parameter(verbose, "maximum factor" , facmax)
end if end if
@ -954,11 +963,14 @@ module evolution
logical :: test logical :: test
integer :: n, l, nrej integer :: n, l, nrej
real(kind=8) :: tm, dtm, ds, umax, utol, emax real(kind=8) :: tm, dtm, ds, umax, utol, emax
real(kind=8) :: fc, fcmn, fcmx
logical , save :: first = .true. logical , save :: first = .true.
real(kind=8), save :: ft, fl, fr, gt, gl real(kind=8), save :: ft, fl, fr, gt, gl
real(kind=8), dimension(:,:,:,:,:), allocatable :: lerr real(kind=8), dimension(:,:,:,:,:), allocatable :: lerr
real(kind=8), parameter :: k1 = -2.9d-01, k2 = 1.05d-01, k3 = -5.0d-02
! !
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
@ -979,6 +991,10 @@ module evolution
end if end if
fc = fac
fcmn = facmin
fcmx = facmax
l = get_dblocks() l = get_dblocks()
#if NDIMS == 3 #if NDIMS == 3
allocate(lerr(l,nf,nc,nc,nc)) allocate(lerr(l,nf,nc,nc,nc))
@ -1105,15 +1121,19 @@ module evolution
errs(2) = errs(1) errs(2) = errs(1)
errs(1) = emax errs(1) = emax
dte = dt * errs(1)**(-0.2d+00) * errs(2)**(0.0505d+00) & dte = dt * min(fcmx, max(fcmn, &
* errs(3)**(-0.025d+00) fc * errs(1)**k1 * errs(2)**k2 * errs(3)**k3))
fcmx = facmax
else else
errs(1) = emax errs(1) = emax
dte = dt * errs(1)**(-0.2d+00) * errs(2)**(0.0505d+00) & dte = dt * min(fcmx, max(fcmn, &
* errs(3)**(-0.025d+00) fc * errs(1)**k1 * errs(2)**k2 * errs(3)**k3))
dt = dte dt = dte
fcmx = fac
nrej = nrej + 1 ! rejection count in the current step nrej = nrej + 1 ! rejection count in the current step
nrejections = nrejections + 1 nrejections = nrejections + 1
end if end if
@ -1944,8 +1964,12 @@ module evolution
logical :: test logical :: test
integer :: i, l, nrej integer :: i, l, nrej
real(kind=8) :: tm, dtm, ds, umax, utol, emax real(kind=8) :: tm, dtm, ds, umax, utol, emax
real(kind=8) :: fc, fcmn, fcmx
real(kind=8), dimension(:,:,:,:,:), allocatable :: lerr real(kind=8), dimension(:,:,:,:,:), allocatable :: lerr
real(kind=8), parameter :: k1 = -5.8d-01 / 3.0d+00, k2 = 7.0d-02, &
k3 = -1.0d-01 / 3.0d+00
! !
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
@ -1973,6 +1997,10 @@ module evolution
end if end if
fc = fac
fcmn = facmin
fcmx = facmax
l = get_dblocks() l = get_dblocks()
#if NDIMS == 3 #if NDIMS == 3
allocate(lerr(l,nf,nc,nc,nc)) allocate(lerr(l,nf,nc,nc,nc))
@ -2179,15 +2207,19 @@ module evolution
errs(2) = errs(1) errs(2) = errs(1)
errs(1) = emax errs(1) = emax
dte = dt * errs(1)**(-0.2d+00) * errs(2)**(0.0505d+00) & dte = dt * min(fcmx, max(fcmn, &
* errs(3)**(-0.025d+00) fc * errs(1)**k1 * errs(2)**k2 * errs(3)**k3))
fcmx = facmax
else else
errs(1) = emax errs(1) = emax
dte = dt * errs(1)**(-0.2d+00) * errs(2)**(0.0505d+00) & dte = dt * min(fcmx, max(fcmn, &
* errs(3)**(-0.025d+00) fc * errs(1)**k1 * errs(2)**k2 * errs(3)**k3))
dt = dte dt = dte
fcmx = fac
nrej = nrej + 1 ! rejection count in the current step nrej = nrej + 1 ! rejection count in the current step
nrejections = nrejections + 1 nrejections = nrejections + 1
end if end if