From 3d13d30e60639b301b8f8db5d56e803ce8475f47 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Fri, 28 Aug 2020 21:30:56 -0300 Subject: [PATCH 01/20] EQUATIONS: Add vector to store integration errors of variables. Signed-off-by: Grzegorz Kowal --- sources/equations.F90 | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/sources/equations.F90 b/sources/equations.F90 index b1315a5..5ef6b33 100644 --- a/sources/equations.F90 +++ b/sources/equations.F90 @@ -221,6 +221,10 @@ module equations ! logical, dimension(:), allocatable :: positive +! the array of maximum errors for conservative variables +! + real(kind=8), dimension(:), allocatable :: errors + ! by default everything is private ! private @@ -250,7 +254,7 @@ module equations public :: eqsys, eos public :: pvars, cvars public :: qpbnd - public :: ivars, positive + public :: ivars, positive, errors !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! @@ -1060,6 +1064,12 @@ module equations end if +! allocate an array for errors +! + allocate(errors(nf), stat = status) + + if (status == 0) errors(:) = 0.0d+00 + ! proceed if everything is fine ! if (status == 0) then @@ -1252,6 +1262,10 @@ module equations ! if (allocated(positive)) deallocate(positive, stat = status) +! deallocate the array for errors +! + if (allocated(errors)) deallocate(errors, stat = status) + #ifdef PROFILE ! stop accounting time for module initialization/finalization ! From 0d47d35104b3b8d5b4f982f01e7848c6b792331b Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Fri, 28 Aug 2020 21:41:04 -0300 Subject: [PATCH 02/20] EVOLUTION: Add time step resulting from the error limit. Signed-off-by: Grzegorz Kowal --- sources/evolution.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/sources/evolution.F90 b/sources/evolution.F90 index 14f9021..5e7f0a0 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -75,6 +75,7 @@ module evolution real(kind=8) , save :: time = 0.0d+00 real(kind=8) , save :: dt = 1.0d+00 real(kind=8) , save :: dtn = 1.0d+00 + real(kind=8) , save :: dte = 0.0d+00 ! the increment array ! @@ -91,7 +92,7 @@ module evolution ! declare public variables ! - public :: step, time, dt, dtn, cfl, glm_alpha + public :: step, time, dt, dtn, dte, cfl, glm_alpha !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! From 51d883b5dbc2874d5a5ca8d226f1a43d279de8b2 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Fri, 28 Aug 2020 21:46:29 -0300 Subject: [PATCH 03/20] INTEGRALS: Store integration errors in a separate file. Signed-off-by: Grzegorz Kowal --- sources/integrals.F90 | 61 +++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 59 insertions(+), 2 deletions(-) diff --git a/sources/integrals.F90 b/sources/integrals.F90 index 8dc78e0..5d344ad 100644 --- a/sources/integrals.F90 +++ b/sources/integrals.F90 @@ -51,16 +51,22 @@ module integrals ! ! funit - a file handler to the integrals file; ! sunit - a file handler to the statistics file; +! eunit - a file handler to the errors file; ! iintd - the number of steps between subsequent intervals storing; ! integer(kind=4), save :: funit = 11 integer(kind=4), save :: sunit = 12 + integer(kind=4), save :: eunit = 13 integer(kind=4), save :: iintd = 1 ! flag indicating if the files are actually written to disk ! logical, save :: stored = .false. +! the format string +! + character(len=32), save :: efmt + ! by default everything is private ! private @@ -99,6 +105,7 @@ module integrals ! import external variables and subroutines ! + use equations , only : pvars, nf use parameters, only : get_parameter ! local variables are not implicit by default @@ -113,7 +120,7 @@ module integrals ! local variables ! - character(len=32) :: fname, append + character(len=32) :: fname, append, stmp logical :: flag ! !------------------------------------------------------------------------------- @@ -239,6 +246,53 @@ module integrals , 'mean(Malf)', 'min(Malf)', 'max(Malf)' write(sunit,"('#')") +! depending on the append parameter, choose the right file initialization for +! the file to store integration errors +! + append = "off" + call get_parameter("errors_append", append) + select case(trim(append)) + case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES") + write(fname, "('errors.dat')") + inquire(file = fname, exist = flag) + case default + write(fname, "('errors_',i2.2,'.dat')") irun + flag = .false. + end select + +! check if the file exists; if not, create a new one, otherwise move to the end +! + if (flag .and. irun > 1) then +#ifdef __INTEL_COMPILER + open(newunit = eunit, file = fname, form = 'formatted', status = 'old' & + , position = 'append', buffered = 'yes') +#else /* __INTEL_COMPILER */ + open(newunit = eunit, file = fname, form = 'formatted', status = 'old' & + , position = 'append') +#endif /* __INTEL_COMPILER */ + write(funit,"('#')") + else +#ifdef __INTEL_COMPILER + open(newunit = eunit, file = fname, form = 'formatted' & + , status = 'replace', buffered = 'yes') +#else /* __INTEL_COMPILER */ + open(newunit = eunit, file = fname, form = 'formatted' & + , status = 'replace') +#endif /* __INTEL_COMPILER */ + end if + +! write the integral file header +! + write(stmp,"(i9)") nf + 4 + efmt = "('#',a8," // trim(adjustl(stmp)) // "(1x,a18))" + write(eunit,efmt) 'step', 'time', 'dt_cfl', 'dt_err', & + 'max_err', pvars(1:nf) + write(eunit,"('#')") + +! prepare format for errors +! + efmt = "(i9," // trim(adjustl(stmp)) // "(1x,1es18.8e3))" + end if ! store #ifdef PROFILE @@ -291,6 +345,7 @@ module integrals if (stored) then close(funit) close(sunit) + close(eunit) end if #ifdef PROFILE @@ -324,7 +379,8 @@ module integrals use equations , only : idn, ipr, ivx, ivz, ibx, iby, ibz, ibp use equations , only : ien, imx, imy, imz use equations , only : magnetized, adiabatic_index, csnd - use evolution , only : step, time, dtn + use equations , only : errors, nf + use evolution , only : step, time, dtn, dte use forcing , only : einj, rinj, arms #ifdef MPI use mpitools , only : reduce_minimum, reduce_maximum, reduce_sum @@ -593,6 +649,7 @@ module integrals , avarr(5), mnarr(5), mxarr(5) & , avarr(6), mnarr(6), mxarr(6) & , avarr(7), mnarr(7), mxarr(7) + write(eunit,efmt) step, time, dtn, dte, maxval(errors(:)), errors(:) end if #ifdef PROFILE From b9e9ce1dafc818bf6a100f077871570a031477e5 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Fri, 28 Aug 2020 22:41:35 -0300 Subject: [PATCH 04/20] INTEGRALS: Remove unused variable from store_integrals(). Signed-off-by: Grzegorz Kowal --- sources/integrals.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sources/integrals.F90 b/sources/integrals.F90 index 5d344ad..dc9806e 100644 --- a/sources/integrals.F90 +++ b/sources/integrals.F90 @@ -379,7 +379,7 @@ module integrals use equations , only : idn, ipr, ivx, ivz, ibx, iby, ibz, ibp use equations , only : ien, imx, imy, imz use equations , only : magnetized, adiabatic_index, csnd - use equations , only : errors, nf + use equations , only : errors use evolution , only : step, time, dtn, dte use forcing , only : einj, rinj, arms #ifdef MPI From 100fcf12f2ca30769daed839b5218c19112ce2a7 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Fri, 28 Aug 2020 22:42:26 -0300 Subject: [PATCH 05/20] EVOLUTION: Add error estimation in evolve_ssprk3_m(). Signed-off-by: Grzegorz Kowal --- sources/evolution.F90 | 208 +++++++++++++++++------------------------- 1 file changed, 83 insertions(+), 125 deletions(-) diff --git a/sources/evolution.F90 b/sources/evolution.F90 index 5e7f0a0..5320a43 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -1853,42 +1853,34 @@ module evolution ! subroutine evolve_ssprk3_m() -! include external variables -! - use blocks , only : block_data, list_data - use equations, only : ibp - use sources , only : update_sources + use blocks , only : block_data, list_data, get_dblocks + use coordinates, only : nc => ncells, nb, ne + use equations , only : errors, nf, ibp +#ifdef MPI + use mpitools , only : reduce_maximum +#endif /* MPI */ + use sources , only : update_sources -! local variables are not implicit by default -! implicit none -! local pointers -! type(block_data), pointer :: pdata -! local saved variables -! logical , save :: first = .true. integer , save :: n, n1, n2, n3, n4 - real(kind=8), save :: r, f1, f2 + real(kind=8), save :: r, f1, f2, g1, g2 -! local variables -! - integer :: i + integer :: i, l, m real(kind=8) :: ds real(kind=8) :: tm, dtm + + real(kind=8), dimension(:,:,:,:,:), allocatable :: lerr ! !------------------------------------------------------------------------------- ! #ifdef PROFILE -! start accounting time for one step update -! call start_timer(imu) #endif /* PROFILE */ -! prepare coefficients -! if (first) then n = int(sqrt(1.0d+00 * stages)) @@ -1902,223 +1894,189 @@ module evolution f2 = (1.0d+00 * (n - 1)) / r r = 1.0d+00 * (stages - n) -! update first flag -! + g1 = 1.0d+00 / ( stages - n) - 1.0d+00 / stages + g2 = 1.0d+00 / (2 * stages - n) - 1.0d+00 / stages + first = .false. end if -!= 1st step: U(1) = [1 + dt/r) L] U(1), for i = 1, ..., n1 -! -! calculate the fractional time step -! ds = dt / r -! integrate the intermediate steps + m = get_dblocks() +#if NDIMS == 3 + allocate(lerr(m,nf,nc,nc,nc), stat = status) +#else /* NDIMS == 3 */ + allocate(lerr(m,nf,nc,nc, 1)) +#endif /* NDIMS == 3 */ + lerr(:,:,:,:,:) = 0.0d+00 + +!= 1st step: U(1) = [1 + dt/r) L] U(1), for i = 1, ..., n1 ! do i = 1, n1 -! prepare times -! tm = time + i * ds dtm = ds -! update fluxes -! call update_fluxes() -! assign pdata with the first block on the data block list -! + l = 1 pdata => list_data - -! iterate over all data blocks -! do while (associated(pdata)) -! calculate the variable increment -! call update_increment(pdata) - -! add the source terms -! call update_sources(pdata, tm, dtm, du(:,:,:,:)) -! update the intermediate solution -! pdata%u0(:,:,:,:) = pdata%u0(:,:,:,:) + ds * du(:,:,:,:) -! assign pdata to the next block -! - pdata => pdata%next +#if NDIMS == 3 + lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g1 * du(:,nb:ne,nb:ne,nb:ne) +#else /* NDIMS == 3 */ + lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g1 * du(:,nb:ne,nb:ne, : ) +#endif /* NDIMS == 3 */ + l = l + 1 + pdata => pdata%next end do ! over data blocks -! update primitive variables -! call update_variables(tm, dtm) end do ! n = 1, n1 != 2nd step: U(2) = U(1) ! -! assign pdata with the first block on the data block list -! - pdata => list_data - ! iterate over all data blocks ! + pdata => list_data do while (associated(pdata)) -! update the final solution -! pdata%u1(:,:,:,:) = pdata%u0(:,:,:,:) -! assign pdata to the next block -! pdata => pdata%next end do ! over data blocks != 3rd step: U(1) = [1 + dt/r) L] U(1), for i = n2, ..., n3 -! -! integrate the intermediate steps ! do i = n2, n3 -! prepare times -! tm = time + i * ds dtm = ds -! update fluxes -! call update_fluxes() -! assign pdata with the first block on the data block list -! + l = 1 pdata => list_data - -! iterate over all data blocks -! do while (associated(pdata)) -! calculate the variable increment -! call update_increment(pdata) - -! add the source terms -! call update_sources(pdata, tm, dtm, du(:,:,:,:)) -! update the intermediate solution -! pdata%u0(:,:,:,:) = pdata%u0(:,:,:,:) + ds * du(:,:,:,:) -! assign pdata to the next block -! - pdata => pdata%next +#if NDIMS == 3 + lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g2 * du(:,nb:ne,nb:ne,nb:ne) +#else /* NDIMS == 3 */ + lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g2 * du(:,nb:ne,nb:ne, : ) +#endif /* NDIMS == 3 */ + l = l + 1 + pdata => pdata%next end do ! over data blocks -! update primitive variables -! call update_variables(tm, dtm) end do ! i = n2, n3 != 4th step: U(1) = n * U(2) + (n - 1) * [1 + dt/r L] U(1)) / (2 * n - 1) -! -! prepare times ! tm = time + dt dtm = ds -! update fluxes -! call update_fluxes() -! assign pdata with the first block on the data block list -! - pdata => list_data - ! iterate over all data blocks ! + l = 1 + pdata => list_data do while (associated(pdata)) -! calculate the variable increment -! call update_increment(pdata) - -! add the source terms -! call update_sources(pdata, tm, dtm, du(:,:,:,:)) -! update the final solution -! - pdata%u0(:,:,:,:) = f1 * pdata%u1(:,:,:,:) + f2 * (pdata%u0(:,:,:,:) & + pdata%u0(:,:,:,:) = f1 * pdata%u1(:,:,:,:) + f2 * (pdata%u0(:,:,:,:) & + ds * du(:,:,:,:)) -! update ψ with its source term -! - if (ibp > 0) pdata%u(ibp,:,:,:) = decay * pdata%u(ibp,:,:,:) +#if NDIMS == 3 + lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g2 * du(:,nb:ne,nb:ne,nb:ne) +#else /* NDIMS == 3 */ + lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g2 * du(:,nb:ne,nb:ne, : ) +#endif /* NDIMS == 3 */ + l = l + 1 -! assign pdata to the next block -! pdata => pdata%next - end do ! over data blocks -! update primitive variables -! call update_variables(tm, dtm) -!= the final step: U(1) = [1 + dt/r) L] U(1), for i = n4, ..., m, U(n+1) = U(1) -! -! integrate the intermediate steps +!= final step: U(1) = [1 + dt/r) L] U(1), for i = n4, ..., m, U(n+1) = U(1) ! do i = n4, stages -! prepare times -! tm = time + (i - n) * ds dtm = ds -! update fluxes -! call update_fluxes() -! assign pdata with the first block on the data block list -! + l = 1 pdata => list_data - -! iterate over all data blocks -! do while (associated(pdata)) -! calculate the variable increment -! call update_increment(pdata) - -! add the source terms -! call update_sources(pdata, tm, dtm, du(:,:,:,:)) -! update the intermediate solution -! pdata%u0(:,:,:,:) = pdata%u0(:,:,:,:) + ds * du(:,:,:,:) -! assign pdata to the next block -! - pdata => pdata%next +#if NDIMS == 3 + lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g1 * du(:,nb:ne,nb:ne,nb:ne) +#else /* NDIMS == 3 */ + lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g1 * du(:,nb:ne,nb:ne, : ) +#endif /* NDIMS == 3 */ + l = l + 1 + pdata => pdata%next end do ! over data blocks -! update primitive variables +! update ψ with its source term in the last step ! + if (i == stages .and. ibp > 0) then + pdata => list_data + do while (associated(pdata)) + pdata%u(ibp,:,:,:) = decay * pdata%u(ibp,:,:,:) + pdata => pdata%next + end do ! over data blocks + end if + call update_variables(tm, dtm) end do ! i = n4, stages +! get the maximum error of each variable +! + do l = 1, nf + errors(l) = dt * maxval(abs(lerr(:,l,:,:,:))) + end do + +#ifdef MPI +! reduce errors +! + call reduce_maximum(errors) +#endif /* MPI */ + +! deallocate the local errors array +! + if (allocated(lerr)) deallocate(lerr) + #ifdef PROFILE ! stop accounting time for one step update ! From 8a23da2802a5b687cf6323d0e460bb62c88496f7 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Fri, 28 Aug 2020 22:44:51 -0300 Subject: [PATCH 06/20] EVOLUTION: Remove few comments from evolve_ssprk3_m(). Signed-off-by: Grzegorz Kowal --- sources/evolution.F90 | 6 ------ 1 file changed, 6 deletions(-) diff --git a/sources/evolution.F90 b/sources/evolution.F90 index 5320a43..9a2476d 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -2068,18 +2068,12 @@ module evolution end do #ifdef MPI -! reduce errors -! call reduce_maximum(errors) #endif /* MPI */ -! deallocate the local errors array -! if (allocated(lerr)) deallocate(lerr) #ifdef PROFILE -! stop accounting time for one step update -! call stop_timer(imu) #endif /* PROFILE */ From 73e2743e076de716abc3ed41f4b1aa3263d3035a Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Fri, 28 Aug 2020 23:07:32 -0300 Subject: [PATCH 07/20] DRIVER: Print maximum error. Signed-off-by: Grzegorz Kowal --- sources/driver.F90 | 32 +++++++++++++++++--------------- 1 file changed, 17 insertions(+), 15 deletions(-) diff --git a/sources/driver.F90 b/sources/driver.F90 index 7156863..60b5f8d 100644 --- a/sources/driver.F90 +++ b/sources/driver.F90 @@ -45,7 +45,7 @@ program amun use domains , only : initialize_domains, finalize_domains use equations , only : initialize_equations, finalize_equations use equations , only : print_equations - use equations , only : nv, nf + use equations , only : nv, nf, errors use evolution , only : initialize_evolution, finalize_evolution use evolution , only : print_evolution use evolution , only : advance, new_time_step @@ -118,6 +118,7 @@ program amun real(kind=8) :: zmin = 0.0d+00, zmax = 1.0d+00 real(kind=8) :: tmax = 0.0d+00, trun = 9.999d+03, tsav = 3.0d+01 real(kind=8) :: dtnext = 0.0d+00 + real(kind=8) :: maxerr = 0.0d+00 ! timer indices ! @@ -613,16 +614,16 @@ program amun ! write(*,*) write(*,"(1x,a)" ) "Evolving the system:" - write(*,'(4x,a4,5x,a4,11x,a2,12x,a6,7x,a3)') 'step', 'time', 'dt' & - , 'blocks', 'ETA' + write(*,"(4x,'step',5x,'time',11x,'timestep',7x,'error',5x," // & + "'blocks',7x,'ETA')") #ifdef __INTEL_COMPILER - write(*,'(i8,2(1x,1es14.6),2x,i8,2x,1i4.1,"d",1i2.2,"h",1i2.2,"m"' // & - ',1i2.2,"s",15x,a1,$)') & - step, time, dt, get_nleafs(), ed, eh, em, es, char(13) + write(*,"(i8,2(1x,1es14.6),1x,1es10.2,2x,i8,2x," // & + "1i4.1,'d',1i2.2,'h',1i2.2,'m',1i2.2,'s',15x,a1,$)") & + step, time, dt, maxerr, get_nleafs(), ed, eh, em, es, char(13) #else /* __INTEL_COMPILER */ - write(*,'(i8,2(1x,1es14.6),2x,i8,2x,1i4.1,"d",1i2.2,"h",1i2.2,"m"' // & - ',1i2.2,"s",15x,a1)',advance="no") & - step, time, dt, get_nleafs(), ed, eh, em, es, char(13) + write(*,"(i8,2(1x,1es14.6),1x,1es10.2,2x,i8,2x," // & + "1i4.1,'d',1i2.2,'h',1i2.2,'m',1i2.2,'s',15x,a1)",advance="no") & + step, time, dt, maxerr, get_nleafs(), ed, eh, em, es, char(13) #endif /* __INTEL_COMPILER */ end if @@ -652,6 +653,7 @@ program amun time = time + dt step = step + 1 nsteps = nsteps + 1 + maxerr = maxval(errors) ! get current time in seconds ! @@ -706,13 +708,13 @@ program amun es = es - 60 * em #ifdef __INTEL_COMPILER - write(*,'(i8,2(1x,1es14.6),2x,i8,2x,1i4.1,"d",1i2.2,"h",1i2.2,"m"' // & - ',1i2.2,"s",15x,a1,$)') & - step, time, dt, get_nleafs(), ed, eh, em, es, char(13) + write(*,"(i8,2(1x,1es14.6),1x,1es10.2,2x,i8,2x," // & + "1i4.1,'d',1i2.2,'h',1i2.2,'m',1i2.2,'s',15x,a1,$)") & + step, time, dt, maxerr, get_nleafs(), ed, eh, em, es, char(13) #else /* __INTEL_COMPILER */ - write(*,'(i8,2(1x,1es14.6),2x,i8,2x,1i4.1,"d",1i2.2,"h",1i2.2,"m"' // & - ',1i2.2,"s",15x,a1)',advance="no") & - step, time, dt, get_nleafs(), ed, eh, em, es, char(13) + write(*,"(i8,2(1x,1es14.6),1x,1es10.2,2x,i8,2x," // & + "1i4.1,'d',1i2.2,'h',1i2.2,'m',1i2.2,'s',15x,a1)",advance="no")& + step, time, dt, maxerr, get_nleafs(), ed, eh, em, es, char(13) #endif /* __INTEL_COMPILER */ ! update the timestamp From 60255efe9f5332d8493ee11253e8eaebaf3cd487 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Fri, 28 Aug 2020 23:36:52 -0300 Subject: [PATCH 08/20] EVOLUTION: Add error estimation in evolve_ssprk2_m(). Signed-off-by: Grzegorz Kowal --- sources/evolution.F90 | 143 ++++++++++++++++-------------------------- 1 file changed, 55 insertions(+), 88 deletions(-) diff --git a/sources/evolution.F90 b/sources/evolution.F90 index 9a2476d..8fae884 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -193,7 +193,7 @@ module evolution stages = max(2, min(9, stages)) write(name_int, "('2nd order SSPRK(',i0,',2)')") stages - evolve => evolve_ssprk2 + evolve => evolve_ssprk2_m cfl = (stages - 1) * cfl case ("rk3", "RK3") @@ -880,8 +880,8 @@ module evolution ! !=============================================================================== ! -! subroutine EVOLVE_SSPRK2: -! ------------------------ +! subroutine EVOLVE_SSPRK2_M: +! -------------------------- ! ! Subroutine advances the solution by one time step using the 2nd order ! m-stage Strong Stability Preserving Runge-Kutta time integration method. @@ -897,182 +897,149 @@ module evolution ! !=============================================================================== ! - subroutine evolve_ssprk2() + subroutine evolve_ssprk2_m() -! include external variables -! - use blocks , only : block_data, list_data - use equations, only : ibp - use sources , only : update_sources + use blocks , only : block_data, list_data, get_dblocks + use coordinates, only : nc => ncells, nb, ne + use equations , only : errors, nf, ibp +#ifdef MPI + use mpitools , only : reduce_maximum +#endif /* MPI */ + use sources , only : update_sources -! local variables are not implicit by default -! implicit none -! local pointers -! type(block_data), pointer :: pdata -! local variables -! - integer :: n + integer :: n, l real(kind=8) :: ds real(kind=8) :: tm, dtm -! local saved variables -! logical , save :: first = .true. - real(kind=8), save :: ft, fl, fr + real(kind=8), save :: ft, fl, fr, gt, gl + + real(kind=8), dimension(:,:,:,:,:), allocatable :: lerr ! !------------------------------------------------------------------------------- ! #ifdef PROFILE -! start accounting time for one step update -! call start_timer(imu) #endif /* PROFILE */ -! prepare things which don't change later -! if (first) then -! calculate integration coefficients -! ft = 1.0d+00 / (stages - 1) fl = 1.0d+00 / stages fr = 1.0d+00 - fl -! update first flag -! + gt = fl - ft + gl = fl + first = .false. end if -! calculate the fractional time step -! ds = ft * dt + l = get_dblocks() +#if NDIMS == 3 + allocate(lerr(l,nf,nc,nc,nc), stat = status) +#else /* NDIMS == 3 */ + allocate(lerr(l,nf,nc,nc, 1)) +#endif /* NDIMS == 3 */ + lerr(:,:,:,:,:) = 0.0d+00 + != 1st step: U(0) = U(n) -! -! assign pdata with the first block on the data block list ! pdata => list_data - -! iterate over all data blocks -! do while (associated(pdata)) -! copy conservative array u0 to u1 -! pdata%u1(:,:,:,:) = pdata%u0(:,:,:,:) -! update the conservative variable pointer -! pdata%u => pdata%u1 -! assign pdata to the next block -! pdata => pdata%next - end do ! over data blocks != 2nd step: U(i) = [1 + dt/(m-1) L] U(i-1), for i = 1, ..., m-1 -! -! integrate the intermediate steps ! do n = 1, stages - 1 -! prepare times -! tm = time + n * ds dtm = ds -! update fluxes -! call update_fluxes() -! assign pdata with the first block on the data block list -! + l = 1 pdata => list_data - -! iterate over all data blocks -! do while (associated(pdata)) -! calculate the variable increment -! call update_increment(pdata) - -! add the source terms -! call update_sources(pdata, tm, dtm, du(:,:,:,:)) -! update the intermediate solution -! pdata%u1(:,:,:,:) = pdata%u1(:,:,:,:) + ds * du(:,:,:,:) -! assign pdata to the next block -! - pdata => pdata%next +#if NDIMS == 3 + lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + gt * du(:,nb:ne,nb:ne,nb:ne) +#else /* NDIMS == 3 */ + lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + gt * du(:,nb:ne,nb:ne, : ) +#endif /* NDIMS == 3 */ + l = l + 1 + pdata => pdata%next end do ! over data blocks -! update primitive variables -! call update_variables(tm, dtm) end do ! n = 1, stages - 1 != the final step: U(n+1) = 1/m U(0) + (m-1)/m [1 + dt/(m-1) L] U(m-1) -! -! prepare times ! tm = time + dt dtm = fl * dt -! update fluxes -! call update_fluxes() -! assign pdata with the first block on the data block list -! + l = 1 pdata => list_data - -! iterate over all data blocks -! do while (associated(pdata)) -! calculate the variable increment -! call update_increment(pdata) - -! add the source terms -! call update_sources(pdata, tm, dtm, du(:,:,:,:)) -! update the final solution -! pdata%u0(:,:,:,:) = fl * pdata%u0(:,:,:,:) & + fr * (pdata%u1(:,:,:,:) + ds * du(:,:,:,:)) -! update the conservative variable pointer -! pdata%u => pdata%u0 +#if NDIMS == 3 + lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + gl * du(:,nb:ne,nb:ne,nb:ne) +#else /* NDIMS == 3 */ + lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + gl * du(:,nb:ne,nb:ne, : ) +#endif /* NDIMS == 3 */ + l = l + 1 + ! update ψ with its source term ! if (ibp > 0) pdata%u(ibp,:,:,:) = decay * pdata%u(ibp,:,:,:) -! assign pointer to the next block -! pdata => pdata%next - end do ! over data blocks -! update primitive variables -! call update_variables(tm, dtm) +! get the maximum error of each variable +! + do l = 1, nf + errors(l) = dt * maxval(abs(lerr(:,l,:,:,:))) + end do + +#ifdef MPI + call reduce_maximum(errors) +#endif /* MPI */ + + if (allocated(lerr)) deallocate(lerr) + #ifdef PROFILE ! stop accounting time for one step update ! @@ -1081,7 +1048,7 @@ module evolution !------------------------------------------------------------------------------- ! - end subroutine evolve_ssprk2 + end subroutine evolve_ssprk2_m ! !=============================================================================== ! From a754fd6a563d72246ce97bdf8b181fae12a559aa Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Fri, 28 Aug 2020 23:38:12 -0300 Subject: [PATCH 09/20] EVOLUTION: Remove one variable from evolve_ssprk3_m(). Signed-off-by: Grzegorz Kowal --- sources/evolution.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/sources/evolution.F90 b/sources/evolution.F90 index 8fae884..ef5c777 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -1836,7 +1836,7 @@ module evolution integer , save :: n, n1, n2, n3, n4 real(kind=8), save :: r, f1, f2, g1, g2 - integer :: i, l, m + integer :: i, l real(kind=8) :: ds real(kind=8) :: tm, dtm @@ -1870,11 +1870,11 @@ module evolution ds = dt / r - m = get_dblocks() + l = get_dblocks() #if NDIMS == 3 - allocate(lerr(m,nf,nc,nc,nc), stat = status) + allocate(lerr(l,nf,nc,nc,nc), stat = status) #else /* NDIMS == 3 */ - allocate(lerr(m,nf,nc,nc, 1)) + allocate(lerr(l,nf,nc,nc, 1)) #endif /* NDIMS == 3 */ lerr(:,:,:,:,:) = 0.0d+00 From d86fea463c2dea6894ff0cdf9fee55e4a18e8d04 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Mon, 31 Aug 2020 11:52:19 -0300 Subject: [PATCH 10/20] EVOLUTION: Remove stat atribute from allocation of lerr(..). Signed-off-by: Grzegorz Kowal --- sources/evolution.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/sources/evolution.F90 b/sources/evolution.F90 index ef5c777..0d93472 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -943,7 +943,7 @@ module evolution l = get_dblocks() #if NDIMS == 3 - allocate(lerr(l,nf,nc,nc,nc), stat = status) + allocate(lerr(l,nf,nc,nc,nc)) #else /* NDIMS == 3 */ allocate(lerr(l,nf,nc,nc, 1)) #endif /* NDIMS == 3 */ @@ -1872,7 +1872,7 @@ module evolution l = get_dblocks() #if NDIMS == 3 - allocate(lerr(l,nf,nc,nc,nc), stat = status) + allocate(lerr(l,nf,nc,nc,nc)) #else /* NDIMS == 3 */ allocate(lerr(l,nf,nc,nc, 1)) #endif /* NDIMS == 3 */ From f24f071540426b2c2de27bcc47d199c92c71b438 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Mon, 31 Aug 2020 18:52:21 -0300 Subject: [PATCH 11/20] BLOCKS: Add additional conserved variable array to data blocks. It is required for high-order low-storage embedded SSPRK integration methods. Signed-off-by: Grzegorz Kowal --- sources/blocks.F90 | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/sources/blocks.F90 b/sources/blocks.F90 index 0fc32b6..12fc0ce 100644 --- a/sources/blocks.F90 +++ b/sources/blocks.F90 @@ -260,10 +260,11 @@ module blocks real(kind=8), dimension(:,:,:,:) , pointer :: u ! an allocatable arrays to store all conserved - ! variables (required two for Runge-Kutta - ! temporal integration methods) + ! variables (requires two additional arrays + ! for embedded low-storage explicit optimal + ! high-order SSPRK integration methods) ! - real(kind=8), dimension(:,:,:,:) , allocatable :: u0, u1 + real(kind=8), dimension(:,:,:,:) , allocatable :: u0, u1, u2 ! an allocatable array to store all primitive ! variables @@ -1178,8 +1179,8 @@ module blocks ! allocate space for conserved and primitive variables, and fluxes ! allocate(pdata%u0(nvars,nx,ny,nz), pdata%u1(nvars,nx,ny,nz), & - pdata%q (nvars,nx,ny,nz), pdata%f(ndims,nflux,nx,ny,nz), & - stat = status) + pdata%u2(nvars,nx,ny,nz), pdata%q (nvars,nx,ny,nz), & + pdata%f(ndims,nflux,nx,ny,nz), stat = status) ! check if allocation succeeded ! @@ -1189,6 +1190,7 @@ module blocks ! pdata%u0 = 0.0d+00 pdata%u1 = 0.0d+00 + pdata%u2 = 0.0d+00 pdata%q = 0.0d+00 pdata%f = 0.0d+00 @@ -1282,6 +1284,7 @@ module blocks ! if (allocated(pdata%u0)) deallocate(pdata%u0, stat = status) if (allocated(pdata%u1)) deallocate(pdata%u1, stat = status) + if (allocated(pdata%u2)) deallocate(pdata%u2, stat = status) ! deallocate primitive variables ! From 931c984dca06e7ab67130433a9d667755bcfaa02 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Mon, 31 Aug 2020 20:10:48 -0300 Subject: [PATCH 12/20] EVOLUTION: Rewrite evolve_ssprk2_m() to use %u2(:,:,:,:). Signed-off-by: Grzegorz Kowal --- sources/evolution.F90 | 35 ++++++++++++++++++++++------------- 1 file changed, 22 insertions(+), 13 deletions(-) diff --git a/sources/evolution.F90 b/sources/evolution.F90 index 0d93472..885b1df 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -949,19 +949,20 @@ module evolution #endif /* NDIMS == 3 */ lerr(:,:,:,:,:) = 0.0d+00 -!= 1st step: U(0) = U(n) +!= 1st step: U(1) = U(n), U(2) = U(n) ! pdata => list_data do while (associated(pdata)) pdata%u1(:,:,:,:) = pdata%u0(:,:,:,:) + pdata%u2(:,:,:,:) = pdata%u0(:,:,:,:) pdata%u => pdata%u1 pdata => pdata%next end do ! over data blocks -!= 2nd step: U(i) = [1 + dt/(m-1) L] U(i-1), for i = 1, ..., m-1 +!= 2nd step: U(1) = U(1) + dt/(m-1) F[U(1)], for i = 1, ..., m-1 ! do n = 1, stages - 1 @@ -993,7 +994,7 @@ module evolution end do ! n = 1, stages - 1 -!= the final step: U(n+1) = 1/m U(0) + (m-1)/m [1 + dt/(m-1) L] U(m-1) +!= the final step: U(1) = (m-1)/m U(1) + 1/m U(2) + dt/m F[U(1)] ! tm = time + dt dtm = fl * dt @@ -1007,10 +1008,8 @@ module evolution call update_increment(pdata) call update_sources(pdata, tm, dtm, du(:,:,:,:)) - pdata%u0(:,:,:,:) = fl * pdata%u0(:,:,:,:) & - + fr * (pdata%u1(:,:,:,:) + ds * du(:,:,:,:)) - - pdata%u => pdata%u0 + pdata%u1(:,:,:,:) = fr * pdata%u1(:,:,:,:) & + + fl * (pdata%u2(:,:,:,:) + dt * du(:,:,:,:)) #if NDIMS == 3 lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + gl * du(:,nb:ne,nb:ne,nb:ne) @@ -1019,15 +1018,9 @@ module evolution #endif /* NDIMS == 3 */ l = l + 1 -! update ψ with its source term -! - if (ibp > 0) pdata%u(ibp,:,:,:) = decay * pdata%u(ibp,:,:,:) - pdata => pdata%next end do ! over data blocks - call update_variables(tm, dtm) - ! get the maximum error of each variable ! do l = 1, nf @@ -1040,6 +1033,22 @@ module evolution if (allocated(lerr)) deallocate(lerr) + pdata => list_data + do while (associated(pdata)) + + pdata%u0(:,:,:,:) = pdata%u1(:,:,:,:) + + pdata%u => pdata%u0 + +! update ψ with its source term +! + if (ibp > 0) pdata%u(ibp,:,:,:) = decay * pdata%u(ibp,:,:,:) + + pdata => pdata%next + end do ! over data blocks + + call update_variables(tm, dtm) + #ifdef PROFILE ! stop accounting time for one step update ! From 1d0cf0758c9d3b26d8eecd896ed9d6fb5e3aa073 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Mon, 31 Aug 2020 22:27:38 -0300 Subject: [PATCH 13/20] EVOLUTION: Rewrite evolve_ssprk3_m() to use %u2(:,:,:,:). Signed-off-by: Grzegorz Kowal --- sources/evolution.F90 | 74 ++++++++++++++++++++++++++++++------------- 1 file changed, 52 insertions(+), 22 deletions(-) diff --git a/sources/evolution.F90 b/sources/evolution.F90 index 885b1df..8371816 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -1887,7 +1887,19 @@ module evolution #endif /* NDIMS == 3 */ lerr(:,:,:,:,:) = 0.0d+00 -!= 1st step: U(1) = [1 + dt/r) L] U(1), for i = 1, ..., n1 +!= 1st step: U(1) = U(n) +! + pdata => list_data + do while (associated(pdata)) + + pdata%u1(:,:,:,:) = pdata%u0(:,:,:,:) + + pdata%u => pdata%u1 + + pdata => pdata%next + end do ! over data blocks + +!= 2ns step: U(1) = U(1) + dt/r F[U(1)], for i = 1, ..., (n-1)*(n-2)/2 ! do i = 1, n1 @@ -1903,7 +1915,7 @@ module evolution call update_increment(pdata) call update_sources(pdata, tm, dtm, du(:,:,:,:)) - pdata%u0(:,:,:,:) = pdata%u0(:,:,:,:) + ds * du(:,:,:,:) + pdata%u1(:,:,:,:) = pdata%u1(:,:,:,:) + ds * du(:,:,:,:) #if NDIMS == 3 lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g1 * du(:,nb:ne,nb:ne,nb:ne) @@ -1919,20 +1931,20 @@ module evolution end do ! n = 1, n1 -!= 2nd step: U(2) = U(1) +!= 3rd step: U(2) = U(1) ! ! iterate over all data blocks ! pdata => list_data do while (associated(pdata)) - pdata%u1(:,:,:,:) = pdata%u0(:,:,:,:) + pdata%u2(:,:,:,:) = pdata%u1(:,:,:,:) pdata => pdata%next end do ! over data blocks -!= 3rd step: U(1) = [1 + dt/r) L] U(1), for i = n2, ..., n3 +!= 4th step: U(1) = U(1) + dt/r F[U(1)], for i = (n-1)*(n-2)/2+1, ..., n*(n+1)/2-1 ! do i = n2, n3 @@ -1948,7 +1960,7 @@ module evolution call update_increment(pdata) call update_sources(pdata, tm, dtm, du(:,:,:,:)) - pdata%u0(:,:,:,:) = pdata%u0(:,:,:,:) + ds * du(:,:,:,:) + pdata%u1(:,:,:,:) = pdata%u1(:,:,:,:) + ds * du(:,:,:,:) #if NDIMS == 3 lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g2 * du(:,nb:ne,nb:ne,nb:ne) @@ -1964,7 +1976,7 @@ module evolution end do ! i = n2, n3 -!= 4th step: U(1) = n * U(2) + (n - 1) * [1 + dt/r L] U(1)) / (2 * n - 1) +!= 5th step: U(1) = n * U(2) + (n - 1) * ( U(1) + dt/r F[U(1)] ) / (2 * n - 1) ! tm = time + dt dtm = ds @@ -1980,8 +1992,8 @@ module evolution call update_increment(pdata) call update_sources(pdata, tm, dtm, du(:,:,:,:)) - pdata%u0(:,:,:,:) = f1 * pdata%u1(:,:,:,:) + f2 * (pdata%u0(:,:,:,:) & - + ds * du(:,:,:,:)) + pdata%u1(:,:,:,:) = f1 * pdata%u2(:,:,:,:) & + + f2 * (pdata%u1(:,:,:,:) + ds * du(:,:,:,:)) #if NDIMS == 3 lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g2 * du(:,nb:ne,nb:ne,nb:ne) @@ -1995,7 +2007,7 @@ module evolution call update_variables(tm, dtm) -!= final step: U(1) = [1 + dt/r) L] U(1), for i = n4, ..., m, U(n+1) = U(1) +!= 6th step: U(1) = U(1) + dt/r F[U(1)], for i = n*(n+1)/2, ..., m ! do i = n4, stages @@ -2011,7 +2023,7 @@ module evolution call update_increment(pdata) call update_sources(pdata, tm, dtm, du(:,:,:,:)) - pdata%u0(:,:,:,:) = pdata%u0(:,:,:,:) + ds * du(:,:,:,:) + pdata%u1(:,:,:,:) = pdata%u1(:,:,:,:) + ds * du(:,:,:,:) #if NDIMS == 3 lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g1 * du(:,nb:ne,nb:ne,nb:ne) @@ -2023,17 +2035,7 @@ module evolution pdata => pdata%next end do ! over data blocks -! update ψ with its source term in the last step -! - if (i == stages .and. ibp > 0) then - pdata => list_data - do while (associated(pdata)) - pdata%u(ibp,:,:,:) = decay * pdata%u(ibp,:,:,:) - pdata => pdata%next - end do ! over data blocks - end if - - call update_variables(tm, dtm) + if (i < stages) call update_variables(tm, dtm) end do ! i = n4, stages @@ -2049,6 +2051,34 @@ module evolution if (allocated(lerr)) deallocate(lerr) +!= final step: U(n+1) = U(1) +! + tm = time + dt + dtm = ds + + pdata => list_data + do while (associated(pdata)) + + pdata%u0(:,:,:,:) = pdata%u1(:,:,:,:) + + pdata%u => pdata%u0 + + pdata => pdata%next + end do ! over data blocks + + call update_variables(tm, dtm) + +! update ψ with its source term in the last step +! + if (i == stages .and. ibp > 0) then + pdata => list_data + do while (associated(pdata)) + pdata%u(ibp,:,:,:) = decay * pdata%u(ibp,:,:,:) + pdata => pdata%next + end do ! over data blocks + end if + + #ifdef PROFILE call stop_timer(imu) #endif /* PROFILE */ From 0c04d4986b15382f05398bc11bd2a674b5d8ce49 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Mon, 31 Aug 2020 22:29:40 -0300 Subject: [PATCH 14/20] EVOLUTION: Add parameters for absolute and relative tolerances. Signed-off-by: Grzegorz Kowal --- sources/evolution.F90 | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/sources/evolution.F90 b/sources/evolution.F90 index 8371816..a61e777 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -77,6 +77,11 @@ module evolution real(kind=8) , save :: dtn = 1.0d+00 real(kind=8) , save :: dte = 0.0d+00 +! absolute and relative tolerances +! + real(kind=8) , save :: atol = 1.0d-04 + real(kind=8) , save :: rtol = 1.0d-04 + ! the increment array ! real(kind=8), dimension(:,:,:,:), allocatable :: du @@ -174,6 +179,11 @@ module evolution call get_parameter("cfl" , cfl ) call get_parameter("glm_alpha" , glm_alpha ) +! get absolute and relative tolerances +! + call get_parameter("atol", atol) + call get_parameter("rtol", rtol) + ! select the integration method, check the correctness of the integration ! parameters and adjust the CFL coefficient if necessary ! From c1fa4d1cffcd62338e85c6663e5ee8549c0cf98d Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Mon, 31 Aug 2020 22:38:46 -0300 Subject: [PATCH 15/20] EVOLUTION: Estimate timestep due to error limitation. Signed-off-by: Grzegorz Kowal --- sources/evolution.F90 | 43 +++++++++++++++++++++++++++++++++++++++---- 1 file changed, 39 insertions(+), 4 deletions(-) diff --git a/sources/evolution.F90 b/sources/evolution.F90 index a61e777..efb2485 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -922,8 +922,7 @@ module evolution type(block_data), pointer :: pdata integer :: n, l - real(kind=8) :: ds - real(kind=8) :: tm, dtm + real(kind=8) :: tm, dtm, ds, umax, utol logical , save :: first = .true. real(kind=8), save :: ft, fl, fr, gt, gl @@ -1031,6 +1030,16 @@ module evolution pdata => pdata%next end do ! over data blocks +! find umax +! + umax = 0.0d+00 + pdata => list_data + do while (associated(pdata)) + umax = max(umax, maxval(abs(pdata%u0(:,:,:,:))), & + maxval(abs(pdata%u1(:,:,:,:)))) + pdata => pdata%next + end do ! over data blocks + ! get the maximum error of each variable ! do l = 1, nf @@ -1038,11 +1047,22 @@ module evolution end do #ifdef MPI + call reduce_maximum(umax) call reduce_maximum(errors) #endif /* MPI */ +! calculate tolerance and time step +! + utol = atol + rtol * umax + dte = dt * (maxval(errors) / utol)**(-2.5d-01) + if (allocated(lerr)) deallocate(lerr) +!= final step: U(n+1) = U(1) +! + tm = time + dt + dtm = ds + pdata => list_data do while (associated(pdata)) @@ -1856,8 +1876,7 @@ module evolution real(kind=8), save :: r, f1, f2, g1, g2 integer :: i, l - real(kind=8) :: ds - real(kind=8) :: tm, dtm + real(kind=8) :: tm, dtm, ds, umax, utol real(kind=8), dimension(:,:,:,:,:), allocatable :: lerr ! @@ -2049,6 +2068,16 @@ module evolution end do ! i = n4, stages +! find umax +! + umax = 0.0d+00 + pdata => list_data + do while (associated(pdata)) + umax = max(umax, maxval(abs(pdata%u0(:,:,:,:))), & + maxval(abs(pdata%u1(:,:,:,:)))) + pdata => pdata%next + end do ! over data blocks + ! get the maximum error of each variable ! do l = 1, nf @@ -2056,9 +2085,15 @@ module evolution end do #ifdef MPI + call reduce_maximum(umax) call reduce_maximum(errors) #endif /* MPI */ +! calculate tolerance and time step +! + utol = atol + rtol * umax + dte = dt * (maxval(errors) / utol)**(-2.5d-01) + if (allocated(lerr)) deallocate(lerr) != final step: U(n+1) = U(1) From f199e231e6674a6ab9ce1d9e39adfbbea651097b Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 1 Sep 2020 16:38:11 -0300 Subject: [PATCH 16/20] EVOLUTION, IO: Add tolerances, interations and rejections counters. Store these variables in restart snapshots. Signed-off-by: Grzegorz Kowal --- sources/driver.F90 | 2 +- sources/evolution.F90 | 50 ++++++++++++++++++++++++++++++++++--------- sources/io.F90 | 39 +++++++++++++++++++++++++++++++-- 3 files changed, 78 insertions(+), 13 deletions(-) diff --git a/sources/driver.F90 b/sources/driver.F90 index 60b5f8d..59ffd15 100644 --- a/sources/driver.F90 +++ b/sources/driver.F90 @@ -773,7 +773,7 @@ program amun "Problem finalizing FORCING module!" end if 2100 continue - call finalize_evolution(status) + call finalize_evolution(master, status) if (check_status(status /= 0) .and. master) then write(error_unit,"('[AMUN::program]: ', a)") & "Problem finalizing EVOLUTION module!" diff --git a/sources/evolution.F90 b/sources/evolution.F90 index efb2485..3a9cde4 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -77,10 +77,19 @@ module evolution real(kind=8) , save :: dtn = 1.0d+00 real(kind=8) , save :: dte = 0.0d+00 -! absolute and relative tolerances +! the absolute and relative tolerances, the maximum number of passes for +! the adaptive step, the total number of integration iterations, +! the number of rejected iterations ! - real(kind=8) , save :: atol = 1.0d-04 - real(kind=8) , save :: rtol = 1.0d-04 + real(kind=8) , save :: atol = 1.0d-04 + real(kind=8) , save :: rtol = 1.0d-04 + integer , save :: mrej = 5 + integer , save :: niterations = 0 + integer , save :: nrejections = 0 + +! errors in three recent steps +! + real(kind=8), dimension(3), save :: errs ! the increment array ! @@ -98,6 +107,7 @@ module evolution ! declare public variables ! public :: step, time, dt, dtn, dte, cfl, glm_alpha + public :: atol, rtol, mrej, niterations, nrejections, errs !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! @@ -179,10 +189,12 @@ module evolution call get_parameter("cfl" , cfl ) call get_parameter("glm_alpha" , glm_alpha ) -! get absolute and relative tolerances +! get the absolute and relative tolerances, the maximum number of passes for +! the adaptive step ! - call get_parameter("atol", atol) - call get_parameter("rtol", rtol) + call get_parameter("absolute_tolerance", atol) + call get_parameter("relative_tolerance", rtol) + call get_parameter("maximum_rejections", mrej) ! select the integration method, check the correctness of the integration ! parameters and adjust the CFL coefficient if necessary @@ -281,6 +293,10 @@ module evolution end if ! status +! reset recent error history +! + errs = 1.0d+00 + #ifdef PROFILE ! stop accounting time for module initialization/finalization ! @@ -304,10 +320,11 @@ module evolution ! !=============================================================================== ! - subroutine finalize_evolution(status) + subroutine finalize_evolution(verbose, status) ! include external procedures ! + use helpers , only : print_section, print_parameter use iso_fortran_env, only : error_unit ! local variables are not implicit by default @@ -316,6 +333,7 @@ module evolution ! subroutine arguments ! + logical, intent(in) :: verbose integer, intent(out) :: status ! local parameters @@ -349,6 +367,14 @@ module evolution end if end if +! print integration summary +! + if (niterations > 0) then + call print_section(verbose, "Integration summary") + call print_parameter(verbose, "Number of iterations", niterations) + call print_parameter(verbose, "Number of rejections", nrejections) + end if + #ifdef PROFILE ! stop accounting time for module initialization/finalization ! @@ -392,11 +418,14 @@ module evolution if (verbose) then call print_section(verbose, "Evolution") - call print_parameter(verbose, "time advance" , name_int) - call print_parameter(verbose, "CFL coefficient" , cfl ) + call print_parameter(verbose, "time advance", name_int) + call print_parameter(verbose, "CFL coefficient", cfl) if (magnetized) then call print_parameter(verbose, "GLM alpha coefficient", glm_alpha) end if + call print_parameter(verbose, "absolute tolerance", atol) + call print_parameter(verbose, "relative tolerance", rtol) + call print_parameter(verbose, "maximum rejections", mrej) end if @@ -620,7 +649,8 @@ module evolution ! round the time ! - if (dtnext > 0.0d+00) dt = min(dtn, dtnext) + if (dte <= 0.0d+00) dte = dtn + if (dtnext > 0.0d+00) dt = min(dtn, dte, dtnext) #ifdef PROFILE ! stop accounting time for new time step estimation diff --git a/sources/io.F90 b/sources/io.F90 index bba455f..0e4b900 100644 --- a/sources/io.F90 +++ b/sources/io.F90 @@ -1261,6 +1261,7 @@ module io use coordinates , only : nn => bcells, ncells use coordinates , only : xmin, xmax, ymin, ymax, zmin, zmax use evolution , only : step, time, dt, dtn + use evolution , only : niterations, nrejections, errs use forcing , only : nmodes, fcoefs, einj use hash , only : xxh64 use iso_fortran_env, only : error_unit @@ -1412,6 +1413,16 @@ module io read(svalue, fmt = *) dt case('dtn') read(svalue, fmt = *) dtn + case('niterations') + read(svalue, fmt = *) niterations + case('nrejections') + read(svalue, fmt = *) nrejections + case('errs(1)') + read(svalue, fmt = *) errs(1) + case('errs(2)') + read(svalue, fmt = *) errs(2) + case('errs(3)') + read(svalue, fmt = *) errs(3) case('fields') il = index(line, 'digest="') + 8 iu = index(line(il:), '"') + il - 2 @@ -2065,7 +2076,8 @@ module io #endif /* NDIMS == 3 */ use coordinates , only : bdims => domain_base_dims use equations , only : eqsys, eos, nv - use evolution , only : step, time, dt, dtn, cfl, glm_alpha + use evolution , only : step, time, dt, dtn, cfl, glm_alpha, errs + use evolution , only : atol, rtol, mrej, niterations, nrejections use forcing , only : nmodes, fcoefs, einj use iso_fortran_env, only : error_unit use mpitools , only : nprocs, nproc @@ -2216,6 +2228,14 @@ module io call write_attribute_xml(lun, "dtn" , dtn) call write_attribute_xml(lun, "cfl" , cfl) call write_attribute_xml(lun, "glm_alpha", glm_alpha) + call write_attribute_xml(lun, "absolute_tolerance", atol) + call write_attribute_xml(lun, "relative_tolerance", rtol) + call write_attribute_xml(lun, "maximum_rejections", mrej) + call write_attribute_xml(lun, "niterations", niterations) + call write_attribute_xml(lun, "nrejections", nrejections) + call write_attribute_xml(lun, "errs(1)", errs(1)) + call write_attribute_xml(lun, "errs(2)", errs(2)) + call write_attribute_xml(lun, "errs(3)", errs(3)) write(lun,"(a)") '' write(lun,"(a)") '' call write_attribute_xml(lun, "nmodes" , nmodes) @@ -3990,7 +4010,8 @@ module io use coordinates , only : xmin, xmax, ymin, ymax, zmin, zmax use coordinates , only : periodic use equations , only : eqsys, eos, adiabatic_index, csnd - use evolution , only : step, time, dt, dtn, cfl, glm_alpha + use evolution , only : step, time, dt, dtn, cfl, glm_alpha, errs + use evolution , only : atol, rtol, mrej, niterations, nrejections use forcing , only : nmodes, einj, fcoefs use hdf5 , only : hid_t use hdf5 , only : h5gcreate_f, h5gclose_f @@ -4078,6 +4099,9 @@ module io call write_attribute(gid, 'step' , step ) call write_attribute(gid, 'isnap' , isnap ) call write_attribute(gid, 'periodic', per(:) ) + call write_attribute(gid, 'niterations', niterations) + call write_attribute(gid, 'nrejections', nrejections) + call write_attribute(gid, 'maximum_rejections', mrej) ! store the real attributes ! @@ -4098,6 +4122,11 @@ module io if (eos == 'iso') then call write_attribute(gid, 'sound_speed', csnd) end if + call write_attribute(gid, 'absolute_tolerance', atol) + call write_attribute(gid, 'relative_tolerance', rtol) + call write_attribute(gid, 'errs(1)', errs(1)) + call write_attribute(gid, 'errs(2)', errs(2)) + call write_attribute(gid, 'errs(3)', errs(3)) ! store the vector attributes ! @@ -4182,6 +4211,7 @@ module io use coordinates , only : ncells use coordinates , only : xmin, xmax, ymin, ymax, zmin, zmax use evolution , only : step, time, dt, dtn + use evolution , only : niterations, nrejections, errs use forcing , only : nmodes, fcoefs use hdf5 , only : hid_t use hdf5 , only : h5gopen_f, h5gclose_f @@ -4242,6 +4272,8 @@ module io call read_attribute(gid, 'nseeds' , lnseeds ) call read_attribute(gid, 'step' , step ) call read_attribute(gid, 'isnap' , isnap ) + call read_attribute(gid, 'niterations', niterations) + call read_attribute(gid, 'nrejections', nrejections) ! restore double precision attributes ! @@ -4254,6 +4286,9 @@ module io call read_attribute(gid, 'time', time) call read_attribute(gid, 'dt' , dt ) call read_attribute(gid, 'dtn' , dtn ) + call read_attribute(gid, 'errs(1)', errs(1)) + call read_attribute(gid, 'errs(2)', errs(2)) + call read_attribute(gid, 'errs(3)', errs(3)) ! check the number of dimensions ! From d9696085410f2fd8c9c0958317ca2ecaf9d4a576 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 1 Sep 2020 16:48:28 -0300 Subject: [PATCH 17/20] EVOLUTION: Implement step repetition in evolve_ssprk3_m(). If the desired tolerance was not achieved, the integration step is repeated. The tolerance is controlled by two parameters: "absolute_tolerance" and "relative_tolerance". The number of allowed repetitions is controlled by parameter "maximum_repetitions". Signed-off-by: Grzegorz Kowal --- sources/evolution.F90 | 249 ++++++++++++++++++++++++------------------ 1 file changed, 140 insertions(+), 109 deletions(-) diff --git a/sources/evolution.F90 b/sources/evolution.F90 index 3a9cde4..6fec4de 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -1905,8 +1905,9 @@ module evolution integer , save :: n, n1, n2, n3, n4 real(kind=8), save :: r, f1, f2, g1, g2 - integer :: i, l - real(kind=8) :: tm, dtm, ds, umax, utol + logical :: test + integer :: i, l, nrej + real(kind=8) :: tm, dtm, ds, umax, utol, emax real(kind=8), dimension(:,:,:,:,:), allocatable :: lerr ! @@ -1936,82 +1937,120 @@ module evolution end if - ds = dt / r - l = get_dblocks() #if NDIMS == 3 allocate(lerr(l,nf,nc,nc,nc)) #else /* NDIMS == 3 */ allocate(lerr(l,nf,nc,nc, 1)) #endif /* NDIMS == 3 */ - lerr(:,:,:,:,:) = 0.0d+00 + + test = .true. + nrej = 0 + + do while(test) + + lerr(:,:,:,:,:) = 0.0d+00 + + ds = dt / r != 1st step: U(1) = U(n) ! - pdata => list_data - do while (associated(pdata)) - - pdata%u1(:,:,:,:) = pdata%u0(:,:,:,:) - - pdata%u => pdata%u1 - - pdata => pdata%next - end do ! over data blocks - -!= 2ns step: U(1) = U(1) + dt/r F[U(1)], for i = 1, ..., (n-1)*(n-2)/2 -! - do i = 1, n1 - - tm = time + i * ds - dtm = ds - - call update_fluxes() - - l = 1 pdata => list_data do while (associated(pdata)) - call update_increment(pdata) - call update_sources(pdata, tm, dtm, du(:,:,:,:)) + pdata%u1(:,:,:,:) = pdata%u0(:,:,:,:) - pdata%u1(:,:,:,:) = pdata%u1(:,:,:,:) + ds * du(:,:,:,:) - -#if NDIMS == 3 - lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g1 * du(:,nb:ne,nb:ne,nb:ne) -#else /* NDIMS == 3 */ - lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g1 * du(:,nb:ne,nb:ne, : ) -#endif /* NDIMS == 3 */ - l = l + 1 + pdata%u => pdata%u1 pdata => pdata%next end do ! over data blocks - call update_variables(tm, dtm) +!= 2ns step: U(1) = U(1) + dt/r F[U(1)], for i = 1, ..., (n-1)*(n-2)/2 +! + do i = 1, n1 - end do ! n = 1, n1 + tm = time + i * ds + dtm = ds + + call update_fluxes() + + l = 1 + pdata => list_data + do while (associated(pdata)) + + call update_increment(pdata) + call update_sources(pdata, tm, dtm, du(:,:,:,:)) + + pdata%u1(:,:,:,:) = pdata%u1(:,:,:,:) + ds * du(:,:,:,:) + +#if NDIMS == 3 + lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g1 * du(:,nb:ne,nb:ne,nb:ne) +#else /* NDIMS == 3 */ + lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g1 * du(:,nb:ne,nb:ne, : ) +#endif /* NDIMS == 3 */ + l = l + 1 + + pdata => pdata%next + end do ! over data blocks + + call update_variables(tm, dtm) + + end do ! n = 1, n1 != 3rd step: U(2) = U(1) ! ! iterate over all data blocks ! - pdata => list_data - do while (associated(pdata)) + pdata => list_data + do while (associated(pdata)) - pdata%u2(:,:,:,:) = pdata%u1(:,:,:,:) + pdata%u2(:,:,:,:) = pdata%u1(:,:,:,:) - pdata => pdata%next + pdata => pdata%next - end do ! over data blocks + end do ! over data blocks != 4th step: U(1) = U(1) + dt/r F[U(1)], for i = (n-1)*(n-2)/2+1, ..., n*(n+1)/2-1 ! - do i = n2, n3 + do i = n2, n3 - tm = time + i * ds + tm = time + i * ds + dtm = ds + + call update_fluxes() + + l = 1 + pdata => list_data + do while (associated(pdata)) + + call update_increment(pdata) + call update_sources(pdata, tm, dtm, du(:,:,:,:)) + + pdata%u1(:,:,:,:) = pdata%u1(:,:,:,:) + ds * du(:,:,:,:) + +#if NDIMS == 3 + lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g2 * du(:,nb:ne,nb:ne,nb:ne) +#else /* NDIMS == 3 */ + lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g2 * du(:,nb:ne,nb:ne, : ) +#endif /* NDIMS == 3 */ + l = l + 1 + + pdata => pdata%next + end do ! over data blocks + + call update_variables(tm, dtm) + + end do ! i = n2, n3 + +!= 5th step: U(1) = n * U(2) + (n - 1) * ( U(1) + dt/r F[U(1)] ) / (2 * n - 1) +! + tm = time + dt dtm = ds call update_fluxes() +! iterate over all data blocks +! l = 1 pdata => list_data do while (associated(pdata)) @@ -2019,7 +2058,8 @@ module evolution call update_increment(pdata) call update_sources(pdata, tm, dtm, du(:,:,:,:)) - pdata%u1(:,:,:,:) = pdata%u1(:,:,:,:) + ds * du(:,:,:,:) + pdata%u1(:,:,:,:) = f1 * pdata%u2(:,:,:,:) & + + f2 * (pdata%u1(:,:,:,:) + ds * du(:,:,:,:)) #if NDIMS == 3 lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g2 * du(:,nb:ne,nb:ne,nb:ne) @@ -2033,96 +2073,87 @@ module evolution call update_variables(tm, dtm) - end do ! i = n2, n3 - -!= 5th step: U(1) = n * U(2) + (n - 1) * ( U(1) + dt/r F[U(1)] ) / (2 * n - 1) -! - tm = time + dt - dtm = ds - - call update_fluxes() - -! iterate over all data blocks -! - l = 1 - pdata => list_data - do while (associated(pdata)) - - call update_increment(pdata) - call update_sources(pdata, tm, dtm, du(:,:,:,:)) - - pdata%u1(:,:,:,:) = f1 * pdata%u2(:,:,:,:) & - + f2 * (pdata%u1(:,:,:,:) + ds * du(:,:,:,:)) - -#if NDIMS == 3 - lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g2 * du(:,nb:ne,nb:ne,nb:ne) -#else /* NDIMS == 3 */ - lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g2 * du(:,nb:ne,nb:ne, : ) -#endif /* NDIMS == 3 */ - l = l + 1 - - pdata => pdata%next - end do ! over data blocks - - call update_variables(tm, dtm) - != 6th step: U(1) = U(1) + dt/r F[U(1)], for i = n*(n+1)/2, ..., m ! - do i = n4, stages + do i = n4, stages - tm = time + (i - n) * ds - dtm = ds + tm = time + (i - n) * ds + dtm = ds - call update_fluxes() + call update_fluxes() - l = 1 - pdata => list_data - do while (associated(pdata)) + l = 1 + pdata => list_data + do while (associated(pdata)) - call update_increment(pdata) - call update_sources(pdata, tm, dtm, du(:,:,:,:)) + call update_increment(pdata) + call update_sources(pdata, tm, dtm, du(:,:,:,:)) - pdata%u1(:,:,:,:) = pdata%u1(:,:,:,:) + ds * du(:,:,:,:) + pdata%u1(:,:,:,:) = pdata%u1(:,:,:,:) + ds * du(:,:,:,:) #if NDIMS == 3 - lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g1 * du(:,nb:ne,nb:ne,nb:ne) + lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g1 * du(:,nb:ne,nb:ne,nb:ne) #else /* NDIMS == 3 */ - lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g1 * du(:,nb:ne,nb:ne, : ) + lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g1 * du(:,nb:ne,nb:ne, : ) #endif /* NDIMS == 3 */ - l = l + 1 + l = l + 1 - pdata => pdata%next - end do ! over data blocks + pdata => pdata%next + end do ! over data blocks - if (i < stages) call update_variables(tm, dtm) + if (i < stages) call update_variables(tm, dtm) - end do ! i = n4, stages + end do ! i = n4, stages ! find umax ! - umax = 0.0d+00 - pdata => list_data - do while (associated(pdata)) - umax = max(umax, maxval(abs(pdata%u0(:,:,:,:))), & - maxval(abs(pdata%u1(:,:,:,:)))) - pdata => pdata%next - end do ! over data blocks + umax = 0.0d+00 + pdata => list_data + do while (associated(pdata)) + umax = max(umax, maxval(abs(pdata%u0(:,:,:,:))), & + maxval(abs(pdata%u1(:,:,:,:)))) + pdata => pdata%next + end do ! over data blocks ! get the maximum error of each variable ! - do l = 1, nf - errors(l) = dt * maxval(abs(lerr(:,l,:,:,:))) - end do + do l = 1, nf + errors(l) = dt * maxval(abs(lerr(:,l,:,:,:))) + end do #ifdef MPI - call reduce_maximum(umax) - call reduce_maximum(errors) + call reduce_maximum(umax) + call reduce_maximum(errors) #endif /* MPI */ ! calculate tolerance and time step ! - utol = atol + rtol * umax - dte = dt * (maxval(errors) / utol)**(-2.5d-01) + utol = atol + rtol * umax + emax = maxval(errors) / utol + + if (emax <= 1.0d+00 .or. nrej >= mrej) then + test = .false. + + errs(3) = errs(2) + errs(2) = errs(1) + errs(1) = emax + + dte = dt * errs(1)**(-0.2d+00) * errs(2)**(0.0505d+00) & + * errs(3)**(-0.025d+00) + else + errs(1) = emax + + dte = dt * errs(1)**(-0.2d+00) * errs(2)**(0.0505d+00) & + * errs(3)**(-0.025d+00) + dt = dte + + nrej = nrej + 1 ! rejection count in the current step + nrejections = nrejections + 1 + end if + + niterations = niterations + 1 + + end do if (allocated(lerr)) deallocate(lerr) From d72c6019b86326ecea2bac93ca51d7ae7f2e2df5 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 1 Sep 2020 16:57:36 -0300 Subject: [PATCH 18/20] EVOLUTION: Implement step repetition in evolve_ssprk2_m(). Signed-off-by: Grzegorz Kowal --- sources/evolution.F90 | 157 +++++++++++++++++++++++++----------------- 1 file changed, 94 insertions(+), 63 deletions(-) diff --git a/sources/evolution.F90 b/sources/evolution.F90 index 6fec4de..aef728d 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -951,8 +951,9 @@ module evolution type(block_data), pointer :: pdata - integer :: n, l - real(kind=8) :: tm, dtm, ds, umax, utol + logical :: test + integer :: n, l, nrej + real(kind=8) :: tm, dtm, ds, umax, utol, emax logical , save :: first = .true. real(kind=8), save :: ft, fl, fr, gt, gl @@ -978,35 +979,71 @@ module evolution end if - ds = ft * dt - l = get_dblocks() #if NDIMS == 3 allocate(lerr(l,nf,nc,nc,nc)) #else /* NDIMS == 3 */ allocate(lerr(l,nf,nc,nc, 1)) #endif /* NDIMS == 3 */ - lerr(:,:,:,:,:) = 0.0d+00 + + test = .true. + nrej = 0 + + do while(test) + + lerr(:,:,:,:,:) = 0.0d+00 + + ds = ft * dt != 1st step: U(1) = U(n), U(2) = U(n) ! - pdata => list_data - do while (associated(pdata)) + pdata => list_data + do while (associated(pdata)) - pdata%u1(:,:,:,:) = pdata%u0(:,:,:,:) - pdata%u2(:,:,:,:) = pdata%u0(:,:,:,:) + pdata%u1(:,:,:,:) = pdata%u0(:,:,:,:) + pdata%u2(:,:,:,:) = pdata%u0(:,:,:,:) - pdata%u => pdata%u1 + pdata%u => pdata%u1 - pdata => pdata%next - end do ! over data blocks + pdata => pdata%next + end do ! over data blocks != 2nd step: U(1) = U(1) + dt/(m-1) F[U(1)], for i = 1, ..., m-1 ! - do n = 1, stages - 1 + do n = 1, stages - 1 - tm = time + n * ds - dtm = ds + tm = time + n * ds + dtm = ds + + call update_fluxes() + + l = 1 + pdata => list_data + do while (associated(pdata)) + + call update_increment(pdata) + call update_sources(pdata, tm, dtm, du(:,:,:,:)) + + pdata%u1(:,:,:,:) = pdata%u1(:,:,:,:) + ds * du(:,:,:,:) + +#if NDIMS == 3 + lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + gt * du(:,nb:ne,nb:ne,nb:ne) +#else /* NDIMS == 3 */ + lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + gt * du(:,nb:ne,nb:ne, : ) +#endif /* NDIMS == 3 */ + l = l + 1 + + pdata => pdata%next + end do ! over data blocks + + call update_variables(tm, dtm) + + end do ! n = 1, stages - 1 + +!= the final step: U(1) = (m-1)/m U(1) + 1/m U(2) + dt/m F[U(1)] +! + tm = time + dt + dtm = fl * dt call update_fluxes() @@ -1017,74 +1054,68 @@ module evolution call update_increment(pdata) call update_sources(pdata, tm, dtm, du(:,:,:,:)) - pdata%u1(:,:,:,:) = pdata%u1(:,:,:,:) + ds * du(:,:,:,:) + pdata%u1(:,:,:,:) = fr * pdata%u1(:,:,:,:) & + + fl * (pdata%u2(:,:,:,:) + dt * du(:,:,:,:)) #if NDIMS == 3 - lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + gt * du(:,nb:ne,nb:ne,nb:ne) + lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + gl * du(:,nb:ne,nb:ne,nb:ne) #else /* NDIMS == 3 */ - lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + gt * du(:,nb:ne,nb:ne, : ) + lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + gl * du(:,nb:ne,nb:ne, : ) #endif /* NDIMS == 3 */ l = l + 1 pdata => pdata%next end do ! over data blocks - call update_variables(tm, dtm) - - end do ! n = 1, stages - 1 - -!= the final step: U(1) = (m-1)/m U(1) + 1/m U(2) + dt/m F[U(1)] -! - tm = time + dt - dtm = fl * dt - - call update_fluxes() - - l = 1 - pdata => list_data - do while (associated(pdata)) - - call update_increment(pdata) - call update_sources(pdata, tm, dtm, du(:,:,:,:)) - - pdata%u1(:,:,:,:) = fr * pdata%u1(:,:,:,:) & - + fl * (pdata%u2(:,:,:,:) + dt * du(:,:,:,:)) - -#if NDIMS == 3 - lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + gl * du(:,nb:ne,nb:ne,nb:ne) -#else /* NDIMS == 3 */ - lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + gl * du(:,nb:ne,nb:ne, : ) -#endif /* NDIMS == 3 */ - l = l + 1 - - pdata => pdata%next - end do ! over data blocks - ! find umax ! - umax = 0.0d+00 - pdata => list_data - do while (associated(pdata)) - umax = max(umax, maxval(abs(pdata%u0(:,:,:,:))), & - maxval(abs(pdata%u1(:,:,:,:)))) - pdata => pdata%next - end do ! over data blocks + umax = 0.0d+00 + pdata => list_data + do while (associated(pdata)) + umax = max(umax, maxval(abs(pdata%u0(:,:,:,:))), & + maxval(abs(pdata%u1(:,:,:,:)))) + pdata => pdata%next + end do ! over data blocks ! get the maximum error of each variable ! - do l = 1, nf - errors(l) = dt * maxval(abs(lerr(:,l,:,:,:))) - end do + do l = 1, nf + errors(l) = dt * maxval(abs(lerr(:,l,:,:,:))) + end do #ifdef MPI - call reduce_maximum(umax) - call reduce_maximum(errors) + call reduce_maximum(umax) + call reduce_maximum(errors) #endif /* MPI */ ! calculate tolerance and time step ! - utol = atol + rtol * umax - dte = dt * (maxval(errors) / utol)**(-2.5d-01) + utol = atol + rtol * umax + emax = maxval(errors) / utol + + if (emax <= 1.0d+00 .or. nrej >= mrej) then + test = .false. + + errs(3) = errs(2) + errs(2) = errs(1) + errs(1) = emax + + dte = dt * errs(1)**(-0.2d+00) * errs(2)**(0.0505d+00) & + * errs(3)**(-0.025d+00) + else + errs(1) = emax + + dte = dt * errs(1)**(-0.2d+00) * errs(2)**(0.0505d+00) & + * errs(3)**(-0.025d+00) + dt = dte + + nrej = nrej + 1 ! rejection count in the current step + nrejections = nrejections + 1 + end if + + niterations = niterations + 1 + + end do if (allocated(lerr)) deallocate(lerr) From 55e6d810d2480107709786acb706f32de2562b9f Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 1 Sep 2020 17:16:18 -0300 Subject: [PATCH 19/20] BOUNDARIES, MESH: Implement parabolic interpolation for prolongation. Signed-off-by: Grzegorz Kowal --- sources/boundaries.F90 | 161 +++++++++++++++++++++++------------------ sources/mesh.F90 | 39 +++++----- 2 files changed, 111 insertions(+), 89 deletions(-) diff --git a/sources/boundaries.F90 b/sources/boundaries.F90 index 080048d..b558726 100644 --- a/sources/boundaries.F90 +++ b/sources/boundaries.F90 @@ -5805,42 +5805,47 @@ module boundaries ! do p = 1, nv + dq(1) = 1.25d-01 * (qn(p,i+1,j,k) - qn(p,i-1,j,k)) + dq(2) = 1.25d-01 * (qn(p,i,j+1,k) - qn(p,i,j-1,k)) + dq(3) = 1.25d-01 * (qn(p,i,j,k+1) - qn(p,i,j,k-1)) + + if (positive(p) .and. qn(p,i,j,k) <= sum(abs(dq(1:NDIMS)))) then + + if (qn(p,i,j,k) <= 0.0d+00) then + write(error_unit,"('[',a,']: ',a,3i4,a)") trim(loc) & + , "Positive variable is not positive at (", i, j, k, " )" + + dq(:) = 0.0d+00 + else + ! calculate limited derivatives in all directions ! - dql = qn(p,i ,j,k) - qn(p,im1,j,k) - dqr = qn(p,ip1,j,k) - qn(p,i ,j,k) - dq(1) = limiter_prol(0.5d+00, dql, dqr) + dql = qn(p,i ,j,k) - qn(p,im1,j,k) + dqr = qn(p,ip1,j,k) - qn(p,i ,j,k) + dq(1) = limiter_prol(2.5d-01, dql, dqr) - dql = qn(p,i,j ,k) - qn(p,i,jm1,k) - dqr = qn(p,i,jp1,k) - qn(p,i,j ,k) - dq(2) = limiter_prol(0.5d+00, dql, dqr) + dql = qn(p,i,j ,k) - qn(p,i,jm1,k) + dqr = qn(p,i,jp1,k) - qn(p,i,j ,k) + dq(2) = limiter_prol(2.5d-01, dql, dqr) - dql = qn(p,i,j,k ) - qn(p,i,j,km1) - dqr = qn(p,i,j,kp1) - qn(p,i,j,k ) - dq(3) = limiter_prol(0.5d+00, dql, dqr) + dql = qn(p,i,j,k ) - qn(p,i,j,km1) + dqr = qn(p,i,j,kp1) - qn(p,i,j,k ) + dq(3) = limiter_prol(2.5d-01, dql, dqr) - if (positive(p)) then - if (qn(p,i,j,k) > 0.0d+00) then do while (qn(p,i,j,k) <= sum(abs(dq(1:NDIMS)))) dq(:) = 0.5d+00 * dq(:) end do - else - write(error_unit,"('[',a,']: ',a,3i4,a)") trim(loc) & - , "Positive variable is not positive at (", i, j, k, " )" - dq(:) = 0.0d+00 end if end if - dq(:) = 0.5d+00 * dq(:) - -! calculate the derivative terms +! calculate derivative terms ! dq1 = dq(1) + dq(2) + dq(3) dq2 = dq(1) - dq(2) - dq(3) dq3 = dq(1) - dq(2) + dq(3) dq4 = dq(1) + dq(2) - dq(3) -! prolong the face region to the output array +! prolong the face region ! qb(p,is,js,ks) = qn(p,i,j,k) - dq1 qb(p,it,js,ks) = qn(p,i,j,k) + dq2 @@ -6114,43 +6119,50 @@ module boundaries ! do p = 1, nv -! calculate limited derivatives in all directions -! - dql = qn(p,i ,j,k) - qn(p,im1,j,k) - dqr = qn(p,ip1,j,k) - qn(p,i ,j,k) - dq(1) = limiter_prol(0.5d+00, dql, dqr) - - dql = qn(p,i,j ,k) - qn(p,i,jm1,k) - dqr = qn(p,i,jp1,k) - qn(p,i,j ,k) - dq(2) = limiter_prol(0.5d+00, dql, dqr) - + dq(1) = 1.25d-01 * (qn(p,i+1,j,k) - qn(p,i-1,j,k)) + dq(2) = 1.25d-01 * (qn(p,i,j+1,k) - qn(p,i,j-1,k)) #if NDIMS == 3 - dql = qn(p,i,j,k ) - qn(p,i,j,km1) - dqr = qn(p,i,j,kp1) - qn(p,i,j,k ) - dq(3) = limiter_prol(0.5d+00, dql, dqr) + dq(3) = 1.25d-01 * (qn(p,i,j,k+1) - qn(p,i,j,k-1)) +#endif /* NDIMS == 3 */ + + if (positive(p) .and. qn(p,i,j,k) <= sum(abs(dq(1:NDIMS)))) then + + if (qn(p,i,j,k) <= 0.0d+00) then + write(error_unit,"('[',a,']: ',a,3i4,a)") trim(loc) & + , "Positive variable is not positive at (", i, j, k, " )" + + dq(:) = 0.0d+00 + else + +! calculate limited derivatives in all directions +! + dql = qn(p,i ,j,k) - qn(p,im1,j,k) + dqr = qn(p,ip1,j,k) - qn(p,i ,j,k) + dq(1) = limiter_prol(2.5d-01, dql, dqr) + + dql = qn(p,i,j ,k) - qn(p,i,jm1,k) + dqr = qn(p,i,jp1,k) - qn(p,i,j ,k) + dq(2) = limiter_prol(2.5d-01, dql, dqr) + +#if NDIMS == 3 + dql = qn(p,i,j,k ) - qn(p,i,j,km1) + dqr = qn(p,i,j,kp1) - qn(p,i,j,k ) + dq(3) = limiter_prol(2.5d-01, dql, dqr) #endif /* NDIMS == 3 */ - if (positive(p)) then - if (qn(p,i,j,k) > 0.0d+00) then do while (qn(p,i,j,k) <= sum(abs(dq(1:NDIMS)))) dq(:) = 0.5d+00 * dq(:) end do - else - write(error_unit,"('[',a,']: ',a,3i4,a)") trim(loc) & - , "Positive variable is not positive at (", i, j, k, " )" - dq(:) = 0.0d+00 end if end if - dq(:) = 0.5d+00 * dq(:) - #if NDIMS == 2 -! calculate the derivative terms +! calculate derivative terms ! dq1 = dq(1) + dq(2) dq2 = dq(1) - dq(2) -! prolong the edge region to the output array +! prolong the edge region ! qb(p,is,js,k ) = qn(p,i,j,k) - dq1 qb(p,it,js,k ) = qn(p,i,j,k) + dq2 @@ -6158,14 +6170,14 @@ module boundaries qb(p,it,jt,k ) = qn(p,i,j,k) + dq1 #endif /* NDIMS == 2 */ #if NDIMS == 3 -! calculate the derivative terms +! calculate derivative terms ! dq1 = dq(1) + dq(2) + dq(3) dq2 = dq(1) - dq(2) - dq(3) dq3 = dq(1) - dq(2) + dq(3) dq4 = dq(1) + dq(2) - dq(3) -! prolong the edge region to the output array +! prolong the edge region ! qb(p,is,js,ks) = qn(p,i,j,k) - dq1 qb(p,it,js,ks) = qn(p,i,j,k) + dq2 @@ -6177,7 +6189,7 @@ module boundaries qb(p,it,jt,kt) = qn(p,i,j,k) + dq1 #endif /* NDIMS == 3 */ - end do ! q = 1, nv + end do ! p = 1, nv end do ! i = il, iu end do ! j = jl, ju @@ -6391,43 +6403,50 @@ module boundaries ! do p = 1, nv -! calculate limited derivatives in all directions -! - dql = qn(p,i ,j,k) - qn(p,im1,j,k) - dqr = qn(p,ip1,j,k) - qn(p,i ,j,k) - dq(1) = limiter_prol(0.5d+00, dql, dqr) - - dql = qn(p,i,j ,k) - qn(p,i,jm1,k) - dqr = qn(p,i,jp1,k) - qn(p,i,j ,k) - dq(2) = limiter_prol(0.5d+00, dql, dqr) - + dq(1) = 1.25d-01 * (qn(p,i+1,j,k) - qn(p,i-1,j,k)) + dq(2) = 1.25d-01 * (qn(p,i,j+1,k) - qn(p,i,j-1,k)) #if NDIMS == 3 - dql = qn(p,i,j,k ) - qn(p,i,j,km1) - dqr = qn(p,i,j,kp1) - qn(p,i,j,k ) - dq(3) = limiter_prol(0.5d+00, dql, dqr) + dq(3) = 1.25d-01 * (qn(p,i,j,k+1) - qn(p,i,j,k-1)) +#endif /* NDIMS == 3 */ + + if (positive(p) .and. qn(p,i,j,k) <= sum(abs(dq(1:NDIMS)))) then + + if (qn(p,i,j,k) <= 0.0d+00) then + write(error_unit,"('[',a,']: ',a,3i4,a)") trim(loc) & + , "Positive variable is not positive at (", i, j, k, " )" + + dq(:) = 0.0d+00 + else + +! calculate limited derivatives in all directions +! + dql = qn(p,i ,j,k) - qn(p,im1,j,k) + dqr = qn(p,ip1,j,k) - qn(p,i ,j,k) + dq(1) = limiter_prol(2.5d-01, dql, dqr) + + dql = qn(p,i,j ,k) - qn(p,i,jm1,k) + dqr = qn(p,i,jp1,k) - qn(p,i,j ,k) + dq(2) = limiter_prol(2.5d-01, dql, dqr) + +#if NDIMS == 3 + dql = qn(p,i,j,k ) - qn(p,i,j,km1) + dqr = qn(p,i,j,kp1) - qn(p,i,j,k ) + dq(3) = limiter_prol(2.5d-01, dql, dqr) #endif /* NDIMS == 3 */ - if (positive(p)) then - if (qn(p,i,j,k) > 0.0d+00) then do while (qn(p,i,j,k) <= sum(abs(dq(1:NDIMS)))) dq(:) = 0.5d+00 * dq(:) end do - else - write(error_unit,"('[',a,']: ',a,3i4,a)") trim(loc) & - , "Positive variable is not positive at (", i, j, k, " )" - dq(:) = 0.0d+00 end if end if - dq(:) = 0.5d+00 * dq(:) - #if NDIMS == 2 -! calculate the derivative terms +! calculate derivative terms ! dq1 = dq(1) + dq(2) dq2 = dq(1) - dq(2) -! prolong the corner region to the output array +! prolong the corner region ! qb(p,is,js,k ) = qn(p,i,j,k) - dq1 qb(p,it,js,k ) = qn(p,i,j,k) + dq2 @@ -6435,14 +6454,14 @@ module boundaries qb(p,it,jt,k ) = qn(p,i,j,k) + dq1 #endif /* NDIMS == 2 */ #if NDIMS == 3 -! calculate the derivative terms +! calculate derivative terms ! dq1 = dq(1) + dq(2) + dq(3) dq2 = dq(1) - dq(2) - dq(3) dq3 = dq(1) - dq(2) + dq(3) dq4 = dq(1) + dq(2) - dq(3) -! prolong the corner region to the output array +! prolong the corner region ! qb(p,is,js,ks) = qn(p,i,j,k) - dq1 qb(p,it,js,ks) = qn(p,i,j,k) + dq2 @@ -6454,7 +6473,7 @@ module boundaries qb(p,it,jt,kt) = qn(p,i,j,k) + dq1 #endif /* NDIMS == 3 */ - end do ! q = 1, nv + end do ! p = 1, nv end do ! i = il, iu end do ! j = jl, ju diff --git a/sources/mesh.F90 b/sources/mesh.F90 index 12fd75f..b3e8f8c 100644 --- a/sources/mesh.F90 +++ b/sources/mesh.F90 @@ -1190,36 +1190,39 @@ module mesh do p = 1, nv - dul = pdata%u(p,i ,j,k) - pdata%u(p,i-1,j,k) - dur = pdata%u(p,i+1,j,k) - pdata%u(p,i ,j,k) - du(1) = limiter_prol(0.5d+00, dul, dur) - - dul = pdata%u(p,i,j ,k) - pdata%u(p,i,j-1,k) - dur = pdata%u(p,i,j+1,k) - pdata%u(p,i,j ,k) - du(2) = limiter_prol(0.5d+00, dul, dur) - + du(1) = 1.25d-01 * (pdata%u(p,i+1,j,k) - pdata%u(p,i-1,j,k)) + du(2) = 1.25d-01 * (pdata%u(p,i,j+1,k) - pdata%u(p,i,j-1,k)) #if NDIMS == 3 - dul = pdata%u(p,i,j,k ) - pdata%u(p,i,j,k-1) - dur = pdata%u(p,i,j,k+1) - pdata%u(p,i,j,k ) - du(3) = limiter_prol(0.5d+00, dul, dur) + du(3) = 1.25d-01 * (pdata%u(p,i,j,k+1) - pdata%u(p,i,j,k-1)) #endif /* NDIMS == 3 */ - if (positive(p)) then + if (positive(p) .and. pdata%u(p,i,j,k) < sum(abs(du(1:NDIMS)))) then if (pdata%u(p,i,j,k) > 0.0d+00) then + + dul = pdata%u(p,i ,j,k) - pdata%u(p,i-1,j,k) + dur = pdata%u(p,i+1,j,k) - pdata%u(p,i ,j,k) + du(1) = limiter_prol(2.5d-01, dul, dur) + + dul = pdata%u(p,i,j ,k) - pdata%u(p,i,j-1,k) + dur = pdata%u(p,i,j+1,k) - pdata%u(p,i,j ,k) + du(2) = limiter_prol(2.5d-01, dul, dur) + +#if NDIMS == 3 + dul = pdata%u(p,i,j,k ) - pdata%u(p,i,j,k-1) + dur = pdata%u(p,i,j,k+1) - pdata%u(p,i,j,k ) + du(3) = limiter_prol(2.5d-01, dul, dur) +#endif /* NDIMS == 3 */ + do while (pdata%u(p,i,j,k) <= sum(abs(du(1:NDIMS)))) du(:) = 0.5d+00 * du(:) end do else - write(error_unit,"('[',a,']: ',a,i9,a,3i4,a)") trim(loc) & - , "Positive variable is not positive for block ID:" & - , pmeta%id, " at ( ", i, j, k, " )" + write(error_unit,"('[',a,']: ',a,3i4,a)") trim(loc) & + , "Positive variable is not positive at (", i, j, k, " )" du(:) = 0.0d+00 - status = 1 end if end if - du(:) = 0.5d+00 * du(:) - #if NDIMS == 2 du1 = du(1) + du(2) du2 = du(1) - du(2) From d82eb99a8a5213c1e5ca9d8ad4bb978dd9867ac2 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 1 Sep 2020 18:12:37 -0300 Subject: [PATCH 20/20] EVOLUTION: Fix variables which may not be initialized. Signed-off-by: Grzegorz Kowal --- sources/evolution.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/sources/evolution.F90 b/sources/evolution.F90 index aef728d..c5b33aa 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -1122,7 +1122,7 @@ module evolution != final step: U(n+1) = U(1) ! tm = time + dt - dtm = ds + dtm = ft * dt pdata => list_data do while (associated(pdata)) @@ -2191,7 +2191,7 @@ module evolution != final step: U(n+1) = U(1) ! tm = time + dt - dtm = ds + dtm = dt / r pdata => list_data do while (associated(pdata)) @@ -2207,7 +2207,7 @@ module evolution ! update ψ with its source term in the last step ! - if (i == stages .and. ibp > 0) then + if (ibp > 0) then pdata => list_data do while (associated(pdata)) pdata%u(ibp,:,:,:) = decay * pdata%u(ibp,:,:,:)