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

Updated mp00ac.f90 to read and output gauge invariant helicity. #187

Open
wants to merge 10 commits into
base: master
Choose a base branch
from
Open

Updated mp00ac.f90 to read and output gauge invariant helicity. #187

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
330d40e
Updated mp00ac.f90 to read and output gauge invariant helicity. I bel…
KseniaAleynikova Mar 18, 2023
722cb1d
retrigger checks
smiet Oct 16, 2023
80a1122
retrigger checks
smiet Oct 16, 2023
d588faa
updated workflows for final re-trigger
smiet Oct 16, 2023
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
fd97ab3
Merge pull request #213 from missing-user/gauge_invariant_helicity
smiet Dec 1, 2024
2c7ae2b
Merge branch 'master' into gauge_invariant_helicity
smiet Dec 1, 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.
43 changes: 41 additions & 2 deletions src/mp00ac.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -117,7 +117,7 @@ subroutine mp00ac( Ndof, Xdof, Fdof, Ddof, Ldfjac, iflag ) ! argument list is fi

!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!

use constants, only : zero, half, one
use constants, only : zero, half, one, pi2

use numerical, only : small, machprec

Expand All @@ -127,7 +127,8 @@ subroutine mp00ac( Ndof, Xdof, Fdof, Ddof, Ldfjac, iflag ) ! argument list is fi
mu, helicity, iota, oita, curtor, curpol, Lrad, Ntor,&
!Lposdef, &
Lconstraint, mupftol, &
Lmatsolver, NiterGMRES, epsGMRES, LGMRESprec, epsILU
Lmatsolver, NiterGMRES, epsGMRES, LGMRESprec, epsILU,&
mpol

use cputiming, only : Tmp00ac

Expand Down Expand Up @@ -196,6 +197,11 @@ subroutine mp00ac( Ndof, Xdof, Fdof, Ddof, Ldfjac, iflag ) ! argument list is fi
REAL, allocatable :: wk(:)
INTEGER,allocatable :: jw(:), iperm(:)

REAL :: cheby(0:Lrad(ivol),0:2), zernike(0:Lrad(1),0:mpol,0:2)
REAL :: dH1, dH2
INTEGER :: ll, id


BEGIN(mp00ac)

!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
Expand Down Expand Up @@ -524,6 +530,39 @@ subroutine mp00ac( Ndof, Xdof, Fdof, Ddof, Ldfjac, iflag ) ! argument list is fi
lABintegral(lvol) = half * sum( solution(1:NN,0) * Ddotx(1:NN) )
endif


! if (GaugeInvariantH) then
dH1 = zero
dH2 = zero
do ll = 0, Lrad(lvol)
! ideriv=0 mn=0
id = Ate(lvol,0 , 1)%i(ll)
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 (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

!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
DALLOCATE( matrix )
DALLOCATE( rhs )
Expand Down