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

EVOLUTION: Treat any unphysical state properly in SSPRK3(2)4 method.

Browse Source 
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
Grzegorz Kowal 2021-11-09 16:09:32 -03:00
parent a52da2830d
commit 6dcd400219
1 changed files with 26 additions and 10 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

36
sources/evolution.F90
View File

@ -2222,11 +2222,12 @@ module evolution
!
subroutine evolve_ssp324()
use blocks , only : block_data, list_data
use boundaries , only : boundary_fluxes
use coordinates, only : nb, ne
use equations , only : errors, ibp, cmax
use sources , only : update_sources
use blocks , only : block_data, list_data
use boundaries , only : boundary_fluxes
use coordinates , only : nb, ne
use equations , only : errors, ibp, cmax
use iso_fortran_env, only : error_unit
use sources , only : update_sources
implicit none
@ -2238,7 +2239,9 @@ module evolution
real(kind=8), parameter :: onethird = 1.0d+00 / 3.0d+00, &
twothird = 2.0d+00 / 3.0d+00
!
character(len=*), parameter :: loc = 'EVOLUTION:evolve_ssp324()'
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
@ -2301,6 +2304,7 @@ module evolution
end do
call update_variables(tm, dtm, status)
if (status /= 0) go to 100
end do ! i = 1, 3
@ -2318,6 +2322,7 @@ module evolution
end do
call update_variables(tm, dtm, status)
if (status /= 0) go to 100
!= 5th step: U(1) = U(1) + ½ dt F[U(1)] <- 3ʳᵈ-order candidate
!
@ -2344,6 +2349,7 @@ module evolution
end do
call update_variables(tm, dtm, status)
if (status /= 0) go to 100
!= 6th step: U(2) = ½ (U(1) + U(2)) <- 2ⁿᵈ-order approximation
!
@ -2385,16 +2391,21 @@ module evolution
fc = 1.0d+00 + chi * atan((fc - 1.0d+00) / chi)
dte = dt * fc
if (fc > fac .or. nrej >= mrej) then
100 continue
if ((fc > fac .or. nrej >= mrej) .and. status == 0) then
test = .false.
errs(3) = errs(2)
errs(2) = errs(1)
else
dt = dte
nrej = nrej + 1 ! rejection count in the current step
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
@ -2409,6 +2420,11 @@ module evolution
end do
call update_variables(tm, dtm, status)
if (status /= 0) then
write(error_unit,"('[', a, ']: ', a)") trim(loc), &
"Possible bug: the revert to the previous state " //&
"found it unphysical."
end if
end if