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 !-------------------------------------------------------------------------------