EVOLUTION: Treat any unphysical state properly in EULER method.

Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
Grzegorz Kowal 2021-11-10 09:18:19 -03:00
parent 3117c19a75
commit d6b4b70aab

View File

@ -230,7 +230,7 @@ module evolution
evolve => evolve_euler
order = 1
registers = 1
registers = 2
name_int = "1st order Euler"
case ("rk2", "RK2")
@ -1494,7 +1494,7 @@ module evolution
type(block_data), pointer :: pdata
integer :: status
real(kind=8) :: tm, dtm
real(kind=8) :: t
!
!-------------------------------------------------------------------------------
!
@ -1502,44 +1502,69 @@ module evolution
call start_timer(imu)
#endif /* PROFILE */
tm = time + dt
dtm = dt
status = 1
pdata => list_data
do while (associated(pdata))
do while(status /= 0)
call update_increment(pdata)
call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:))
t = time + dt
pdata => list_data
do while (associated(pdata))
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1)
pdata%u => pdata%uu(:,:,:,:,2)
pdata => pdata%next
end do
pdata => list_data
do while (associated(pdata))
call update_increment(pdata)
call update_sources(pdata, t, dt, pdata%du(:,:,:,:))
pdata => pdata%next
end do
call boundary_fluxes()
pdata => list_data
do while (associated(pdata))
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + dt * 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)
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(:,:,:,:,1) + dt * 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 */