EVOLUTION: Update subroutine to calculate errors using L2-norm.

Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>

sources/evolution.F90

@@ -4304,13 +4304,16 @@ module evolution
 ! subroutine UPDATE_ERRORS_L2:
 ! ---------------------------
 !
-!   Subroutine updates the errors with respect to the absolute and relative
-!   tolerances using the L2 norm. The errors are calculated per variable.
+! Description:
+!
+!   The subroutine is designed to update errors with respect to absolute and
+!   relative tolerances using the L2 norm. The errors are calculated for
+!   each variable.
 !
 ! Arguments:
 !
-!     nh - the index of the higher order solution
-!     nl - the index of the lower order solution
+!   nh - the index of the higher order solution;
+!   nl - the index of the lower order solution;
 !
 !===============================================================================
 !
@@ -4318,62 +4321,101 @@ module evolution
 
     use blocks     , only : block_data, data_blocks
     use blocks     , only : get_nleafs, get_dblocks
-    use coordinates, only : ncells, nb, ne
+    use coordinates, only : nc => ncells, nb, ne
     use equations  , only : nf, errors
+    use helpers    , only : print_message
 #ifdef MPI
     use mpitools   , only : reduce_sum
 #endif /* MPI */
+    use workspace  , only : resize_workspace, work, work_in_use
 
     implicit none
 
     integer, intent(in) :: nh, nl
 
+    character(len=*), parameter :: loc = 'EVOLUTION:update_errors_l2()'
+
     type(block_data), pointer :: pdata
 
-    integer      :: l, n, p
-    real(kind=8) :: err
+    integer            :: nblks, n, p
 
     logical     , save :: first = .true.
+    integer     , save :: nt    = 0
     real(kind=8), save :: fnorm = 1.0d+00
 
+    real(kind=8), dimension(nf)                   :: lerrors
+    real(kind=8), dimension(:,:,:), pointer, save :: uhi, ulo, div
+
+!$omp threadprivate(first,nt,uhi,ulo,div)
+!$  integer :: omp_get_thread_num
+
 !-------------------------------------------------------------------------------
 !
+    errors(:) = 0.0d+00
+
+    nblks = get_dblocks()
+
+!$omp parallel default(shared) private(pdata,lerrors,n,p)
     if (first) then
-      fnorm = 1.0d+00 / ncells**NDIMS
+!$    nt = omp_get_thread_num()
+      n  = nc**NDIMS
+
+#if NDIMS == 3
+      uhi(1:nc,1:nc,1:nc) => work(    1:  n,nt)
+      ulo(1:nc,1:nc,1:nc) => work(  n+1:2*n,nt)
+      div(1:nc,1:nc,1:nc) => work(2*n+1:3*n,nt)
+#else /* NDIMS == 3 */
+      uhi(1:nc,1:nc,1: 1) => work(    1:  n,nt)
+      ulo(1:nc,1:nc,1: 1) => work(  n+1:2*n,nt)
+      div(1:nc,1:nc,1: 1) => work(2*n+1:3*n,nt)
+#endif /* NDIMS == 3 */
+
+!$omp single
+      fnorm = 1.0d+00 / n
+!$omp end single
 
       first = .false.
     end if
 
-    errors(:) = 0.0d+00
+    if (work_in_use(nt)) &
+      call print_message(loc, "Workspace is being used right now! " // &
+                              "Corruptions can occur!")
 
-    n = get_dblocks()
+    work_in_use(nt) = .true.
+!$omp barrier
+
+!$omp do reduction(+:errors)
+    do n = 1, nblks
+      pdata => data_blocks(n)%ptr
 
-!$omp parallel do default(shared) private(pdata,p,err)
-    do l = 1, n
-      pdata => data_blocks(l)%ptr
       do p = 1, nf
 #if NDIMS == 3
-        err = sum((abs(pdata%uu(p,nb:ne,nb:ne,nb:ne,nh) &
-                     - pdata%uu(p,nb:ne,nb:ne,nb:ne,nl)) / (atol + rtol * &
-               max(abs(pdata%uu(p,nb:ne,nb:ne,nb:ne,nh)), &
-                   abs(pdata%uu(p,nb:ne,nb:ne,nb:ne,nl)))))**2)
+        uhi = pdata%uu(p,nb:ne,nb:ne,nb:ne,nh)
+        ulo = pdata%uu(p,nb:ne,nb:ne,nb:ne,nl)
 #else /* NDIMS == 3 */
-        err = sum((abs(pdata%uu(p,nb:ne,nb:ne,  :  ,nh) &
-                     - pdata%uu(p,nb:ne,nb:ne,  :  ,nl)) / (atol + rtol * &
-               max(abs(pdata%uu(p,nb:ne,nb:ne,  :  ,nh)), &
-                   abs(pdata%uu(p,nb:ne,nb:ne,  :  ,nl)))))**2)
+        uhi = pdata%uu(p,nb:ne,nb:ne,  :  ,nh)
+        ulo = pdata%uu(p,nb:ne,nb:ne,  :  ,nl)
 #endif /* NDIMS == 3 */
-!$omp atomic update
-        errors(p) = errors(p) + fnorm * err
+
+        div = atol + rtol * max(abs(uhi), abs(ulo))
+
+        lerrors(p) = fnorm * sum((abs(uhi - ulo) / div)**2)
       end do
+
+      errors = errors + lerrors
     end do
-!$omp end parallel do
+!$omp end do nowait
+
+    work_in_use(nt) = .false.
+!$omp end parallel
+
+    errors(:) = errors(:) / get_nleafs()
 
 #ifdef MPI
     call reduce_sum(errors)
 #endif /* MPI */
 
-    errors(:) = sqrt(errors(:) / get_nleafs())
+    errors(:) = sqrt(errors(:))
 
 !-------------------------------------------------------------------------------
 !