diff --git a/sources/evolution.F90 b/sources/evolution.F90 index e806006..e0cfbcb 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") @@ -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)" @@ -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)" @@ -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 @@ -1463,7 +1494,7 @@ module evolution type(block_data), pointer :: pdata integer :: status - real(kind=8) :: tm, dtm + real(kind=8) :: t ! !------------------------------------------------------------------------------- ! @@ -1471,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 */ @@ -1545,7 +1601,7 @@ module evolution type(block_data), pointer :: pdata integer :: status - real(kind=8) :: tm, dtm + real(kind=8) :: t, ds ! !------------------------------------------------------------------------------- ! @@ -1553,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 */ @@ -1645,7 +1719,6 @@ module evolution ! "Numerical Recipes in Fortran", ! Cambridge University Press, Cambridge, 1992 ! -! !=============================================================================== ! subroutine evolve_rk3() @@ -1660,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 @@ -1672,105 +1744,125 @@ module evolution call start_timer(imu) #endif /* PROFILE */ + status = 1 + + do while(status /= 0) + != 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 */ @@ -1809,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 ! !------------------------------------------------------------------------------- ! @@ -1818,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 */ @@ -1976,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 @@ -1996,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 */ @@ -2185,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 @@ -2202,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 @@ -2278,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 */