Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[Gauge invariant Helicity] Out of bounds indexing corrected #213

Merged
merged 4 commits into from
Dec 1, 2024
Merged

[Gauge invariant Helicity] Out of bounds indexing corrected #213

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
fd12aba
out of bounds indexing corrected
missing-user Nov 23, 2024
55d0c17
initialize dH2 and dH1 to zero for computing gauge
missing-user Nov 24, 2024
2a07c77
of course dH1 must be reset OUTSIDE of the loop....
missing-user Nov 26, 2024
62ca200
modified reference runs for modified helicity
missing-user Nov 26, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Binary file modified ci/G1V03L2Fi/compare.h5
View file
Open in desktop
Binary file not shown.
Binary file modified ci/G2V32L1Fi/compare.h5
View file
Open in desktop
Binary file not shown.
Binary file modified ci/G3V08L3Fi/compare.h5
View file
Open in desktop
Binary file not shown.
Binary file modified ci/G3V08L3Fr/compare.h5
View file
Open in desktop
Binary file not shown.
Binary file modified ci/toroidal_freeboundary_vacuum/compare.h5
View file
Open in desktop
Binary file not shown.
34 changes: 20 additions & 14 deletions src/mp00ac.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
Loading