diff --git a/sources/evolution.F90 b/sources/evolution.F90 index a9b4d3c..725d37c 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -70,6 +70,7 @@ module evolution logical , save :: error_control = .false. character(len=255), save :: name_int = "" character(len=255), save :: name_norm = "" + integer , save :: order = 2 integer , save :: stages = 2 integer , save :: registers = 1 real(kind=8) , save :: cfl = 5.0d-01 @@ -118,7 +119,7 @@ module evolution ! declare public subroutines ! public :: initialize_evolution, finalize_evolution, print_evolution - public :: advance, new_time_step + public :: initialize_time_step, advance, new_time_step ! declare public variables ! @@ -222,12 +223,14 @@ module evolution case ("euler", "EULER", "rk1", "RK1") evolve => evolve_euler + order = 1 registers = 1 name_int = "1st order Euler" case ("rk2", "RK2") evolve => evolve_rk2 + order = 2 registers = 2 name_int = "2nd order Runge-Kutta" @@ -235,6 +238,7 @@ module evolution error_control = .true. evolve => evolve_ssprk2_m + order = 2 registers = 3 stages = max(2, min(9, stages)) cfl = (stages - 1) * cfl @@ -243,12 +247,14 @@ module evolution case ("rk3", "RK3") evolve => evolve_rk3 + order = 3 registers = 2 name_int = "3rd order Runge-Kutta" case ("rk3.4", "RK3.4", "ssprk(4,3)", "SSPRK(4,3)") evolve => evolve_ssprk34 + order = 3 registers = 2 cfl = 2.0d+00 * cfl name_int = "3rd order SSPRK(4,3)" @@ -256,6 +262,7 @@ module evolution case ("rk3.5", "RK3.5", "ssprk(5,3)", "SSPRK(5,3)") evolve => evolve_ssprk35 + order = 3 registers = 2 cfl = 2.65062919143939d+00 * cfl name_int = "3rd order SSPRK(5,3)" @@ -264,6 +271,7 @@ module evolution error_control = .true. evolve => evolve_ssp324 + order = 3 registers = 3 stages = 4 cfl = 2.0d+00 * cfl @@ -281,6 +289,7 @@ module evolution error_control = .true. evolve => evolve_ssprk32m + order = 3 registers = 4 n = 2 do while(n**2 <= stages) @@ -312,6 +321,7 @@ module evolution case ("rk4.10", "RK4.10", "ssprk(10,4)", "SSPRK(10,4)") evolve => evolve_ssprk4_10 + order = 4 registers = 2 cfl = 6.0d+00 * cfl name_int = "Optimal 4th order SSPRK(10,4)" @@ -591,6 +601,197 @@ module evolution ! !=============================================================================== ! +! subroutine INITIALIZE_TIME_STEP: +! ------------------------------- +! +! Subroutine estimates the initial time step. +! +! References: +! +! [1] Courant, R., Friedrichs, K., Lewy, H., +! "Über die partiellen Differenzengleichungen der mathematischen Physik", +! Mathematische Annalen, 1928, volume 100, pages 32–74, +! DOI: 10.1007/BF01448839 +! [2] Hairer, E., Nørsett, S. P., Wanner, G., +! "Solving Ordinary Differential Equations I: Nonstiff Problems", +! Springer-Verlag Berlin Heidelberg, 1993, +! ISBN: 978-3-540-56670-0, DOI: 10.1007/978-3-540-78862-1 +! +!=============================================================================== +! + subroutine initialize_time_step() + + 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 +#endif /* NDIMS == 3 */ + use equations , only : nf, maxspeed, cmax, cmax2 +#ifdef MPI + use mpitools , only : reduce_maximum, reduce_sum +#endif /* MPI */ + use sources , only : viscosity, resistivity, update_sources + + implicit none + + type(block_data), pointer :: pdata + + integer :: mlev + real(kind=8) :: dx_min, fnorm, h0, h1 + real(kind=8), dimension(3) :: d + +#if NDIMS == 3 + real(kind=8), dimension(nf,nc,nc,nc) :: sc, df +#else /* NDIMS == 3 */ + real(kind=8), dimension(nf,nc,nc, 1) :: sc, df +#endif /* NDIMS == 3 */ + + real(kind=8), parameter :: eps = tiny(cmax) +! +!------------------------------------------------------------------------------- +! +#ifdef PROFILE + call start_timer(imi) +#endif /* PROFILE */ + + cmax = eps + mlev = 1 + + pdata => list_data + do while (associated(pdata)) + + mlev = max(mlev, pdata%meta%level) + cmax = max(cmax, maxspeed(pdata%q(:,:,:,:))) + + pdata => pdata%next + end do + +#ifdef MPI + call reduce_maximum(cmax) + call reduce_maximum(mlev) +#endif /* MPI */ + + cmax2 = cmax * cmax + +#if NDIMS == 2 + dx_min = min(adx(mlev), ady(mlev)) +#endif /* NDIMS == 2 */ +#if NDIMS == 3 + dx_min = min(adx(mlev), ady(mlev), adz(mlev)) +#endif /* NDIMS == 3 */ + +! determine the time step due to the stability +! + dth = dx_min / max(cmax & + + 2.0d+00 * max(viscosity, resistivity) / dx_min, eps) + dt = cfl * dth + +! determine the time step due to the error +! + if (error_control) then + + fnorm = 1.0d+00 / (get_nleafs() * nc**NDIMS) + + d(:) = 0.0d+00 + + pdata => list_data + do while (associated(pdata)) + + call update_increment(pdata) + call update_sources(pdata, 0.0d+00, 0.0d+00, pdata%du(:,:,:,:)) + +#if NDIMS == 3 + sc(:,:,:,:) = atol + rtol * abs(pdata%uu(:,nb:ne,nb:ne,nb:ne,1)) + + d(1) = d(1) + sum((pdata%uu(:,nb:ne,nb:ne,nb:ne,1) / sc(:,:,:,:))**2) + d(2) = d(2) + sum((pdata%du(:,nb:ne,nb:ne,nb:ne) / sc(:,:,:,:))**2) +#else /* NDIMS == 3 */ + sc(:,:,:,:) = atol + rtol * abs(pdata%uu(:,nb:ne,nb:ne, : ,1)) + + d(1) = d(1) + sum((pdata%uu(:,nb:ne,nb:ne, : ,1) / sc(:,:,:,:))**2) + d(2) = d(2) + sum((pdata%du(:,nb:ne,nb:ne, : ) / sc(:,:,:,:))**2) +#endif /* NDIMS == 3 */ + + pdata => pdata%next + end do + +#ifdef MPI + call reduce_sum(d(1:2)) +#endif /* MPI */ + + d(1:2) = sqrt(fnorm * d(1:2)) + + if (minval(d(1:2)) >= 1.0d-05) then + h0 = 1.0d-02 * d(1) / d(2) + else + h0 = 1.0d-06 + end if + + 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) + + 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)) +#else /* NDIMS == 3 */ + 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(:,:,:,:)) + +#if NDIMS == 3 + df(:,:,:,:) = df(:,:,:,:) - pdata%du(:,nb:ne,nb:ne,nb:ne) +#else /* NDIMS == 3 */ + df(:,:,:,:) = df(:,:,:,:) - pdata%du(:,nb:ne,nb:ne, : ) +#endif /* NDIMS == 3 */ + d(3) = d(3) + sum((df(:,:,:,:) / sc(:,:,:,:))**2) + + pdata%u => pdata%uu(:,:,:,:,1) + + pdata => pdata%next + end do + +#ifdef MPI + call reduce_sum(d(3:3)) +#endif /* MPI */ + + 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) + + end if + +#ifdef PROFILE + call stop_timer(imi) +#endif /* PROFILE */ + +!------------------------------------------------------------------------------- +! + end subroutine initialize_time_step +! +!=============================================================================== +! ! subroutine NEW_TIME_STEP: ! ------------------------ !