Grid virtual casing

src/bnorml.f90

@@ -85,18 +85,21 @@ subroutine bnorml( mn, Ntz, efmn, ofmn )
 
   use fileunits, only : ounit, lunit
 
-  use inputlist, only : Wmacros, Wbnorml, Igeometry, Lcheck, vcasingtol, vcasingper, Lrad
+  use inputlist, only : Wmacros, Wbnorml, Igeometry, Lcheck, vcasingtol, vcasingper, Lrad, Lvcgrid  
 
   use cputiming, only : Tbnorml
 
+  use inputlist, only : vcNz, vcNt
+
   use allglobal, only : ncpu, myid, cpus, MPI_COMM_SPEC, pi2nfp, Mvol, &
                         Nt, Nz, &
                         Rij, Zij, guvij, sg, TT, &
                         NOTstellsym, Lcoordinatesingularity, &
                         im, in, Ate, Aze, Ato, Azo, &
                         Nt, Nz, cfmn, sfmn, &
                         ijreal, ijimag, jireal, jiimag, &
-                        globaljk, tetazeta, virtualcasingfactor, gteta, gzeta, Dxyz, Nxyz
+                        globaljk, virtualcasingfactor, gteta, gzeta, Dxyz, Nxyz, &
+                        Jxyz, Pbxyz
 
 !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
 
@@ -105,12 +108,12 @@ subroutine bnorml( mn, Ntz, efmn, ofmn )
   INTEGER, intent(in)  :: mn, Ntz
   REAL   , intent(out) :: efmn(1:mn), ofmn(1:mn)
 
-  INTEGER              :: lvol, Lcurvature, Lparallel, ii, jj, kk, jk, ll, kkmodnp, jkmodnp, ifail, id01daf, nvccalls, icasing, ideriv
-  REAL                 :: lss, zeta, teta, cszeta(0:1), tetalow, tetaupp, absacc, gBn
-  REAL                 :: Jxyz(1:Ntz,1:3), Bxyz(1:Ntz,1:3), dAt(1:Ntz), dAz(1:Ntz), distance(1:Ntz)
+  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                 :: accuracyestimate, resulth, resulth2, resulth4, deltah4h2, deltah2h 
+  INTEGER              :: vcstep
 
- !REAL                 :: vcintegrand, zetalow, zetaupp
-! external             :: vcintegrand, zetalow, zetaupp
 
   BEGIN(bnorml)
 
@@ -121,10 +124,6 @@ subroutine bnorml( mn, Ntz, efmn, ofmn )
 !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
 
   ijreal(1:Ntz) = zero ! normal plasma field; 15 Oct 12;
- !ijimag(1:Ntz) = zero
-
- !jireal(1:Ntz) = zero
- !jiimag(1:Ntz) = zero
 
 !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
 
@@ -138,12 +137,62 @@ subroutine bnorml( mn, Ntz, efmn, ofmn )
 
 !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
 
-  do kk = 0, Nz-1 ; zeta = kk * pi2nfp / Nz
+#ifdef COMPARECASING
+! When comparing results, both methods should always run
+if (.true.) then
+#else
+if ( Lvcgrid.eq.1 ) then
+#endif
+  ! Precompute Jxyz(1:Ntz,1:3) and the corresponding positions on the high resolution plasma boundary
+
+  !$OMP PARALLEL DO SHARED(Pbxyz, Jxyz) PRIVATE(jk, teta, zeta) COLLAPSE(2)
+  do kk = 0, vcNz-1 ; 
+    do jj = 0, vcNt-1 ; 
+      zeta = kk * pi2nfp / vcNz
+      teta = jj * pi2  / vcNt ; 
+      jk = 1 + jj + kk*vcNt
+
+      call dvcfieldimpl( teta, zeta, Pbxyz(jk,1:3), Jxyz(jk,1:3) )
+
+    enddo
+  enddo
 
-   if( Igeometry.eq.3 ) then ; cszeta(0:1) = (/ cos(zeta), sin(zeta) /)
-   endif
+  ! iterate over resolutions of the virtual casing grid to get an estimate of the accuracy
+  do vcstep = 3, 0, -1
+    !$OMP PARALLEL DO SHARED(Dxyz, Pbxyz, Jxyz, ijimag) PRIVATE(jk, gBn)
+    do jk = 1, Ntz
+      call casing2( Dxyz(:,jk), Nxyz(:,jk), Pbxyz, Jxyz, 2**vcstep,  gBn)
+
+      ijreal(jk) = ijimag(jk) ! previous solution (lower resolution)
+      ijimag(jk) = gBn ! current solution (higher resolution)
 
