EVOLUTION: Treat any unphysical state properly in RK3 method.

Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
Grzegorz Kowal 2021-11-10 09:31:30 -03:00
parent 3ffbf1a3c2
commit 9d3123eeed

View File

@ -1719,7 +1719,6 @@ module evolution
! "Numerical Recipes in Fortran",
! Cambridge University Press, Cambridge, 1992
!
!
!===============================================================================
!
subroutine evolve_rk3()
@ -1734,8 +1733,7 @@ module evolution
type(block_data), pointer :: pdata
integer :: status
real(kind=8) :: ds
real(kind=8) :: tm, dtm
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
@ -1746,105 +1744,127 @@ module evolution
call start_timer(imu)
#endif /* PROFILE */
status = 1
do while(status /= 0)
ds = 5.0d-01 * dt
!= 1st substep: U(1) = U(n) + dt F[U(n)]
!
ds = dt
tm = time + ds
dtm = ds
t = time + dt
pdata => list_data
do while (associated(pdata))
pdata => list_data
do while (associated(pdata))
call update_increment(pdata)
call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:))
call update_increment(pdata)
call update_sources(pdata, t, 0.0d+00, pdata%du(:,:,:,:))
pdata => pdata%next
end do
pdata => pdata%next
end do
call boundary_fluxes()
call boundary_fluxes()
pdata => list_data
do while (associated(pdata))
pdata => list_data
do while (associated(pdata))
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) + ds * pdata%du(:,:,:,:)
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) + dt * pdata%du(:,:,:,:)
pdata%u => pdata%uu(:,:,:,:,2)
pdata%u => pdata%uu(:,:,:,:,2)
pdata => pdata%next
end do
pdata => pdata%next
end do
call update_variables(tm, dtm, status)
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
tm = time + 0.5d+00 * dt
dtm = ds
ds = f22 * dt
t = time + 0.5d+00 * dt
pdata => list_data
do while (associated(pdata))
pdata => list_data
do while (associated(pdata))
call update_increment(pdata)
call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:))
call update_increment(pdata)
call update_sources(pdata, t, ds, pdata%du(:,:,:,:))
pdata => pdata%next
end do
pdata => pdata%next
end do
call boundary_fluxes()
call boundary_fluxes()
pdata => list_data
do while (associated(pdata))
pdata => list_data
do while (associated(pdata))
pdata%uu(:,:,:,:,2) = f21 * pdata%uu(:,:,:,:,1) &
+ f22 * pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:)
pdata%uu(:,:,:,:,2) = f21 * pdata%uu(:,:,:,:,1) &
+ f22 * pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:)
pdata => pdata%next
end do
pdata => pdata%next
end do
call update_variables(tm, dtm, status)
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
tm = time + dt
dtm = ds
ds = f32 * dt
t = time + dt
pdata => list_data
do while (associated(pdata))
pdata => list_data
do while (associated(pdata))
call update_increment(pdata)
call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:))
call update_increment(pdata)
call update_sources(pdata, t, ds, pdata%du(:,:,:,:))
pdata => pdata%next
end do
call boundary_fluxes()
pdata => list_data
do while (associated(pdata))
pdata%uu(:,:,:,:,2) = f31 * pdata%uu(:,:,:,:,1) &
+ f32 * pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:)
pdata%u => pdata%uu(:,:,:,:,1)
pdata => pdata%next
end do
if (ibp > 0) then
adecay(:) = exp(aglm(:) * cmax * dt)
pdata => list_data
do while (associated(pdata))
pdata%u(ibp,:,:,:) = adecay(pdata%meta%level) * pdata%u(ibp,:,:,:)
pdata => pdata%next
end do
end if
call update_variables(t, dt, status)
100 continue
if (status /= 0) dt = 2.5d-01 * dt
pdata => pdata%next
end do
call boundary_fluxes()
pdata => list_data
do while (associated(pdata))
pdata%uu(:,:,:,:,1) = f31 * pdata%uu(:,:,:,:,1) &
+ f32 * pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:)
pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2)
pdata%u => pdata%uu(:,:,:,:,1)
pdata => pdata%next
end do
if (ibp > 0) then
adecay(:) = exp(aglm(:) * cmax * dt)
pdata => list_data
do while (associated(pdata))
pdata%u(ibp,:,:,:) = adecay(pdata%meta%level) * pdata%u(ibp,:,:,:)
pdata => pdata%next
end do
end if
call update_variables(tm, dtm, status)
#ifdef PROFILE
call stop_timer(imu)
#endif /* PROFILE */