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

EVOLUTION: Improve the conservation in the SSPRK3 and SSPRK(4,3) methods.

Browse Source 
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
Grzegorz Kowal 2023-08-29 10:43:42 -03:00
parent 91d3583a78
commit dfe7310408
1 changed files with 13 additions and 17 deletions

30
sources/evolution.F90

@ -1727,9 +1727,6 @@ module evolution
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
@ -1771,8 +1768,8 @@ module evolution
!= 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
ds = dt / 4.0d+00
t = time + 5.0d-01 * dt
!$omp parallel do default(shared) private(pdata)
do l = 1, n
@ -1791,8 +1788,9 @@ module evolution
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,2) = f21 * pdata%uu(:,:,:,:,1) &
+ f22 * pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:)
pdata%uu(:,:,:,:,2) = (3.0d+00 * pdata%uu(:,:,:,:,1) &
+ pdata%uu(:,:,:,:,2) &
+ ds * pdata%du(:,:,:,:)) / 4.0d+00
end do
!$omp end parallel do
@ -1802,7 +1800,7 @@ module evolution
!= 3rd step: U(n+1) = 1/3 U(n) + 2/3 U(2) + 2/3 dt F[U(2)]
!
ds = f32 * dt
ds = 2.0d+00 * dt / 3.0d+00
t = time + dt
!$omp parallel do default(shared) private(pdata)
@ -1822,8 +1820,9 @@ module evolution
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,2) = f31 * pdata%uu(:,:,:,:,1) &
+ f32 * pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:)
pdata%uu(:,:,:,:,2) = ( pdata%uu(:,:,:,:,1) &
+ 2.0d+00 * (pdata%uu(:,:,:,:,2) &
+ ds * pdata%du(:,:,:,:))) / 3.0d+00
pdata%u => pdata%uu(:,:,:,:,1)
end do
@ -1845,7 +1844,7 @@ module evolution
100 continue
if (status /= 0) dt = 2.5d-01 * dt
if (status /= 0) dt = dt / 4.0d+00
end do
@ -1895,9 +1894,6 @@ module evolution
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
@ -1988,9 +1984,9 @@ module evolution
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,2) = b1 * pdata%uu(:,:,:,:,1) &
+ b2 * pdata%uu(:,:,:,:,2) &
+ b3 * dt * pdata%du(:,:,:,:)
pdata%uu(:,:,:,:,2) = (2.0d+00 * pdata%uu(:,:,:,:,1) &
+ pdata%uu(:,:,:,:,2) &
+ 5.0d-01 * dt * pdata%du(:,:,:,:)) / 3.0d+00
end do
!$omp end parallel do