From 3117c19a75bbbed283a67df5f7b09d2e0359b3b0 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Wed, 10 Nov 2021 09:04:05 -0300 Subject: [PATCH 1/7] EVOLUTION: Improve initial time step estimation. In case the estimation returns an unphysical state, its time step is reduced and the estimation is repeated. If it fails for 10 consecutive times, the estimation is abandoned and the CFL time step is used. Signed-off-by: Grzegorz Kowal --- sources/evolution.F90 | 105 +++++++++++++++++++++++++++--------------- 1 file changed, 68 insertions(+), 37 deletions(-) diff --git a/sources/evolution.F90 b/sources/evolution.F90 index e806006..b4b8173 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -1172,22 +1172,24 @@ module evolution ! subroutine initialize_time_step() - use blocks , only : block_data, list_data, get_nleafs - use coordinates, only : nc => ncells, nb, ne, adx, ady + use blocks , only : block_data, list_data, get_nleafs + use coordinates , only : nc => ncells, nb, ne, adx, ady #if NDIMS == 3 - use coordinates, only : adz + use coordinates , only : adz #endif /* NDIMS == 3 */ - use equations , only : nf, maxspeed, cmax, cmax2 + use equations , only : nf, maxspeed, cmax, cmax2 + use iso_fortran_env, only : error_unit #ifdef MPI - use mpitools , only : reduce_maximum, reduce_sum + use mpitools , only : reduce_maximum, reduce_sum #endif /* MPI */ - use sources , only : viscosity, resistivity, update_sources + use sources , only : viscosity, resistivity, update_sources implicit none - type(block_data), pointer :: pdata + type(block_data), pointer :: pdata - integer :: mlev, status + logical :: flag + integer :: mlev, n, status real(kind=8) :: dx_min, fnorm, h0, h1 real(kind=8), dimension(3) :: d @@ -1198,6 +1200,8 @@ module evolution #endif /* NDIMS == 3 */ real(kind=8), parameter :: eps = tiny(cmax) + + character(len=*), parameter :: loc = 'EVOLUTION:initialize_time_step()' ! !------------------------------------------------------------------------------- ! @@ -1278,59 +1282,86 @@ module evolution h0 = 1.0d-06 end if - pdata => list_data - do while (associated(pdata)) + n = 10 + flag = .true. - pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) + h0 * pdata%du(:,:,:,:) + do while(flag) - pdata%u => pdata%uu(:,:,:,:,2) + pdata => list_data + do while (associated(pdata)) + + pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) + h0 * pdata%du(:,:,:,:) + + pdata%u => pdata%uu(:,:,:,:,2) + + pdata => pdata%next + end do + + call update_variables(time + h0, h0, status) + + if (status /= 0) then + h0 = 2.5d-01 * h0 + flag = n > 0 + n = n - 1 + else + flag = .false. + end if - pdata => pdata%next end do - call update_variables(time + h0, h0, status) + if (status /= 0) then - pdata => list_data - do while (associated(pdata)) + write(error_unit,"('[', a, ']: ', a)") trim(loc), & + "Could not estimate the initial time step due to " // & + "the error. Setting it to the CFL time step." + + dte = dt + + else + + pdata => list_data + do while (associated(pdata)) #if NDIMS == 3 - df(:,:,:,:) = pdata%du(:,nb:ne,nb:ne,nb:ne) - sc(:,:,:,:) = atol + rtol * abs(pdata%uu(:,nb:ne,nb:ne,nb:ne,1)) + df(:,:,:,:) = pdata%du(:,nb:ne,nb:ne,nb:ne) + sc(:,:,:,:) = atol + rtol * abs(pdata%uu(:,nb:ne,nb:ne,nb:ne,1)) #else /* NDIMS == 3 */ - df(:,:,:,:) = pdata%du(:,nb:ne,nb:ne, : ) - sc(:,:,:,:) = atol + rtol * abs(pdata%uu(:,nb:ne,nb:ne, : ,1)) + df(:,:,:,:) = pdata%du(:,nb:ne,nb:ne, : ) + sc(:,:,:,:) = atol + rtol * abs(pdata%uu(:,nb:ne,nb:ne, : ,1)) #endif /* NDIMS == 3 */ - call update_increment(pdata) - call update_sources(pdata, time + h0, 0.0d+00, pdata%du(:,:,:,:)) + call update_increment(pdata) + call update_sources(pdata, time + h0, 0.0d+00, pdata%du(:,:,:,:)) #if NDIMS == 3 - df(:,:,:,:) = df(:,:,:,:) - pdata%du(:,nb:ne,nb:ne,nb:ne) + df(:,:,:,:) = df(:,:,:,:) - pdata%du(:,nb:ne,nb:ne,nb:ne) #else /* NDIMS == 3 */ - df(:,:,:,:) = df(:,:,:,:) - pdata%du(:,nb:ne,nb:ne, : ) + df(:,:,:,:) = df(:,:,:,:) - pdata%du(:,nb:ne,nb:ne, : ) #endif /* NDIMS == 3 */ - d(3) = d(3) + sum((df(:,:,:,:) / sc(:,:,:,:))**2) + d(3) = d(3) + sum((df(:,:,:,:) / sc(:,:,:,:))**2) - pdata%u => pdata%uu(:,:,:,:,1) + pdata%u => pdata%uu(:,:,:,:,1) - pdata => pdata%next - end do + pdata => pdata%next + end do #ifdef MPI - call reduce_sum(d(3:3)) + call reduce_sum(d(3:3)) #endif /* MPI */ - d(3) = sqrt(fnorm * d(3)) / h0 + d(3) = sqrt(fnorm * d(3)) / h0 + + if (max(d(2), d(3)) > 1.0d-15) then + h1 = (1.0d-02 / max(d(2), d(3)))**(1.0d+00 / (order + 1)) + else + h1 = max(1.0d-06, 1.0d-03 * h0) + end if + + dte = min(1.0d+02 * h0, h1) + dt = min(dt, dte) - if (max(d(2), d(3)) > 1.0d-15) then - h1 = (1.0d-02 / max(d(2), d(3)))**(1.0d+00 / (order + 1)) - else - h1 = max(1.0d-06, 1.0d-03 * h0) end if - dte = min(1.0d+02 * h0, h1) - dt = min(dt, dte) - end if #ifdef PROFILE From d6b4b70aaba76351e6f17dc49cc6dac128116984 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Wed, 10 Nov 2021 09:18:19 -0300 Subject: [PATCH 2/7] EVOLUTION: Treat any unphysical state properly in EULER method. Signed-off-by: Grzegorz Kowal --- 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 */ From 3ffbf1a3c2b854610bbbb4cbc800321767895f72 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Wed, 10 Nov 2021 09:25:31 -0300 Subject: [PATCH 3/7] EVOLUTION: Treat any unphysical state properly in RK2 method. Signed-off-by: Grzegorz Kowal --- sources/evolution.F90 | 104 +++++++++++++++++++++++++----------------- 1 file changed, 61 insertions(+), 43 deletions(-) diff --git a/sources/evolution.F90 b/sources/evolution.F90 index cc03629..ba99d4e 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -1601,7 +1601,7 @@ module evolution type(block_data), pointer :: pdata integer :: status - real(kind=8) :: tm, dtm + real(kind=8) :: t, ds ! !------------------------------------------------------------------------------- ! @@ -1609,76 +1609,94 @@ module evolution call start_timer(imu) #endif /* PROFILE */ + status = 1 + + do while(status /= 0) + + ds = 5.0d-01 * dt + != 1st step: U(1) = U(n) + dt * F[U(n)] ! - tm = time + dt - dtm = 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, 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) + dt * 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(n+1) = 1/2 U(n) + 1/2 U(1) + 1/2 dt * F[U(1)] ! - tm = time + dt - dtm = 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, t, ds, pdata%du(:,:,:,:)) - call update_increment(pdata) - call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) + pdata => pdata%next + end do + + call boundary_fluxes() + + pdata => list_data + do while (associated(pdata)) + + pdata%uu(:,:,:,:,2) = 5.0d-01 * (pdata%uu(:,:,:,:,1) & + + 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, ds, 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) = 0.5d+00 * (pdata%uu(:,:,:,:,1) & - + pdata%uu(:,:,:,:,2) & - + 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 */ From 9d3123eeed1efc5c621766ba45323676714aba75 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Wed, 10 Nov 2021 09:31:30 -0300 Subject: [PATCH 4/7] EVOLUTION: Treat any unphysical state properly in RK3 method. Signed-off-by: Grzegorz Kowal --- sources/evolution.F90 | 146 ++++++++++++++++++++++++------------------ 1 file changed, 83 insertions(+), 63 deletions(-) diff --git a/sources/evolution.F90 b/sources/evolution.F90 index ba99d4e..d2c78bf 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -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 */ From c42310b3cb1d7cb771cf6f7edecc7820fb091ddc Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Wed, 10 Nov 2021 09:37:21 -0300 Subject: [PATCH 5/7] 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 */ From 1d22adb2fababad8f574f55f43a418ab6930f3b9 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Wed, 10 Nov 2021 09:51:06 -0300 Subject: [PATCH 6/7] EVOLUTION: Treat any unphysical state properly in RK3.5 method. Signed-off-by: Grzegorz Kowal --- sources/evolution.F90 | 236 +++++++++++++++++++++++------------------- 1 file changed, 130 insertions(+), 106 deletions(-) diff --git a/sources/evolution.F90 b/sources/evolution.F90 index 3ff2890..f313e06 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -259,7 +259,7 @@ module evolution evolve => evolve_ssprk35 order = 3 - registers = 2 + registers = 3 cfl = 2.65062919143939d+00 * cfl name_int = "3rd order SSPRK(5,3)" @@ -2090,8 +2090,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 :: b1 = 3.77268915331368d-01 real(kind=8), parameter :: b3 = 2.42995220537396d-01 @@ -2110,159 +2109,184 @@ module evolution call start_timer(imu) #endif /* PROFILE */ + status = 1 + + do while(status /= 0) + != 1st step: U(1) = U(n) + b1 dt F[U(n)] ! - ds = b1 * dt - tm = time + ds - dtm = ds + ds = b1 * 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) + ds * 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(1) + b1 dt F[U(1)] ! - tm = time + 2.0d+00 * ds - dtm = ds + t = time + 2.0d+00 * 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) = pdata%uu(:,:,:,:,2) + ds * 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) = a31 U(n) + a33 U(2) + b3 dt F[U(2)] ! - ds = b3 * dt - tm = time + (2.0d+00 * a33 * b1 + b3) * dt - dtm = ds + ds = b3 * dt + t = time + (2.0d+00 * a33 * b1 + b3) * 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) = a31 * pdata%uu(:,:,:,:,1) & - + a33 * pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:) + pdata%uu(:,:,:,:,2) = a31 * pdata%uu(:,:,:,:,1) & + + a33 * 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 != 4th step: U(4) = a41 U(n) + a44 U(3) + b4 dt F[U(3)] ! - ds = b4 * dt - tm = time + ((2.0d+00 * b1 * a33 + b3) * a44 + b4) * dt - dtm = ds + ds = b4 * dt + t = time + ((2.0d+00 * b1 * a33 + b3) * a44 + b4) * 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(:,:,:,:,3) = a41 * pdata%uu(:,:,:,:,1) & + + a44 * pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:) + + pdata%u => pdata%uu(:,:,:,:,3) + + pdata => pdata%next + end do + + call update_variables(t, ds, status) + + if (status /= 0) go to 100 + +!= the final step: U(n+1) = a53 U(2) + a55 U(4) + b5 dt F[U(4)] +! + ds = b5 * dt + t = time + dt + + pdata => list_data + do while (associated(pdata)) + + 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) = a53 * pdata%uu(:,:,:,:,2) & + + a55 * pdata%uu(:,:,:,:,3) + ds * pdata%du(:,:,:,:) + + pdata%u => pdata%uu(:,:,:,:,2) + + 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) = a41 * pdata%uu(:,:,:,:,1) & - + a44 * pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:) + pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2) pdata%u => pdata%uu(:,:,:,:,1) pdata => pdata%next end do - call update_variables(tm, dtm, status) - -!= the final step: U(n+1) = a53 U(2) + a55 U(4) + b5 dt F[U(4)] -! - ds = b5 * dt - tm = time + dt - dtm = ds - - pdata => list_data - do while (associated(pdata)) - - call update_increment(pdata) - call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) - - pdata => pdata%next - end do - - call boundary_fluxes() - - pdata => list_data - do while (associated(pdata)) - - pdata%uu(:,:,:,:,1) = a53 * pdata%uu(:,:,:,:,2) & - + a55 * pdata%uu(:,:,:,:,1) + 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(tm, dtm, status) - #ifdef PROFILE call stop_timer(imu) #endif /* PROFILE */ From 7578206f83b4cadaf1a7af1086ba5427995048a4 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Wed, 10 Nov 2021 10:01:34 -0300 Subject: [PATCH 7/7] EVOLUTION: Treat any unphysical state properly in RK4.10 method. Signed-off-by: Grzegorz Kowal --- sources/evolution.F90 | 168 ++++++++++++++++++++++++------------------ 1 file changed, 96 insertions(+), 72 deletions(-) diff --git a/sources/evolution.F90 b/sources/evolution.F90 index f313e06..e0cfbcb 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -267,7 +267,7 @@ module evolution evolve => evolve_ssprk4_10 order = 4 - registers = 2 + registers = 3 cfl = 6.0d+00 * cfl name_int = "Optimal 4th order SSPRK(10,4)" @@ -2323,8 +2323,8 @@ module evolution type(block_data), pointer :: pdata - integer :: n, status - real(kind=8) :: tm, dtm + integer :: i, status + real(kind=8) :: t real(kind=8), dimension(9) :: ds @@ -2340,73 +2340,109 @@ module evolution ds(:) = c(:) * dt + status = 1 + + do while(status /= 0) + != 1st step: U(2) = U(1) ! - pdata => list_data - do while (associated(pdata)) + pdata => list_data + do while (associated(pdata)) - pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) + pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) + pdata%uu(:,:,:,:,3) = pdata%uu(:,:,:,:,1) - pdata => pdata%next - end do + pdata%u => pdata%uu(:,:,:,:,2) + + pdata => pdata%next + end do != 2nd step: U(1) = [1 + dt/6 L] U(1), for i = 1, ..., 5 ! - do n = 1, 5 + do i = 1, 5 - tm = time + ds(n) - dtm = ds(1) + t = time + ds(i) - 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(:,:,:,:,1) = pdata%uu(:,:,:,:,1) + dtm * pdata%du(:,:,:,:) + pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + ds(1) * pdata%du(:,:,:,:) - pdata => pdata%next - end do + pdata => pdata%next + end do - call update_variables(tm, dtm, status) + call update_variables(t, ds(1), status) - end do ! n = 1, 5 + if (status /= 0) go to 100 + + end do ! i = 1, 5 != 3rd step: U(2) = U(2)/25 + 9/25 U(1), U(1) = 15 U(2) - 5 U(1) ! - pdata => list_data - do while (associated(pdata)) + pdata => list_data + do while (associated(pdata)) - pdata%uu(:,:,:,:,2) = (pdata%uu(:,:,:,:,2) & - + 9.0d+00 * pdata%uu(:,:,:,:,1)) / 2.5d+01 - pdata%uu(:,:,:,:,1) = 1.5d+01 * pdata%uu(:,:,:,:,2) & - - 5.0d+00 * pdata%uu(:,:,:,:,1) + pdata%uu(:,:,:,:,3) = (pdata%uu(:,:,:,:,3) & + + 9.0d+00 * pdata%uu(:,:,:,:,2)) / 2.5d+01 + pdata%uu(:,:,:,:,2) = 1.5d+01 * pdata%uu(:,:,:,:,3) & + - 5.0d+00 * pdata%uu(:,:,:,:,2) - pdata => pdata%next - end do + pdata => pdata%next + end do != 4th step: U(1) = [1 + dt/6 L] U(1), for i = 6, ..., 9 ! ! integrate the intermediate steps ! - do n = 6, 9 + do i = 6, 9 - tm = time + ds(n) - dtm = ds(1) + t = time + ds(i) + + pdata => list_data + do while (associated(pdata)) + + call update_increment(pdata) + call update_sources(pdata, t, ds(1), pdata%du(:,:,:,:)) + + pdata => pdata%next + end do + + call boundary_fluxes() + + pdata => list_data + do while (associated(pdata)) + + pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + ds(1) * pdata%du(:,:,:,:) + + pdata => pdata%next + end do + + call update_variables(t, ds(1), status) + + if (status /= 0) go to 100 + + end do ! i = 6, 9 + +!= the final step: U(n+1) = U(2) + 3/5 U(1) + 1/10 dt F[U(1)] +! + t = time + dt pdata => list_data do while (associated(pdata)) call update_increment(pdata) - call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) + call update_sources(pdata, t, 1.0d-01 * dt, pdata%du(:,:,:,:)) pdata => pdata%next end do @@ -2416,55 +2452,43 @@ module evolution pdata => list_data do while (associated(pdata)) - pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,1) + dtm * pdata%du(:,:,:,:) + pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,3) & + + 6.0d-01 * pdata%uu(:,:,:,:,2) & + + 1.0d-01 * dt * pdata%du(:,:,:,:) pdata => pdata%next end do - call update_variables(tm, dtm, status) + if (ibp > 0) then + adecay(:) = exp(aglm(:) * cmax * dt) - end do ! n = 6, 9 + pdata => list_data + do while (associated(pdata)) -!= the final step: U(n+1) = U(2) + 3/5 U(1) + 1/10 dt F[U(1)] -! - tm = time + dt - dtm = dt / 1.0d+01 + 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 + + end do pdata => list_data do while (associated(pdata)) - call update_increment(pdata) - call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) + pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2) + + pdata%u => pdata%uu(:,:,:,:,1) pdata => pdata%next end do - call boundary_fluxes() - - pdata => list_data - do while (associated(pdata)) - - pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2) & - + 6.0d-01 * pdata%uu(:,:,:,:,1) & - + dtm * 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(tm, dtm, status) - #ifdef PROFILE call stop_timer(imu) #endif /* PROFILE */