From 27fd2fb6e141b6fea2fe3af70e613e810d7adf96 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Thu, 16 Dec 2021 09:18:32 -0300 Subject: [PATCH] EVOLUTION: Improve the initial time step estimation. First of all, the estimation is of the third order, so use the cubic square in the calculation of h1. Secondly, fix the normalization factor fnorm. Finally, use only the variables corresponding to the fluxes in the increment norms sc(..) and df(..). Signed-off-by: Grzegorz Kowal --- sources/evolution.F90 | 53 ++++++++++++++++++++++++------------------- 1 file changed, 30 insertions(+), 23 deletions(-) diff --git a/sources/evolution.F90 b/sources/evolution.F90 index 28220d4..613b0a3 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -1240,7 +1240,7 @@ module evolution allocate(s(n,3)) s(:,:) = 0.0d+00 - fnorm = 1.0d+00 / (get_nleafs() * nc**NDIMS) + fnorm = 1.0d+00 / (get_nleafs() * size(sc)) d(:) = 0.0d+00 @@ -1252,17 +1252,18 @@ module evolution 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)) + sc(:,:,:,:) = atol + rtol * abs(pdata%uu(1:nf,nb:ne,nb:ne,nb:ne,1)) - s(l,1) = sum((pdata%uu(:,nb:ne,nb:ne,nb:ne,1) / sc(:,:,:,:))**2) - s(l,2) = sum((pdata%du(:,nb:ne,nb:ne,nb:ne) / sc(:,:,:,:))**2) + s(l,1) = sum((pdata%uu(1:nf,nb:ne,nb:ne,nb:ne,1) / sc(:,:,:,:))**2) + s(l,2) = sum((pdata%du(1:nf,nb:ne,nb:ne,nb:ne) / sc(:,:,:,:))**2) #else /* NDIMS == 3 */ - sc(:,:,:,:) = atol + rtol * abs(pdata%uu(:,nb:ne,nb:ne, : ,1)) + sc(:,:,:,:) = atol + rtol * abs(pdata%uu(1:nf,nb:ne,nb:ne, : ,1)) - s(l,1) = sum((pdata%uu(:,nb:ne,nb:ne, : ,1) / sc(:,:,:,:))**2) - s(l,2) = sum((pdata%du(:,nb:ne,nb:ne, : ) / sc(:,:,:,:))**2) + s(l,1) = sum((pdata%uu(1:nf,nb:ne,nb:ne, : ,1) / sc(:,:,:,:))**2) + s(l,2) = sum((pdata%du(1:nf,nb:ne,nb:ne, : ) / sc(:,:,:,:))**2) #endif /* NDIMS == 3 */ + pdata%u => pdata%uu(:,:,:,:,2) end do !$omp end parallel do d(1) = sum(s(:,1)) @@ -1290,8 +1291,6 @@ module evolution pdata => data_blocks(l)%ptr pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) + h0 * pdata%du(:,:,:,:) - - pdata%u => pdata%uu(:,:,:,:,2) end do !$omp end parallel do @@ -1317,41 +1316,36 @@ module evolution pdata => data_blocks(l)%ptr #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(1:nf,nb:ne,nb:ne,nb:ne) + sc(:,:,:,:) = atol + rtol * abs(pdata%uu(1:nf,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(1:nf,nb:ne,nb:ne, : ) + sc(:,:,:,:) = atol + rtol * abs(pdata%uu(1:nf,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) + df(:,:,:,:) = df(:,:,:,:) - pdata%du(1:nf,nb:ne,nb:ne,nb:ne) #else /* NDIMS == 3 */ - df(:,:,:,:) = df(:,:,:,:) - pdata%du(:,nb:ne,nb:ne, : ) + df(:,:,:,:) = df(:,:,:,:) - pdata%du(1:nf,nb:ne,nb:ne, : ) #endif /* NDIMS == 3 */ s(l,3) = sum((df(:,:,:,:) / sc(:,:,:,:))**2) - - pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) - - pdata%u => pdata%uu(:,:,:,:,1) end do !$omp end parallel do d(3) = sum(s(:,3)) deallocate(s) - call update_variables(time, 0.0d+00, status) - #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)) + h1 = max(d(2), d(3)) + if (h1 > 1.0d-15) then + h1 = (1.0d-02 / h1)**(1.0d+00 / 3.0d+00) else h1 = max(1.0d-06, 1.0d-03 * h0) end if @@ -1361,6 +1355,19 @@ module evolution end if +! revert the initial state +! +!$omp parallel do default(shared) private(pdata,df,sc) + do l = 1, n + pdata => data_blocks(l)%ptr + pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) + + pdata%u => pdata%uu(:,:,:,:,1) + end do +!$omp end parallel do + + call update_variables(time, 0.0d+00, status) + end if !-------------------------------------------------------------------------------