Merge pull request #213 from missing-user/gauge_invariant_helicity

[Gauge invariant Helicity] Out of bounds indexing corrected LGTM for merge into this branch. Test fails because branch is out-of-date with install changes due to numpy 2.0.
PrincetonUniversity · Dec 1, 2024 · fd97ab3 · fd97ab3
2 parents d588faa + 62ca200
commit fd97ab3
Show file tree

Hide file tree

Showing 6 changed files with 20 additions and 14 deletions.
diff --git a/ci/G1V03L2Fi/compare.h5 b/ci/G1V03L2Fi/compare.h5
diff --git a/ci/G2V32L1Fi/compare.h5 b/ci/G2V32L1Fi/compare.h5
diff --git a/ci/G3V08L3Fi/compare.h5 b/ci/G3V08L3Fi/compare.h5
diff --git a/ci/G3V08L3Fr/compare.h5 b/ci/G3V08L3Fr/compare.h5
diff --git a/ci/toroidal_freeboundary_vacuum/compare.h5 b/ci/toroidal_freeboundary_vacuum/compare.h5
diff --git a/src/mp00ac.f90 b/src/mp00ac.f90
@@ -532,28 +532,34 @@ subroutine mp00ac( Ndof, Xdof, Fdof, Ddof, Ldfjac, iflag ) ! argument list is fi
 
 
   ! if (GaugeInvariantH) then
+  dH1 = zero
+  dH2 = zero
   do ll = 0, Lrad(lvol)
     !             ideriv=0   mn=0
     id = Ate(lvol,0        , 1)%i(ll)
-    if (Lcoordinatesingularity) then
-      call get_zernike_d2(small, Lrad(lvol), mpol, zernike)
-      dH1 = dH1 + pi2 * solution(id,0) * zernike(ll,0,0)
-    else
-      call get_cheby_d2(-one+small, Lrad(lvol), cheby(0:Lrad(lvol),0:2))
-      dH1 = dH1 + pi2 * solution(id,0) * cheby(ll,0)
+    if (id/=0) then ! Indirection maps some elements to 0, to "discard" them during matrix construction. Ignore them here.
+      if (Lcoordinatesingularity) then
+        call get_zernike_d2(small, Lrad(lvol), mpol, zernike)
+        dH1 = dH1 + pi2 * solution(id,0) * zernike(ll,0,0)
+      else
+        call get_cheby_d2(-one+small, Lrad(lvol), cheby(0:Lrad(lvol),0:2))
+        dH1 = dH1 + pi2 * solution(id,0) * cheby(ll,0)
+      endif
     endif
-
+    
     !             ideriv=0   mn=0
     id = Aze(lvol,0        , 1)%i(ll)
-    if (Lcoordinatesingularity) then
-      call get_zernike_d2(one, Lrad(lvol), mpol, zernike)
-      dH2 = dH2 + pi2 * solution(id,0) * zernike(ll,0,0)
-    else
-      call get_cheby_d2(one, Lrad(lvol), cheby(0:Lrad(lvol),0:2))
-      dH2 = dH2 + pi2 * solution(id,0) * cheby(ll,0)
+    if (id/=0) then ! Indirection maps some elements to 0, to "discard" them during matrix construction. Ignore them here.
+      if (Lcoordinatesingularity) then
+        call get_zernike_d2(one, Lrad(lvol), mpol, zernike)
+        dH2 = dH2 + pi2 * solution(id,0) * zernike(ll,0,0)
+      else
+        call get_cheby_d2(one, Lrad(lvol), cheby(0:Lrad(lvol),0:2))
+        dH2 = dH2 + pi2 * solution(id,0) * cheby(ll,0)
+      endif
     endif
   enddo
-
+  
   lABintegral(lvol) = lABintegral(lvol) - dH1 - dH2
   ! endif