documented, benchmarked, increased accuracy casinggrid

src/bnorml.f90

@@ -85,7 +85,7 @@ subroutine bnorml( mn, Ntz, efmn, ofmn )
 
   use fileunits, only : ounit, lunit
 
-  use inputlist, only : Wmacros, Wbnorml, Igeometry, Lcheck, vcasingtol, vcasingper, Lrad, Lvcgrid  
+  use inputlist, only : Wmacros, Wbnorml, Igeometry, Lcheck, vcasingtol, vcasingits, vcasingper, Lrad, Lvcgrid  
 
   use cputiming, only : Tbnorml
 
@@ -111,7 +111,7 @@ subroutine bnorml( mn, Ntz, efmn, ofmn )
   INTEGER              :: lvol, Lparallel, ii, jj, kk, jk, ll, kkmodnp, jkmodnp, ifail, id01daf, nvccalls, icasing
   REAL                 :: zeta, teta, gBn
   REAL                 :: Bxyz(1:Ntz,1:3), distance(1:Ntz)
-  REAL                 :: vcgriderr, resulth, resulth2, resulth4, deltah4h2, deltah2h 
+  REAL                 :: absvcerr, relvcerr, resulth, resulth2, resulth4, deltah4h2, deltah2h
   INTEGER              :: vcstride, Nzwithsym 
 
 
@@ -150,27 +150,27 @@ subroutine bnorml( mn, Ntz, efmn, ofmn )
 #else
 if ( Lvcgrid.eq.1 ) then
 #endif
-  ! Precompute Jxyz(1:Ntz,1:3) and the corresponding positions on the high resolution plasma boundary
+  ! Precompute Jxyz(:,1:3) and the corresponding positions on the high resolution plasma boundary
 
   !$OMP PARALLEL DO SHARED(Pbxyz, Jxyz) PRIVATE(jk, jj, kk, teta, zeta) COLLAPSE(2)
   do kk = 0, vcNz-1 ; 
     do jj = 0, vcNt-1 ; 
-      zeta = kk * pi2nfp / vcNz
+      zeta = kk * pi2 / vcNz ! zeta = kk * pi2nfp / vcNz
       teta = jj * pi2  / vcNt ; 
       jk = 1 + jj + kk*vcNt
 
       ! Each MPI rank only computes every a 1/ncpu surfacecurrent() calls 
       select case( Lparallel ) 
       case( 0 ) ! Lparallel = 0 
-       if( myid.ne.modulo(kk,ncpu) ) cycle
+      if( myid.ne.modulo(kk,ncpu) ) cycle
       case( 1 ) ! Lparallel = 1 
-       if( myid.ne.modulo(jk-1,ncpu) ) cycle
+      if( myid.ne.modulo(jk-1,ncpu) ) cycle
       case default ! Lparallel
         FATAL( bnorml, .true., invalid Lparallel in parallelization loop )
       end select 
 
       ! pbxyz and jxyz are  both [out] parameters
-      call surfacecurrent( teta, zeta, Pbxyz(jk,1:3), Jxyz(jk,1:3) )
+      call surfacecurrent( teta, zeta, Pbxyz(jk,1:3), Jxyz(jk,1:3)  )
     enddo
   enddo
 
@@ -186,8 +186,13 @@ subroutine bnorml( mn, Ntz, efmn, ofmn )
   endif
 
   ! iterate over resolutions of the virtual casing grid to get an estimate of the accuracy. Write the result into ijimag