-   do jj = 0, Nt-1 ; teta = jj * pi2    / Nt ; jk = 1 + jj + kk*Nt
+    enddo
+    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
+    accuracyestimate = deltah2h / sum(abs(ijimag))
+  enddo
+
+  write(ounit, '("bnorml : ", 10x ," : accuracyestimate = ",es13.5," ; vcasingtol = ",es13.5)') accuracyestimate, vcasingtol
+  ! if (accuracyestimate.gt.vcasingtol) then
+  !   FATAL( bnorml, .true., virtual casing accuracy is too low, increase vcNt and vcNz )
+  ! endif
+#ifdef COMPARECASING
+endif
+! When comparing results, both methods should always run
+if (.true.) then
+#else
+else ! if not Lvcgrid  
+#endif
+
+  do kk = 0, Nz-1 ; 
+   zeta = kk * pi2 / Nz
+
+   do jj = 0, Nt-1 ; 
+    teta = jj * pi2 / Nt ; jk = 1 + jj + kk*Nt
 
     globaljk = jk ! this is global; passed through to vcintegrand & casing;
 
@@ -156,35 +205,18 @@ subroutine bnorml( mn, Ntz, efmn, ofmn )
      FATAL( bnorml, .true., invalid Lparallel in parallelization loop )
     end select ! end of select case( Lparallel ) ; 09 Mar 17;
 
-    tetazeta(1:2) = (/ teta, zeta /) ! this is global; passed through to zetalow & zetaupp; 14 Apr 17;
-
-!#ifdef COMPARECASING
-!
-!    tetalow = tetazeta(1) - vcasingper * pi ; tetaupp = tetazeta(1) + vcasingper * pi ; absacc = vcasingtol
-!
-!    id01daf = 1 ; call D01DAF( tetalow, tetaupp, zetalow, zetaupp, vcintegrand, absacc, gBn, nvccalls, id01daf ) ! 04 May 17;
-!
-!    ijimag(jk) = gBn
-!
-!#endif
-
-    WCALL( bnorml, casing, ( teta, zeta, gBn, icasing ) ) ! tetazeta is global; 26 Apr 17;
-
+    WCALL( bnorml, casing, ( teta, zeta, gBn, icasing ) )
     ijreal(jk) = gBn
 
-!#ifdef COMPARECASING
-!    write(ounit,1000) myid, zeta, teta, ijreal(jk), ijimag(jk), ijreal(jk)-ijimag(jk)
-!#endif
-!
-!#ifdef CASING
-!    write(ounit,1000) myid, zeta, teta, ijreal(jk), ijimag(jk), ijreal(jk)-ijimag(jk)
-!#endif
-
-1000 format("bnorml : ", 10x ," : myid=",i3," : \z =",f6.3," ; \t =",f6.3," ; B . x_t x x_z =",2f22.15," ; ":"err =",es13.5," ;")
+#ifdef COMPARECASING
+  ! write(ounit,1000) myid, zeta, teta, ijreal(jk), ijimag(jk), ijreal(jk)-ijimag(jk)
+  write(ounit,1000) myid, zeta, teta, ijreal(jk), ijimag(jk), ijreal(jk)-ijimag(jk)
+  1000 format("bnorml : ", 10x ," : myid=",i3," : \z =",f6.3," ; \t =",f6.3," ; B . x_t x x_z =",2f22.15," ; ":"err =",es13.5," ;")
+#endif
 
    enddo ! end of do jj;
   enddo ! end of do kk;
-
+endif ! end of if (Lvcgrid  )
 !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
 
 1001 format("bnorml : ", 10x ," : "a1" : (t,z) = ("f8.4","f8.4" ) ; gBn=",f23.15," ; ":" error =",f23.15" ;")
@@ -206,11 +238,6 @@ subroutine bnorml( mn, Ntz, efmn, ofmn )
    case( 0 ) ! Lparallel = 0 ; 09 Mar 17;
 
     RlBCAST(ijreal(1+kk*Nt:Nt+kk*Nt),Nt,kkmodnp) ! plasma; 03 Apr 13;
-   !RlBCAST(ijimag(1+kk*Nt:Nt+kk*Nt),Nt,kkmodnp)
-
-   !RlBCAST(jireal(1+kk*Nt:Nt+kk*Nt),Nt,kkmodnp)
-   !RlBCAST(jiimag(1+kk*Nt:Nt+kk*Nt),Nt,kkmodnp)
-
    case( 1 ) ! Lparallel = 1 ; 09 Mar 17;
 
     do jj = 0, Nt-1
@@ -220,10 +247,6 @@ subroutine bnorml( mn, Ntz, efmn, ofmn )
      jkmodnp = modulo(jk-1,ncpu)
 
      RlBCAST(ijreal(jk),1,jkmodnp) ! plasma; 03 Apr 13;
-    !RlBCAST(ijimag(jk),1,jkmodnp)
-
-    !RlBCAST(jireal(jk),1,jkmodnp)
-    !RlBCAST(jiimag(jk),1,jkmodnp)
 
     enddo