EVOLUTION: Treat any unphysical state properly in RK3.4 method.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
parent
9d3123eeed
commit
c42310b3cb
@ -1748,8 +1748,6 @@ module evolution
|
||||
|
||||
do while(status /= 0)
|
||||
|
||||
ds = 5.0d-01 * dt
|
||||
|
||||
!= 1st substep: U(1) = U(n) + dt F[U(n)]
|
||||
!
|
||||
t = time + dt
|
||||
@ -1903,8 +1901,7 @@ module evolution
|
||||
type(block_data), pointer :: pdata
|
||||
|
||||
integer :: status
|
||||
real(kind=8) :: ds
|
||||
real(kind=8) :: tm, dtm
|
||||
real(kind=8) :: t, ds
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
@ -1912,126 +1909,149 @@ module evolution
|
||||
call start_timer(imu)
|
||||
#endif /* PROFILE */
|
||||
|
||||
status = 1
|
||||
|
||||
do while(status /= 0)
|
||||
|
||||
!= 1st step: U(1) = U(n) + 1/2 dt F[U(n)]
|
||||
!
|
||||
ds = dt / 2.0d+00
|
||||
tm = time + ds
|
||||
dtm = ds
|
||||
ds = 5.0d-01 * dt
|
||||
t = time + ds
|
||||
|
||||
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) + dtm * pdata%du(:,:,:,:)
|
||||
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) + ds * 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, ds, status)
|
||||
|
||||
if (status /= 0) go to 100
|
||||
|
||||
!= 2nd step: U(2) = U(2) + 1/2 dt F[U(1)]
|
||||
!
|
||||
tm = time + 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
|
||||
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(:,:,:,:,2) + dtm * pdata%du(:,:,:,:)
|
||||
pdata%uu(:,:,:,:,2) = 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(3) = 2/3 U(n) + 1/3 (U(2) + 1/2 dt F[U(2)])
|
||||
!
|
||||
tm = time + ds
|
||||
t = time + ds
|
||||
|
||||
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) = (2.0d+00 * pdata%uu(:,:,:,:,1) &
|
||||
+ pdata%uu(:,:,:,:,2) &
|
||||
+ dtm * pdata%du(:,:,:,:)) / 3.0d+00
|
||||
pdata%uu(:,:,:,:,2) = (2.0d+00 * pdata%uu(:,:,:,:,1) &
|
||||
+ pdata%uu(:,:,:,:,2) &
|
||||
+ ds * pdata%du(:,:,:,:)) / 3.0d+00
|
||||
|
||||
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
|
||||
|
||||
!= the final step: U(n+1) = U(3) + 1/2 dt F[U(3)]
|
||||
!
|
||||
tm = time + 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) = pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:)
|
||||
|
||||
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) = 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 */
|
||||
|
Loading…
x
Reference in New Issue
Block a user