-  do vcstride = 3, 0, -1
-    !$OMP PARALLEL DO SHARED(Dxyz, Nxyz, Pbxyz, Jxyz, ijreal, ijimag) PRIVATE(jk, gBn) COLLAPSE(2)
+  deltah2h = 1e12
+  !> The surface currents get precomputed at all points, but maybe we don't have to integrate over the whole grid. 
+  !> Start with a large stride over the plasmaboundary (= how many points to skip at each step), then get progressively finer.
+  !> Do at least vcasingits iterations (sqrt(vcasingits) resolution per dimension) and halve the stride until the accuracy is reached. 
+  !> At least one step in the resolution cascade is required to determine an estimate of the current accuracy.
+  do vcstride = INT(log(real(vcNt/sqrt(real(vcasingits))))/log(2.0)), 0, -1
+    !$OMP PARALLEL DO SHARED(Dxyz, Nxyz, Pbxyz, Jxyz, ijreal, ijimag) FIRSTPRIVATE(jk, gBn) COLLAPSE(2)
     do kk = 0, Nzwithsym-1 ; 
       do jj = 0, Nt-1 ; 
         jk = 1 + jj + kk*Nt
@@ -206,28 +211,30 @@ subroutine bnorml( mn, Ntz, efmn, ofmn )
         FATAL( bnorml, .true., invalid Lparallel in parallelization loop )
         end select ! end of select case( Lparallel ) 
 
-        call casinggrid( Dxyz(:,jk), Nxyz(:,jk), Pbxyz, Jxyz, 2**vcstride,  gBn)
+        gBn = zero
+        globaljk = jk ! only to check against the precomputed values against dvcfield
+        call casinggrid( Dxyz(1:3,jk), Nxyz(1:3,jk), Pbxyz, Jxyz, 2**vcstride,  gBn)
 
         ijimag(jk) = ijreal(jk) ! previous solution (lower resolution)
         ijreal(jk) = gBn ! current solution (higher resolution)
-      enddo
-    enddo
+      enddo ! end of do jj
+    enddo ! end of do kk
 
     deltah4h2 = deltah2h
     deltah2h =  sum(abs(ijimag - ijreal)) ! mean delta between the h and h/2 solutions
 
     ! Order of the integration method: log(deltah4h2/deltah2h)/log(2.0) = 1
-    ! relative error: deltah2h/abs(ijimag(jk))
-    ! absolute error: delta2h/Ntz
-    vcgriderr = deltah2h / sum(abs(ijreal))
-  enddo
-
-  if (myid.eq.0) then
-    write(ounit, '("bnorml : ", 10x ," : vcgriderr = ",es13.5," ; vcasingtol = ",es13.5s)') vcgriderr, vcasingtol
-  endif
-  ! if (vcgriderr.gt.vcasingtol) then
-  !   FATAL( bnorml, .true., virtual casing accuracy is too low, increase vcNt and vcNz )
-  ! endif
+    absvcerr = deltah2h / sum(abs(ijreal))
+    relvcerr = deltah2h / Ntz
+    if (myid.eq.0) then
+      write(ounit, '("bnorml : ", 10x ," : relvcerr = ",es13.5," ; absvcerr = ",es13.5," ; vcasingtol = ",es13.5s)') relvcerr, absvcerr, vcasingtol
+      ! print *, "Convergence order: ",  log(deltah4h2/deltah2h)/log(2.0)
+    endif
+    ! Tolerance was already reached, exit the loop. (Different threads may exit at different times)
+    if (absvcerr.lt.vcasingtol .and. relvcerr.lt.vcasingtol) then
+      exit
+    endif
+  enddo ! end of do vcstride
 #ifdef COMPARECASING
   ! To compare with the casing() implementation, copy the results into ijimag
   if(Lvcgrid.eq.1) then
@@ -332,7 +339,6 @@ subroutine bnorml( mn, Ntz, efmn, ofmn )
       enddo
     enddo
   endif
-
 !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
 
   ijreal(1:Ntz) = ijreal(1:Ntz) * virtualcasingfactor

src/casing.f90

@@ -224,36 +224,41 @@ end subroutine casing
 !> with \f${\bf r}\f$ approximately the distance vector from an evaluation point of the current surface \f${\bf x}(\theta, \zeta)\f$ to the 
 !> external point \f$(x,y,z)\f$, and the exact formula includes a regularization factor \f$\epsilon\f$  Eqn.\f$(\ref{eq:vcasing_distance})\f$.
 !>
