From d6b4b70aaba76351e6f17dc49cc6dac128116984 Mon Sep 17 00:00:00 2001
From: Grzegorz Kowal <grzegorz@amuncode.org>
Date: Wed, 10 Nov 2021 09:18:19 -0300
Subject: [PATCH] EVOLUTION: Treat any unphysical state properly in EULER
method.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
---
sources/evolution.F90 | 77 ++++++++++++++++++++++++++++---------------
1 file changed, 51 insertions(+), 26 deletions(-)
diff --git a/sources/evolution.F90 b/sources/evolution.F90
index b4b8173..cc03629 100644
--- a/sources/evolution.F90
+++ b/sources/evolution.F90
@@ -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 */