From c42310b3cb1d7cb771cf6f7edecc7820fb091ddc Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Wed, 10 Nov 2021 09:37:21 -0300 Subject: [PATCH] EVOLUTION: Treat any unphysical state properly in RK3.4 method. Signed-off-by: Grzegorz Kowal --- sources/evolution.F90 | 168 +++++++++++++++++++++++------------------- 1 file changed, 94 insertions(+), 74 deletions(-) diff --git a/sources/evolution.F90 b/sources/evolution.F90 index d2c78bf..3ff2890 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -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 */