+!> Numerically it uses a Kahan Summation algorithm to accumulate the contributions of all inner points, because truncation error quickly dominates the computation when thousands 
+!> of points are used in both dimensions.
+!> 
 !> @param[in] xyz \f$(x,y,z)\f$ cartesian coordinates on the computational boundary
-!> @param[in] nxyz cartesian normal vector on the computational boundary
+!> @param[in] ntz cartesian normal vector on the computational boundary
 !> @param[in] Pbxyz  array of cartesian coordinates on the plasma boundary
 !> @param[in] Jxyz array of surface currents \f${\bf B}_{Plasma} \cdot {\bf e}_\theta \times {\bf e}_\zeta \;\f$ on the plasma boundary
 !> @param[in] vcstride integer stride for the fixed resolution grid. vcstride = 1 computes the integral at full resolution, vcstride = 2 computes the integral at half resolution, etc.
 !> @param[out] gBn normal field \f${\bf B}_{Plasma} \cdot n \;\f$ on the computational boundary
-subroutine casinggrid( xyz, nxyz,  Pbxyz, Jxyz, vcstride, gBn)
+subroutine casinggrid( xyz, ntz,  Pbxyz, Jxyz, vcstride, gBn)
   use constants, only : zero, one, three, pi2
 
   use fileunits, only : ounit, vunit
 
   use inputlist, only : vcasingeps, vcNz, vcNt
 
-  use allglobal, only : myid, ncpu, cpus, MPI_COMM_SPEC !, Pbxyz, Jxyz
+  use allglobal, only : myid, ncpu, cpus, MPI_COMM_SPEC,  Dxyz, Nxyz,globaljk,pi2nfp
 
   LOCALS
 
-  REAL, intent(in)     :: xyz(3) ! arbitrary location; Cartesian;
-  REAL, intent(in)     :: nxyz(3) ! surface normal on the computational boundary; Cartesian;
+  REAL, intent(in)     :: xyz(1:3) ! arbitrary location; Cartesian;
+  REAL, intent(in)     :: ntz(1:3) ! surface normal on the computational boundary; Cartesian;
   REAL, intent(in)     :: Pbxyz(1:vcNz*vcNt, 1:3), Jxyz(1:vcNz*vcNt, 1:3) 
   INTEGER, intent(in)  :: vcstride
   REAL, intent(out)    :: gBn ! B.n on the computational boundary;
 
   REAL :: rr(1:3),  distance(1:3), jj(1:3), Bxyz(1:3), accumulator, firstorderfactor
+  REAL :: gBnlocal ! Local accumulator so we don't directly sum into the destination pointer (Avoid cache thrashing due to false share) 
+  REAL :: c, t, y ! Helper variables for the Kahan summation
   INTEGER :: plasmaNtz, jk
-
 !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
   ! plasmaNtz = SIZE(Pbxyz, 1)
   plasmaNtz = vcNz*vcNt
-  gBn = zero
+  gBnlocal = zero
+  c = zero
   ! loop over the high resolution plasma boundary (inner boundary for virtual casing)
   do jk = 1, plasmaNtz, vcstride ;
       ! distance vector between point on computational boundary and point on plasma boundary
@@ -275,9 +280,13 @@ subroutine casinggrid( xyz, nxyz,  Pbxyz, Jxyz, vcstride, gBn)
                      jj(1) * rr(2) - jj(2) * rr(1)  /)
 
       ! Accumulate B.r/r^3 contributions 
-      gBn = gBn + sum( Bxyz * nxyz ) * firstorderfactor
+      ! Kahan Summation algorithm
+      y = sum( Bxyz * ntz ) * firstorderfactor - c  ! Correct the next value with the compensation
+      t = gBnlocal + y                                  ! Add the corrected value to the sum
+      c = (t - gBnlocal) - y                            ! Update compensation
+      gBnlocal = t                                      ! Update the running total      
  enddo
- gBn = gBn * pi2 * pi2 / (plasmaNtz / vcstride) 
+ gBn = gBnlocal * pi2 * pi2 / (plasmaNtz / vcstride) 
  return
 
 end subroutine casinggrid