From 0cb35b3e7c66856bbac5529a3aa4aae0d1e7aae0 Mon Sep 17 00:00:00 2001 From: ejpaul Date: Wed, 9 Sep 2020 07:22:46 -0400 Subject: [PATCH 01/20] Initial commit. --- brcast.f90 | 2 +- curent.f90 | 46 ++++++++++++++++++++++++++++++++++++---------- global.f90 | 6 +++--- ma02aa.f90 | 38 +++++++++++++++++++++++++++----------- mp00ac.f90 | 50 ++++++++++++++++++++++++++++++++++++++++++++++---- 5 files changed, 113 insertions(+), 29 deletions(-) diff --git a/brcast.f90 b/brcast.f90 index e3da6a32..ca6ee5ba 100644 --- a/brcast.f90 +++ b/brcast.f90 @@ -133,7 +133,7 @@ subroutine brcast( lvol ) !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! ! if( lvol.gt.Nvol .and. Lconstraint.eq.-1 .and. Wcurent ) then ! 27 Feb 17; - if( lvol.gt.Nvol .and. Wcurent ) then ! 27 Feb 17; + if( (lvol.gt.Nvol .or. Lconstraint .eq.-2) .and. Wcurent ) then ! 27 Feb 17; !write(ounit,'("brcast : " 10x " : myid="i3" ; broadcasting : curtor="es13.5" ; curpol="es13.5" ;")') myid, curtor, curpol RlBCAST( curtor, 1, llmodnp ) RlBCAST( curpol, 1, llmodnp ) diff --git a/curent.f90 b/curent.f90 index 77b0733c..0ca99377 100644 --- a/curent.f90 +++ b/curent.f90 @@ -58,13 +58,13 @@ subroutine curent( lvol, mn, Nt, Nz, iflag, ldItGp ) !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! - use constants, only : zero, one, two, pi2 + use constants, only : zero, one, two, pi2, half use numerical, only : use fileunits, only : ounit - use inputlist, only : Wmacros, Wcurent, Lrad + use inputlist, only : Wmacros, Wcurent, Lrad, Lconstraint use cputiming, only : Tcurent @@ -74,7 +74,7 @@ subroutine curent( lvol, mn, Nt, Nz, iflag, ldItGp ) sg, guvij, & Ntz, ijreal, ijimag, jireal, jiimag, & efmn, ofmn, cfmn, sfmn, evmn, odmn, comn, simn, & - Ate, Aze, Ato, Azo, TT + Ate, Aze, Ato, Azo, TT, Lcoordinatesingularity, regumm !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! @@ -84,9 +84,9 @@ subroutine curent( lvol, mn, Nt, Nz, iflag, ldItGp ) REAL , intent(out) :: ldItGp(0:1,-1:2) INTEGER :: innout, lvol, ideriv, ii, ll, Lcurvature, ifail - REAL :: lss - REAL :: lAte(1:mn,-1:2), lAze(1:mn,-1:2), lAto(1:mn,-1:2), lAzo(1:mn,-1:2) - REAL :: Bsupt(1:Nt*Nz,-1:2), Bsupz(1:Nt*Nz,-1:2) + REAL :: lss, mfactor + REAL :: lAte(1:mn,-1:2), lAze(1:mn,-1:2), lAto(1:mn,-1:2), lAzo(1:mn,-1:2), tAze(1:mn,-1:2), zAte(1:mn,-1:2) + REAL :: Bsupt(1:Nt*Nz,-1:2), Bsupz(1:Nt*Nz,-1:2), Bsups(1:Nt*Nz,-1:2), Bsups_2(1:Nt*Nz,-1:2) BEGIN(curent) @@ -108,6 +108,8 @@ subroutine curent( lvol, mn, Nt, Nz, iflag, ldItGp ) lAte(1:mn,-1:2) = zero ! radial derivatives of vector potential evaluated at interfaces; 20 Apr 13; lAze(1:mn,-1:2) = zero + zAte(1:mn,-1:2) = zero + tAze(1:mn,-1:2) = zero ! if( NOTstellsym ) then lAto(1:mn,-1:2) = zero ! this is used below and needs to be assigned a (trivial) value ; 26 Jan 16; lAzo(1:mn,-1:2) = zero ! this is used below and needs to be assigned a (trivial) value ; 26 Jan 16; @@ -115,6 +117,11 @@ subroutine curent( lvol, mn, Nt, Nz, iflag, ldItGp ) !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! + if (lconstraint .eq. -2) then + innout = 1.0 + lss = 1.0 + end if + do ideriv = -1, 2 ! labels derivative of magnetic field wrt enclosed fluxes; 20 Apr 13; if( iflag.eq. 1 .and. ideriv.ne.0 ) cycle ! derivatives of currents are not required; 20 Jun 14; @@ -124,11 +131,17 @@ subroutine curent( lvol, mn, Nt, Nz, iflag, ldItGp ) if( YESstellsym ) then do ii = 1, mn ! loop over Fourier harmonics; 20 Apr 13; + + if( Lcoordinatesingularity ) then ; mfactor = regumm(ii) * half ! include radial regularization factor near coordinate origin; 21 Apr 13; + else ; mfactor = zero + endif do ll = 0, Lrad(lvol) ! loop over Chebyshev polynomials; 20 Apr 13; - lAte(ii,ideriv) = lAte(ii,ideriv) + Ate(lvol,ideriv,ii)%s(ll) * TT(ll,innout,1) ! compute radial derivative of vector potential; 20 Apr 13; - lAze(ii,ideriv) = lAze(ii,ideriv) - Aze(lvol,ideriv,ii)%s(ll) * TT(ll,innout,1) ! note inclusion of sign factor ; 26 Jan 16; + lAte(ii,ideriv) = lAte(ii,ideriv) + Ate(lvol,ideriv,ii)%s(ll) * (TT(ll,innout,1) + mfactor) ! compute radial derivative of vector potential; 20 Apr 13; + lAze(ii,ideriv) = lAze(ii,ideriv) - Aze(lvol,ideriv,ii)%s(ll) * (TT(ll,innout,1) + mfactor) ! note inclusion of sign factor ; 26 Jan 16; + tAze(ii,ideriv) = tAze(ii,ideriv) - im(ii)*Aze(lvol,ideriv,ii)%s(ll) * TT(ll,innout,0) ! theta derivative of Az + ZAte(ii,ideriv) = zAte(ii,ideriv) - in(ii)*Ate(lvol,ideriv,ii)%s(ll) * TT(ll,innout,0) ! zeta derivative of At enddo ! end of do ll; @@ -155,6 +168,14 @@ subroutine curent( lvol, mn, Nt, Nz, iflag, ldItGp ) call invfft( mn, im(1:mn), in(1:mn), lAte(1:mn,ideriv), lAto(1:mn,ideriv), lAze(1:mn,ideriv), lAzo(1:mn,ideriv), & Nt, Nz, Bsupz(1:Ntz,ideriv), Bsupt(1:Ntz,ideriv) ) ! map to real space; + if (Lconstraint .eq. -2) then + call invfft( mn, im(1:mn), in(1:mn), lAto(1:mn,ideriv), tAze(1:mn,ideriv), lAzo(1:mn,ideriv), zAte(1:mn,ideriv), & + Nt, Nz, Bsups(1:Ntz,ideriv), Bsups_2(1:Ntz,ideriv)) + Bsups(1:Ntz,ideriv) = Bsups(1:Ntz,ideriv) + Bsups_2(1:ntz,ideriv) + else + Bsups = zero + end if + enddo ! end of do ideriv; 31 Jan 13; !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! @@ -175,8 +196,13 @@ subroutine curent( lvol, mn, Nt, Nz, iflag, ldItGp ) if( iflag.eq. 2 .and. ideriv.lt.0 ) cycle ! derivatives of currents wrt geometry is not required; 20 Jun 14; if( iflag.eq.-1 .and. ideriv.gt.0 ) cycle ! derivatives of currents wrt enclosed toroidal and enclosed poloidal flux are not required; 20 Jun 14; - ijreal(1:Ntz) = ( Bsupt(1:Ntz,ideriv) * guvij(1:Ntz,2,2,0) + Bsupz(1:Ntz,ideriv) * guvij(1:Ntz,2,3,0) ) / sg(1:Ntz,0) - ijimag(1:Ntz) = ( Bsupt(1:Ntz,ideriv) * guvij(1:Ntz,2,3,0) + Bsupz(1:Ntz,ideriv) * guvij(1:Ntz,3,3,0) ) / sg(1:Ntz,0) + if (Lconstraint .eq. -2) then + ijreal(1:Ntz) = (Bsups(1:Ntz,ideriv) * guvij(1:Ntz,2,1,0) + Bsupt(1:Ntz,ideriv) * guvij(1:Ntz,2,2,0) + Bsupz(1:Ntz,ideriv) * guvij(1:Ntz,2,3,0) ) / sg(1:Ntz,0) + ijimag(1:Ntz) = (Bsups(1:Ntz,ideriv) * guvij(1:Ntz,1,3,0) + Bsupt(1:Ntz,ideriv) * guvij(1:Ntz,2,3,0) + Bsupz(1:Ntz,ideriv) * guvij(1:Ntz,3,3,0) ) / sg(1:Ntz,0) + else + ijreal(1:Ntz) = (Bsupt(1:Ntz,ideriv) * guvij(1:Ntz,2,2,0) + Bsupz(1:Ntz,ideriv) * guvij(1:Ntz,2,3,0) ) / sg(1:Ntz,0) + ijimag(1:Ntz) = (Bsupt(1:Ntz,ideriv) * guvij(1:Ntz,2,3,0) + Bsupz(1:Ntz,ideriv) * guvij(1:Ntz,3,3,0) ) / sg(1:Ntz,0) + end if if( ideriv.eq.-1 ) then ! add derivatives of metrics with respect to interface geometry; 15 Sep 16; ijreal(1:Ntz) = ijreal(1:Ntz) + ( Bsupt(1:Ntz, 0) * guvij(1:Ntz,2,2,1) + Bsupz(1:Ntz, 0) * guvij(1:Ntz,2,3,1) ) / sg(1:Ntz,0) ijimag(1:Ntz) = ijimag(1:Ntz) + ( Bsupt(1:Ntz, 0) * guvij(1:Ntz,2,3,1) + Bsupz(1:Ntz, 0) * guvij(1:Ntz,3,3,1) ) / sg(1:Ntz,0) diff --git a/global.f90 b/global.f90 index 4a655b11..50b7be10 100644 --- a/global.f90 +++ b/global.f90 @@ -1557,7 +1557,7 @@ subroutine readin FATAL( readin, mupftol.le.zero, mupftol is too small ) FATAL( readin, abs(one+gamma).lt.vsmall, 1+gamma appears in denominator in dforce ) ! Please check this; SRH: 27 Feb 18; FATAL( readin, abs(one-gamma).lt.vsmall, 1-gamma appears in denominator in fu00aa ) ! Please check this; SRH: 27 Feb 18; - FATAL( readin, Lconstraint.lt.-1 .or. Lconstraint.gt.2, illegal Lconstraint ) + FATAL( readin, Lconstraint.lt.-2 .or. Lconstraint.gt.2, illegal Lconstraint ) FATAL( readin, Igeometry.eq.1 .and. rpol.lt.vsmall, poloidal extent of slab too small or negative ) FATAL( readin, Igeometry.eq.1 .and. rtor.lt.vsmall, toroidal extent of slab too small or negative ) @@ -2118,6 +2118,8 @@ subroutine readin iZbc(ii,Mvol) = zero endif + endif ! matches if( Lfreebound.eq.1 ) ; + iVns(ii ) = ( Vns( kk, mm) - Vns(-kk,-mm) ) * jj iBns(ii ) = ( Bns( kk, mm) - Bns(-kk,-mm) ) * jj if( NOTstellsym ) then @@ -2128,8 +2130,6 @@ subroutine readin iBnc(ii ) = zero endif - endif ! matches if( Lfreebound.eq.1 ) ; - endif ! end of if( mm.eq.0 .and. nn.eq.0 ) ; enddo ! end of do ii = 1, mn; diff --git a/ma02aa.f90 b/ma02aa.f90 index 3435b959..371c8f1b 100644 --- a/ma02aa.f90 +++ b/ma02aa.f90 @@ -415,14 +415,18 @@ subroutine ma02aa( lvol, NN ) if( LBlinear ) then ! assume Beltrami field is parameterized by helicity multiplier (and poloidal flux); - lastcpu = GETTIME - + lastcpu = GETTIME if( Lplasmaregion ) then Xdof(1:2) = xoffset + (/ mu(lvol), dpflux(lvol) /) ! initial guess for degrees of freedom; offset from zero so that relative error is small; + + if (Lconstraint .eq. -2) then + Xdof(1:1) = xoffset + (/ dtflux(lvol) /) + end if select case( Lconstraint ) + case( -2 ) ; ; Nxdof = 1 ! toroidal flux IS varied to match linking linking current; case( -1 ) ; ; Nxdof = 0 ! multiplier & poloidal flux NOT varied ; case( 0 ) ; ; Nxdof = 0 ! multiplier & poloidal flux NOT varied ; case( 1 ) ; if( Lcoordinatesingularity ) then ; Nxdof = 1 ! multiplier IS varied to match outer transform; @@ -468,14 +472,24 @@ subroutine ma02aa( lvol, NN ) WK(1:Ndof,1), WK(1:Ndof,2), WK(1:Ndof,3), WK(1:Ndof,4) ) ) if( Lplasmaregion ) then - - select case( ihybrj ) - case( 0: ) ; mu(lvol) = Xdof(1) - xoffset - ; ; dpflux(lvol) = Xdof(2) - xoffset - case( :-1) ; Xdof(1) = mu(lvol) + xoffset ! mu and dpflux have been updated in mp00ac; early termination; - ; ; Xdof(2) = dpflux(lvol) + xoffset ! mu and dpflux have been updated in mp00ac; early termination; - end select - + + + if (Lconstraint .eq. -2) then + if (ihybrj .le. -1) then + Xdof(1) = dtflux(lvol) + xoffset + elseif (ihybrj .ge. 0) then + dtflux(lvol) = Xdof(1) - xoffset + end if + print *,"dtflux(lvol): ", dtflux(lvol) + else + select case( ihybrj ) + case( 0: ) ; mu(lvol) = Xdof(1) - xoffset + ; ; dpflux(lvol) = Xdof(2) - xoffset + case( :-1) ; Xdof(1) = mu(lvol) + xoffset ! mu and dpflux have been updated in mp00ac; early termination; + ; ; Xdof(2) = dpflux(lvol) + xoffset ! mu and dpflux have been updated in mp00ac; early termination; + end select + end if + else ! Lvacuumregion; select case( ihybrj ) @@ -491,7 +505,7 @@ subroutine ma02aa( lvol, NN ) !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! - if( Lconstraint.eq.1 .or. ( Lvacuumregion .and. Lconstraint.eq.0 ) ) then + if( Lconstraint.eq.1 .or. ( Lvacuumregion .and. Lconstraint.eq.0 ) .or. Lconstraint.eq.-2) then iflag = 2 ; Ldfjac = Ndof ! call mp00ac: tr00ab/curent to ensure the derivatives of B, transform, currents, wrt mu/dtflux & dpflux are calculated; @@ -527,6 +541,8 @@ subroutine ma02aa( lvol, NN ) end select endif ! end of if( LBlinear ) then; + + print *,"End of ma02aa.f90- dtflux(lvol): ", dtflux(lvol) !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!- diff --git a/mp00ac.f90 b/mp00ac.f90 index 3b8b702a..07266e7d 100644 --- a/mp00ac.f90 +++ b/mp00ac.f90 @@ -198,6 +198,17 @@ subroutine mp00ac( Ndof, Xdof, Fdof, Ddof, Ldfjac, iflag ) ! argument list is fi if( Ndof.eq.2 ) then ; dpf = Xdof(2) - xoffset else ; dpf = dpflux(lvol) endif + + + ! only poloidal flux is modified + if (Lconstraint .eq. -2) then + lmu = mu(lvol) + dpf = dpflux(lvol) + dtf = Xdof(1) - xoffset + print *,"dtf: ", dtf + print *,"dpf: ", dpf + print *,"lmu: ", lmu + end if else ! Lvacuumregion; @@ -263,12 +274,22 @@ subroutine mp00ac( Ndof, Xdof, Fdof, Ddof, Ldfjac, iflag ) ! argument list is fi ! !;case( 1 ) ; rhs(1:NN,1) = - matmul( - one * dME(1:NN,1:2), dpsi(1:2) ) - matmul( - one * dMD(1:NN,1:NN), solution(1:NN,0) ) ! !; ; ; rhs(1:NN,2) = - matmul( dMB(1:NN,1:2) - lmu * dME(1:NN,1:2), ppsi(1:2) ) ! !;end select - + ;select case( ideriv ) - ;case( 0 ) ; rhs(1:NN,0) = - matmul( dMB(1:NN,1:2) , dpsi(1:2) ) + ;case( 0 ) ; rhs(1:NN,0) = - dMG(1:NN) - matmul( dMB(1:NN,1:2) , dpsi(1:2) ) ;case( 1 ) ; rhs(1:NN,1) = - matmul( - one * dMD(1:NN,1:NN), solution(1:NN,0) ) ; ; ; rhs(1:NN,2) = - matmul( dMB(1:NN,1:2) , ppsi(1:2) ) ;end select + + if (Lconstraint .eq. -2) then + select case(ideriv) + case (0) + rhs(1:NN,0) = - dMG(1:NN) - matmul( dMB(1:NN,1:2), dpsi(1:2) ) + case (1) + rhs(1:NN,1) = - matmul( dMB(1:NN,1:2), tpsi(1:2) ) + ! rhs(1:NN,2) = - matmul( - one * dMD(1:NN,1:NN),solution(1:NN,0) ) + end select + end if else ! .not.Lcoordinatesingularity; @@ -409,7 +430,22 @@ subroutine mp00ac( Ndof, Xdof, Fdof, Ddof, Ldfjac, iflag ) ! argument list is fi !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! select case( Lconstraint ) - + + case( -2 ) ! Lconstraint=-2 + + WCALL( mp00ac, curent,( lvol, mn, Nt, Nz, iflag, dItGpdxtp(0:1,-1:2,lvol)) ) + + ! Constraint on curpol only + if( iflag.eq.1 ) Fdof(1 ) = dItGpdxtp(1,0,lvol) - curpol + ! Derivative w.r.t. toroidal flux only + if( iflag.eq.2 ) Ddof(1,1) = dItGpdxtp(1,1,lvol) + if (iflag .eq. 1) then + print *,"curtor (computed): ", dItGpdxtp(0,0,lvol) + print *,"curpol (computed): ", dItGpdxtp(1,0,lvol) + print *,"curpol (requested): ", curpol + print *,"Fdof: ", Fdof(1) + end if + case( -1 ) ! Lconstraint=-1; if( Lplasmaregion ) then @@ -531,13 +567,19 @@ subroutine mp00ac( Ndof, Xdof, Fdof, Ddof, Ldfjac, iflag ) ! argument list is fi if( sum( abs( Fdof(1:Ndof) ) ) / Ndof .lt. mupftol ) then ! satisfactory; - if ( Lplasmaregion ) then ; mu(lvol) = lmu ; ; dpflux(lvol) = dpf + if ( Lplasmaregion .and. Lconstraint .ne. 2) then ; mu(lvol) = lmu ; ; dpflux(lvol) = dpf #ifdef FORCEFREEVACUUM else ; mu(lvol) = lmu ; dtflux(lvol) = dtf ; dpflux(lvol) = dpf #else else ; mu(lvol) = zero ; dtflux(lvol) = dtf ; dpflux(lvol) = dpf #endif endif + + if ( Lconstraint .eq. -2) then + !mu(lvol) = lmu + dtflux(lvol) = dtf + print *,"mu(lvol): ", mu(lvol) + end if iflag = -2 ! return "acceptance" flag through to ma02aa via ifail; early termination; From b5d0e8954462eb2e0159211a3f8d12a1dc17cc57 Mon Sep 17 00:00:00 2001 From: Richard Nies Date: Mon, 17 May 2021 17:39:50 -0400 Subject: [PATCH 02/20] Changes in SPEC + further merge conflict resolutions. --- .../py_spec/input/spec_namelist.py | 5 +- .../pythontools/py_spec/output/_processing.py | 71 +++++++++++++++++-- Utilities/pythontools/py_spec/output/spec.py | 2 + curent.f90 | 11 ++- global.f90 | 26 +++---- mp00ac.f90 | 1 + 6 files changed, 93 insertions(+), 23 deletions(-) diff --git a/Utilities/pythontools/py_spec/input/spec_namelist.py b/Utilities/pythontools/py_spec/input/spec_namelist.py index 344113cc..b9a9ec9f 100644 --- a/Utilities/pythontools/py_spec/input/spec_namelist.py +++ b/Utilities/pythontools/py_spec/input/spec_namelist.py @@ -216,7 +216,10 @@ def run(self, spec_command="./xspec", filename="spec.sp", force=False, quiet=Fal if run_result.returncode == 0: # the run is successful if not quiet: print("SPEC runs successfully.") - return SPECout(filename + ".h5") + try: + return SPECout(filename + ".h5") + except: + return SPECout(filename[:-3] + ".h5") else: print("SPEC runs unsuccessfully, check terminal output.") return None diff --git a/Utilities/pythontools/py_spec/output/_processing.py b/Utilities/pythontools/py_spec/output/_processing.py index 3247f6d0..99ba3e62 100644 --- a/Utilities/pythontools/py_spec/output/_processing.py +++ b/Utilities/pythontools/py_spec/output/_processing.py @@ -165,7 +165,6 @@ def metric( _, _, _, g = get_grid_and_jacobian_and_metric(self, lvol, sarr, tarr, zarr) return g - def get_B( self, lvol=0, @@ -174,6 +173,28 @@ def get_B( tarr=np.linspace(0, 0, 1), zarr=np.linspace(0, 0, 1), ): + Bcontra_and_s_der = get_B_and_s_der(self, lvol, jacobian, sarr, tarr, zarr) + return Bcontra_and_s_der[0:3] + +def get_s_der_B( + self, + lvol=0, + jacobian=None, + sarr=np.linspace(0, 0, 1), + tarr=np.linspace(0, 0, 1), + zarr=np.linspace(0, 0, 1), +): + Bcontra_and_s_der = get_B_and_s_der(self, lvol, jacobian, sarr, tarr, zarr) + return Bcontra_and_s_der[3:6] + +def get_B_and_s_der( + self, + lvol=0, + jacobian=None, + sarr=np.linspace(0, 0, 1), + tarr=np.linspace(0, 0, 1), + zarr=np.linspace(0, 0, 1), +): if jacobian is None: R, Z, jacobian, g = get_grid_and_jacobian_and_metric( @@ -207,7 +228,7 @@ def get_B( # Zernike polynomial being used from ._processing import _get_zernike - zernike, dzernike = _get_zernike(sarr, Lrad, Mpol) + zernike, dzernike, d2zernike = _get_zernike(sarr, Lrad, Mpol) c = ( im[nax, :, nax, nax] * Azo.T[:, :, nax, nax] @@ -220,12 +241,14 @@ def get_B( ] Bs = np.sum(zernike[:, :, im, None, None] * c[None, :, :, :, :], axis=(1, 2)) + dsBs = np.sum(dzernike[:, :, im, None, None] * c[None, :, :, :, :], axis=(1, 2)) c1 = ( Aze.T[:, :, nax, nax] * cosa[nax, :, :, :] + Azo.T[:, :, nax, nax] * sina[nax, :, :, :] ) Bt = -np.sum(dzernike[:, :, im, None, None] * c1[None, :, :, :, :], axis=(1, 2)) + dsBt = -np.sum(d2zernike[:, :, im, None, None] * c1[None, :, :, :, :], axis=(1, 2)) c2 = ( Ate.T[:, :, nax, nax] * cosa[nax, :, :, :] @@ -233,6 +256,7 @@ def get_B( ) Bz = np.sum(dzernike[:, :, im, None, None] * c2[None, :, :, :, :], axis=(1, 2)) + dsBz = np.sum(dzernike[:, :, im, None, None] * c2[None, :, :, :, :], axis=(1, 2)) else: # Chebyshev being used @@ -262,6 +286,7 @@ def get_B( ] Bs = np.rollaxis(np.sum(Cheb.chebval(sarr, c), axis=0), 2) + dsBs = np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c)), axis=0), 2) c1 = ( Aze.T[:, :, nax, nax] * cosa[nax, :, :, :] @@ -274,6 +299,13 @@ def get_B( ), 2, ) + dsBt = -np.rollaxis( + np.sum( + Cheb.chebval(sarr, Cheb.chebder(c1,m=2)), + axis=0, + ), + 2, + ) c2 = ( Ate.T[:, :, nax, nax] * cosa[nax, :, :, :] @@ -286,9 +318,16 @@ def get_B( ), 2, ) + dsBz = np.rollaxis( + np.sum( + Cheb.chebval(sarr, Cheb.chebder(c2,m=2)), + axis=0, + ), + 2, + ) - Bcontrav = np.array([Bs, Bt, Bz]) / jacobian - return Bcontrav + Bcontrav_and_s_der = np.array([Bs, Bt, Bz, dsBs, dsBt, dsBz]) / jacobian + return Bcontrav_and_s_der # Bcontrav = get_B(s,lvol=lvol,jacobian=jacobian,sarr=sarr,tarr=tarr,zarr=zarr) @@ -307,7 +346,7 @@ def get_B_covariant(self, Bcontrav, g): def _get_zernike(sarr, lrad, mpol): """ - Get the value of the zernike polynomials and their derivatives + Get the value of the zernike polynomials, their first and second derivatives Adapted from basefn.f90 """ @@ -319,17 +358,22 @@ def _get_zernike(sarr, lrad, mpol): zernike = np.zeros([ns, lrad + 1, mpol + 1], dtype=np.float64) dzernike = np.zeros_like(zernike) + d2zernike = np.zeros_like(zernike) for m in range(mpol + 1): if lrad >= m: zernike[:, m, m] = rm dzernike[:, m, m] = m * rm1 + d2zernike[:, m, m] = m * (m-1) * rm1/r if lrad >= m + 2: zernike[:, m + 2, m] = float(m + 2) * rm * r ** 2 - float(m + 1) * rm dzernike[:, m + 2, m] = ( float(m + 2) ** 2 * rm * r - float((m + 1) * m) * rm1 ) + d2zernike[:, m + 2, m] = ( + float(m + 2) ** 2 * (m + 1) * rm - float((m + 1) * m * (m - 1)) * rm1 / r + ) for n in range(m + 4, lrad + 1, 2): factor1 = float(n) / float(n ** 2 - m ** 2) @@ -348,6 +392,13 @@ def _get_zernike(sarr, lrad, mpol): + (factor2 * r ** 2 - factor3) * dzernike[:, n - 2, m] - factor4 * dzernike[:, n - 4, m] ) + d2zernike[:, n, m] = factor1 * ( + 2 * factor2 * zernike[:, n - 2, m] + + 2 * factor2 * r * dzernike[:, n - 2, m] + + (factor2 * 2 * r ) * dzernike[:, n - 2, m] + + (factor2 * r ** 2 - factor3) * d2zernike[:, n - 2, m] + - factor4 * d2zernike[:, n - 4, m] + ) rm1 = rm rm = rm * r @@ -363,12 +414,18 @@ def _get_zernike(sarr, lrad, mpol): dzernike[:, n, 1] = dzernike[:, n, 1] - (-1) ** ((n - 1) / 2) * float( (n + 1) / 2 ) + d2zernike[:, n, 1] = d2zernike[:, n, 1] - (-1) ** ((n - 1) / 2) * float( + (n + 1) / 2 / r + ) + for m in range(mpol + 1): for n in range(m, lrad + 1, 2): zernike[:, n, m] = zernike[:, n, m] / float(n + 1) dzernike[:, n, m] = dzernike[:, n, m] / float(n + 1) + d2zernike[:, n, m] = d2zernike[:, n, m] / float(n + 1) - dzernike = dzernike * 0.5 # to account for the factor of half in sbar = (1+s)/2 + dzernike = dzernike * 0.5 # to account for the factor of half in sbar = (1+s)/2 + d2zernike = d2zernike * 0.25 # to account for the factor of half in sbar = (1+s)/2 - return zernike, dzernike \ No newline at end of file + return zernike, dzernike, d2zernike diff --git a/Utilities/pythontools/py_spec/output/spec.py b/Utilities/pythontools/py_spec/output/spec.py index 7eaf7fc2..3a2cb3ac 100644 --- a/Utilities/pythontools/py_spec/output/spec.py +++ b/Utilities/pythontools/py_spec/output/spec.py @@ -35,6 +35,8 @@ class SPECout: jacobian, metric, get_B, + get_s_der_B, + get_B_and_s_der, get_modB, get_B_covariant ) diff --git a/curent.f90 b/curent.f90 index b630f46a..78423989 100644 --- a/curent.f90 +++ b/curent.f90 @@ -125,13 +125,18 @@ subroutine curent( lvol, mn, Nt, Nz, iflag, ldItGp ) call build_vector_potential(lvol, innout, ideriv, 0) do ii = 1, mn ! loop over Fourier harmonics; 17 May 21; - efmn(ii) = -im(ii)*efmn(ii) ! theta derivative of Aze - cfmn(ii) = -in(ii)*cfmn(ii) ! zeta derivative of Ate + efmn(ii) = -in(ii)*efmn(ii) ! zeta derivative of Ate + ofmn(ii) = -in(ii)*ofmn(ii) ! zeta derivative of Ato + cfmn(ii) = -im(ii)*cfmn(ii) ! theta derivative of Aze + sfmn(ii) = -im(ii)*sfmn(ii) ! theta derivative of Azo enddo ! end of do ii; - call invfft( mn, im(1:mn), in(1:mn), ofmn(1:mn), efmn(1:mn), -sfmn(1:mn), cfmn(1:mn), & + call invfft( mn, im(1:mn), in(1:mn), ofmn(1:mn), cfmn(1:mn), sfmn(1:mn), efmn(1:mn), & Nt, Nz, Bsups(1:Ntz,ideriv), Bsups_2(1:Ntz,ideriv)) Bsups(1:Ntz,ideriv) = Bsups(1:Ntz,ideriv) + Bsups_2(1:ntz,ideriv) + print *,"Bsups: ", Bsups(1,ideriv) + print *,"Bsupt: ", Bsupt(1,ideriv) + print *,"Bsupz: ", Bsupz(1,ideriv) else Bsups = zero end if diff --git a/global.f90 b/global.f90 index 67ab8081..ca6cc215 100644 --- a/global.f90 +++ b/global.f90 @@ -829,21 +829,23 @@ subroutine build_vector_potential(lvol, iocons, aderiv, tderiv) if( Lcoordinatesingularity ) then mi = im(ii) - do ll = mi, Lrad(lvol),2 ! loop over Zernike polynomials; Lrad is the radial resolution; 01 Jul 19; if( tderiv .eq. 1) then - ; ; efmn(ii) = efmn(ii) + Ate(lvol,aderiv,ii)%s(ll) * RTT(ll,mi,iocons,1) * half - ; ; cfmn(ii) = cfmn(ii) + Aze(lvol,aderiv,ii)%s(ll) * RTT(ll,mi,iocons,1) * half - if( NOTstellsym ) then ; ofmn(ii) = ofmn(ii) + Ato(lvol,aderiv,ii)%s(ll) * RTT(ll,mi,iocons,1) * half - ; ; sfmn(ii) = sfmn(ii) + Azo(lvol,aderiv,ii)%s(ll) * RTT(ll,mi,iocons,1) * half - endif + do ll = mi, Lrad(lvol),2 ! loop over Zernike polynomials; Lrad is the radial resolution; 01 Jul 19; + ; ; efmn(ii) = efmn(ii) + Ate(lvol,aderiv,ii)%s(ll) * RTT(ll,mi,iocons,1) * half + ; ; cfmn(ii) = cfmn(ii) + Aze(lvol,aderiv,ii)%s(ll) * RTT(ll,mi,iocons,1) * half + if( NOTstellsym ) then ; ofmn(ii) = ofmn(ii) + Ato(lvol,aderiv,ii)%s(ll) * RTT(ll,mi,iocons,1) * half + ; ; sfmn(ii) = sfmn(ii) + Azo(lvol,aderiv,ii)%s(ll) * RTT(ll,mi,iocons,1) * half + endif + enddo ! end of do ll; 20 Feb 13; else - ; ; efmn(ii) = efmn(ii) + Ate(lvol,aderiv,ii)%s(ll) * RTT(ll,mi,iocons,0) - ; ; cfmn(ii) = cfmn(ii) + Aze(lvol,aderiv,ii)%s(ll) * RTT(ll,mi,iocons,0) - if( NOTstellsym ) then ; ofmn(ii) = ofmn(ii) + Ato(lvol,aderiv,ii)%s(ll) * RTT(ll,mi,iocons,0) - ; ; sfmn(ii) = sfmn(ii) + Azo(lvol,aderiv,ii)%s(ll) * RTT(ll,mi,iocons,0) - endif + do ll = 0, Lrad(lvol) ! loop over Zernike polynomials; Lrad is the radial resolution; 01 Jul 19; + ; ; efmn(ii) = efmn(ii) + Ate(lvol,aderiv,ii)%s(ll) * RTT(ll,mi,iocons,0) + ; ; cfmn(ii) = cfmn(ii) + Aze(lvol,aderiv,ii)%s(ll) * RTT(ll,mi,iocons,0) + if( NOTstellsym ) then ; ofmn(ii) = ofmn(ii) + Ato(lvol,aderiv,ii)%s(ll) * RTT(ll,mi,iocons,0) + ; ; sfmn(ii) = sfmn(ii) + Azo(lvol,aderiv,ii)%s(ll) * RTT(ll,mi,iocons,0) + endif + enddo endif - enddo ! end of do ll; 20 Feb 13; else do ll = 0, Lrad(lvol) ! loop over Chebyshev polynomials; Lrad is the radial resolution; ; ; efmn(ii) = efmn(ii) + Ate(lvol,aderiv,ii)%s(ll) * TT(ll,iocons,tderiv) ! aderiv labels deriv. wrt mu, pflux; diff --git a/mp00ac.f90 b/mp00ac.f90 index 5a6cfc44..dbf1ea03 100644 --- a/mp00ac.f90 +++ b/mp00ac.f90 @@ -350,6 +350,7 @@ subroutine mp00ac( Ndof, Xdof, Fdof, Ddof, Ldfjac, iflag ) ! argument list is fi else ! Matrix free version rhs(1:NN,1) = Ddotx(1:NN) ! rhs(1:NN,2) = - matmul( - one * dMD(1:NN,1:NN),solution(1:NN,0) ) + end if end select end if From 9352b8c9ca0386aef92edb062be88d25256bd95f Mon Sep 17 00:00:00 2001 From: Richard Nies Date: Mon, 24 May 2021 13:03:43 -0400 Subject: [PATCH 03/20] Normal field B.C. now correctly read in. --- .../pythontools/py_spec/output/_processing.py | 94 ++--- curent.f90 | 10 +- global.f90 | 342 +----------------- preset.f90 | 5 +- 4 files changed, 44 insertions(+), 407 deletions(-) diff --git a/Utilities/pythontools/py_spec/output/_processing.py b/Utilities/pythontools/py_spec/output/_processing.py index 99ba3e62..72c9998f 100644 --- a/Utilities/pythontools/py_spec/output/_processing.py +++ b/Utilities/pythontools/py_spec/output/_processing.py @@ -169,7 +169,7 @@ def get_B( self, lvol=0, jacobian=None, - sarr=np.linspace(0, 0, 1), + sarr=np.linspace(1, 1, 1), tarr=np.linspace(0, 0, 1), zarr=np.linspace(0, 0, 1), ): @@ -180,7 +180,7 @@ def get_s_der_B( self, lvol=0, jacobian=None, - sarr=np.linspace(0, 0, 1), + sarr=np.linspace(1, 1, 1), tarr=np.linspace(0, 0, 1), zarr=np.linspace(0, 0, 1), ): @@ -191,7 +191,7 @@ def get_B_and_s_der( self, lvol=0, jacobian=None, - sarr=np.linspace(0, 0, 1), + sarr=np.linspace(1, 1, 1), tarr=np.linspace(0, 0, 1), zarr=np.linspace(0, 0, 1), ): @@ -202,6 +202,7 @@ def get_B_and_s_der( ) # Lrad = s.input.physics.Lrad[lvol] + Ate = self.vector_potential.Ate[lvol] Aze = self.vector_potential.Aze[lvol] Ato = self.vector_potential.Ato[lvol] @@ -228,7 +229,7 @@ def get_B_and_s_der( # Zernike polynomial being used from ._processing import _get_zernike - zernike, dzernike, d2zernike = _get_zernike(sarr, Lrad, Mpol) + zernike, dzernike, d2zernike = _get_zernike(sbar, Lrad, Mpol) c = ( im[nax, :, nax, nax] * Azo.T[:, :, nax, nax] @@ -240,23 +241,23 @@ def get_B_and_s_der( nax, :, :, : ] - Bs = np.sum(zernike[:, :, im, None, None] * c[None, :, :, :, :], axis=(1, 2)) + Bs = np.sum( zernike[:, :, im, None, None] * c[None, :, :, :, :], axis=(1, 2)) dsBs = np.sum(dzernike[:, :, im, None, None] * c[None, :, :, :, :], axis=(1, 2)) c1 = ( - Aze.T[:, :, nax, nax] * cosa[nax, :, :, :] + Aze.T[:, :, nax, nax] * cosa[nax, :, :, :] + Azo.T[:, :, nax, nax] * sina[nax, :, :, :] ) - Bt = -np.sum(dzernike[:, :, im, None, None] * c1[None, :, :, :, :], axis=(1, 2)) + Bt = -np.sum( dzernike[:, :, im, None, None] * c1[None, :, :, :, :], axis=(1, 2)) dsBt = -np.sum(d2zernike[:, :, im, None, None] * c1[None, :, :, :, :], axis=(1, 2)) c2 = ( - Ate.T[:, :, nax, nax] * cosa[nax, :, :, :] + Ate.T[:, :, nax, nax] * cosa[nax, :, :, :] + Ato.T[:, :, nax, nax] * sina[nax, :, :, :] ) - Bz = np.sum(dzernike[:, :, im, None, None] * c2[None, :, :, :, :], axis=(1, 2)) - dsBz = np.sum(dzernike[:, :, im, None, None] * c2[None, :, :, :, :], axis=(1, 2)) + Bz = np.sum( dzernike[:, :, im, None, None] * c2[None, :, :, :, :], axis=(1, 2)) + dsBz = np.sum(d2zernike[:, :, im, None, None] * c2[None, :, :, :, :], axis=(1, 2)) else: # Chebyshev being used @@ -285,46 +286,23 @@ def get_B_and_s_der( nax, :, :, : ] - Bs = np.rollaxis(np.sum(Cheb.chebval(sarr, c), axis=0), 2) + Bs = np.rollaxis(np.sum(Cheb.chebval(sarr, c), axis=0), 2) dsBs = np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c)), axis=0), 2) c1 = ( - Aze.T[:, :, nax, nax] * cosa[nax, :, :, :] + Aze.T[:, :, nax, nax] * cosa[nax, :, :, :] + Azo.T[:, :, nax, nax] * sina[nax, :, :, :] ) - Bt = -np.rollaxis( - np.sum( - Cheb.chebval(sarr, Cheb.chebder(c1)), - axis=0, - ), - 2, - ) - dsBt = -np.rollaxis( - np.sum( - Cheb.chebval(sarr, Cheb.chebder(c1,m=2)), - axis=0, - ), - 2, - ) + Bt = -np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c1)), axis=0),2) + dsBt = -np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c1,m=2)), axis=0),2) c2 = ( - Ate.T[:, :, nax, nax] * cosa[nax, :, :, :] + Ate.T[:, :, nax, nax] * cosa[nax, :, :, :] + Ato.T[:, :, nax, nax] * sina[nax, :, :, :] ) - Bz = np.rollaxis( - np.sum( - Cheb.chebval(sarr, Cheb.chebder(c2)), - axis=0, - ), - 2, - ) - dsBz = np.rollaxis( - np.sum( - Cheb.chebval(sarr, Cheb.chebder(c2,m=2)), - axis=0, - ), - 2, - ) + Bz = np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c2)), axis=0),2) + dsBz = np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c2,m=2)), axis=0),2) + Bcontrav_and_s_der = np.array([Bs, Bt, Bz, dsBs, dsBt, dsBz]) / jacobian return Bcontrav_and_s_der @@ -355,6 +333,7 @@ def _get_zernike(sarr, lrad, mpol): r = (sarr + 1.0) / 2 rm = np.ones_like(r) # r to the power of m'th rm1 = np.zeros_like(r) # r to the power of m-1'th + rm2 = np.zeros_like(r) # r to the power of m-2'th zernike = np.zeros([ns, lrad + 1, mpol + 1], dtype=np.float64) dzernike = np.zeros_like(zernike) @@ -362,25 +341,19 @@ def _get_zernike(sarr, lrad, mpol): for m in range(mpol + 1): if lrad >= m: - zernike[:, m, m] = rm - dzernike[:, m, m] = m * rm1 - d2zernike[:, m, m] = m * (m-1) * rm1/r + zernike[:, m, m] = rm + dzernike[:, m, m] = m * rm1 + d2zernike[:, m, m] = m * (m-1) * rm2 if lrad >= m + 2: - zernike[:, m + 2, m] = float(m + 2) * rm * r ** 2 - float(m + 1) * rm - dzernike[:, m + 2, m] = ( - float(m + 2) ** 2 * rm * r - float((m + 1) * m) * rm1 - ) - d2zernike[:, m + 2, m] = ( - float(m + 2) ** 2 * (m + 1) * rm - float((m + 1) * m * (m - 1)) * rm1 / r - ) + zernike[:, m + 2, m] = float(m + 2) * rm * r**2 - float(m + 1) * rm + dzernike[:, m + 2, m] = (float(m + 2)**2 * rm * r - float((m + 1) * m) * rm1) + d2zernike[:, m + 2, m] = (float((m + 2)**2 * (m + 1)) * rm - float((m + 1) * m * (m - 1)) * rm2) for n in range(m + 4, lrad + 1, 2): factor1 = float(n) / float(n ** 2 - m ** 2) factor2 = float(4 * (n - 1)) - factor3 = float((n - 2 + m) ** 2) / float(n - 2) + float( - (n - m) ** 2 - ) / float(n) + factor3 = float((n - 2 + m) ** 2) / float(n - 2) + float((n - m) ** 2) / float(n) factor4 = float((n - 2) ** 2 - m ** 2) / float(n - 2) zernike[:, n, m] = factor1 * ( @@ -394,14 +367,14 @@ def _get_zernike(sarr, lrad, mpol): ) d2zernike[:, n, m] = factor1 * ( 2 * factor2 * zernike[:, n - 2, m] - + 2 * factor2 * r * dzernike[:, n - 2, m] - + (factor2 * 2 * r ) * dzernike[:, n - 2, m] + + 4 * factor2 * r * dzernike[:, n - 2, m] + (factor2 * r ** 2 - factor3) * d2zernike[:, n - 2, m] - factor4 * d2zernike[:, n - 4, m] ) + rm2 = rm1 rm1 = rm - rm = rm * r + rm = rm * r for n in range(2, lrad + 1, 2): zernike[:, n, 0] = zernike[:, n, 0] - (-1) ** (n / 2) @@ -414,15 +387,12 @@ def _get_zernike(sarr, lrad, mpol): dzernike[:, n, 1] = dzernike[:, n, 1] - (-1) ** ((n - 1) / 2) * float( (n + 1) / 2 ) - d2zernike[:, n, 1] = d2zernike[:, n, 1] - (-1) ** ((n - 1) / 2) * float( - (n + 1) / 2 / r - ) for m in range(mpol + 1): for n in range(m, lrad + 1, 2): - zernike[:, n, m] = zernike[:, n, m] / float(n + 1) - dzernike[:, n, m] = dzernike[:, n, m] / float(n + 1) + zernike[:, n, m] = zernike[:, n, m] / float(n + 1) + dzernike[:, n, m] = dzernike[:, n, m] / float(n + 1) d2zernike[:, n, m] = d2zernike[:, n, m] / float(n + 1) dzernike = dzernike * 0.5 # to account for the factor of half in sbar = (1+s)/2 diff --git a/curent.f90 b/curent.f90 index 78423989..05f29074 100644 --- a/curent.f90 +++ b/curent.f90 @@ -105,7 +105,7 @@ subroutine curent( lvol, mn, Nt, Nz, iflag, ldItGp ) !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! - if (lconstraint .eq. -2) then + if (Lconstraint .eq. -2) then innout = 1.0 lss = 1.0 end if @@ -126,17 +126,15 @@ subroutine curent( lvol, mn, Nt, Nz, iflag, ldItGp ) do ii = 1, mn ! loop over Fourier harmonics; 17 May 21; efmn(ii) = -in(ii)*efmn(ii) ! zeta derivative of Ate - ofmn(ii) = -in(ii)*ofmn(ii) ! zeta derivative of Ato + ofmn(ii) = in(ii)*ofmn(ii) ! zeta derivative of Ato cfmn(ii) = -im(ii)*cfmn(ii) ! theta derivative of Aze - sfmn(ii) = -im(ii)*sfmn(ii) ! theta derivative of Azo + sfmn(ii) = im(ii)*sfmn(ii) ! theta derivative of Azo enddo ! end of do ii; call invfft( mn, im(1:mn), in(1:mn), ofmn(1:mn), cfmn(1:mn), sfmn(1:mn), efmn(1:mn), & Nt, Nz, Bsups(1:Ntz,ideriv), Bsups_2(1:Ntz,ideriv)) Bsups(1:Ntz,ideriv) = Bsups(1:Ntz,ideriv) + Bsups_2(1:ntz,ideriv) - print *,"Bsups: ", Bsups(1,ideriv) - print *,"Bsupt: ", Bsupt(1,ideriv) - print *,"Bsupz: ", Bsupz(1,ideriv) + else Bsups = zero end if diff --git a/global.f90 b/global.f90 index ca6cc215..10f5ed2d 100644 --- a/global.f90 +++ b/global.f90 @@ -829,23 +829,21 @@ subroutine build_vector_potential(lvol, iocons, aderiv, tderiv) if( Lcoordinatesingularity ) then mi = im(ii) - if( tderiv .eq. 1) then - do ll = mi, Lrad(lvol),2 ! loop over Zernike polynomials; Lrad is the radial resolution; 01 Jul 19; + do ll = mi, Lrad(lvol),2 ! loop over Zernike polynomials; Lrad is the radial resolution; 01 Jul 19; + if( tderiv .eq. 1) then ; ; efmn(ii) = efmn(ii) + Ate(lvol,aderiv,ii)%s(ll) * RTT(ll,mi,iocons,1) * half ; ; cfmn(ii) = cfmn(ii) + Aze(lvol,aderiv,ii)%s(ll) * RTT(ll,mi,iocons,1) * half if( NOTstellsym ) then ; ofmn(ii) = ofmn(ii) + Ato(lvol,aderiv,ii)%s(ll) * RTT(ll,mi,iocons,1) * half ; ; sfmn(ii) = sfmn(ii) + Azo(lvol,aderiv,ii)%s(ll) * RTT(ll,mi,iocons,1) * half endif - enddo ! end of do ll; 20 Feb 13; - else - do ll = 0, Lrad(lvol) ! loop over Zernike polynomials; Lrad is the radial resolution; 01 Jul 19; + else ; ; efmn(ii) = efmn(ii) + Ate(lvol,aderiv,ii)%s(ll) * RTT(ll,mi,iocons,0) ; ; cfmn(ii) = cfmn(ii) + Aze(lvol,aderiv,ii)%s(ll) * RTT(ll,mi,iocons,0) if( NOTstellsym ) then ; ofmn(ii) = ofmn(ii) + Ato(lvol,aderiv,ii)%s(ll) * RTT(ll,mi,iocons,0) ; ; sfmn(ii) = sfmn(ii) + Azo(lvol,aderiv,ii)%s(ll) * RTT(ll,mi,iocons,0) endif - enddo - endif + endif + enddo ! end of do ll; 20 Feb 13; else do ll = 0, Lrad(lvol) ! loop over Chebyshev polynomials; Lrad is the radial resolution; ; ; efmn(ii) = efmn(ii) + Ate(lvol,aderiv,ii)%s(ll) * TT(ll,iocons,tderiv) ! aderiv labels deriv. wrt mu, pflux; @@ -1428,336 +1426,6 @@ subroutine broadcast_inputs LlBCAST( Wreadin, 1, 0 ) LlBCAST( Wwrtend, 1, 0 ) LlBCAST( Wmacros, 1, 0 ) -!TODO<<<<<<< HEAD -!======= -! -!!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! -! -!! set internal parameters that depend on physicslist; -! -!!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! -! -! select case( Istellsym ) -! case( 0 ) ; YESstellsym = .false. ; NOTstellsym = .true. -! case( 1 ) ; YESstellsym = .true. ; NOTstellsym = .false. -! case default ; -! FATAL( readin, .true., illegal Istellsym ) -! end select -! -!!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! -! -!!latex \subsubsection{\type{Mvol} : total number of volumes} -!!latex \begin{enumerate} -!!latex \item The number of plasma volumes is \internal{Mvol}=\inputvar{Nvol}+\inputvar{Lfreebound}; -!!latex \end{enumerate} -! -! FATAL( readin, Lfreebound.lt.0 .or. Lfreebound.gt.1, illegal Lfreebound ) -! -! Mvol = Nvol + Lfreebound -! -!!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! -! -! SALLOCATE( beltramierror,(1:Mvol,1:3), zero) -! -!!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! -! -!!latex \subsubsection{\type{mn}, \type{im(1:mn)} and \type{in(1:mn)} : Fourier mode identification} -!!latex \begin{enumerate} -!!latex \item The Fourier description of even periodic functions is -!!latex \be f(\t,\z) = \sum_{n=0}^{N} f_{0,n} \cos(-n\z) + \sum_{m=1}^{M}\sum_{n=-N}^{N} f_{m,n} \cos(m\t-n\z), -!!latex \ee -!!latex where the resolution is given on input, $M\equiv $ \inputvar{ Mpol} and $N\equiv $ \inputvar{ Ntor}. -!!latex \item For convenience, the Fourier summations are written as -!!latex \be f(\s,\t,\z) &=& \sum_j f_j(s) \cos( m_j \t - n_j \z ), -!!latex \ee -!!latex for $j=1,$ \type{mn}, where \type{mn}$ = N + 1 + M ( 2 N + 1 )$. -!!latex \item The integer arrays \type{im(1:mn)} and \type{in(1:mn)} contain the $m_j$ and $n_j$. -!!latex \item The array \type{in} includes the \type{Nfp} factor. -!!latex \end{enumerate} -! -! mn = 1 + Ntor + Mpol * ( 2 * Ntor + 1 ) ! Fourier resolution of interface geometry & vector potential; -! -! SALLOCATE( im, (1:mn), 0 ) -! SALLOCATE( in, (1:mn), 0 ) -! -! call gi00ab( Mpol, Ntor, Nfp, mn, im(1:mn), in(1:mn) ) ! this sets the im and in mode identification arrays; -! -!!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! -! -!!latex \subsubsection{\type{halfmm(1:mn}, regumm(1:mn) : regularization factor} -!!latex \begin{enumerate} -!!latex \item The ``regularization'' factor, \type{halfmm(1:mn)} = \type{im(1:mn)} * \type{half}, is real. -!!latex \item This is used in \link{lforce}, \link{bfield}, \link{stzxyz}, \link{coords}, \link{jo00aa}, \link{ma00aa}, \link{sc00aa} and \link{tr00ab}. -!!latex \end{enumerate} -! -! SALLOCATE( halfmm, (1:mn), im(1:mn) * half ) -! SALLOCATE( regumm, (1:mn), im(1:mn) * half ) -! -! if( Mregular.ge.2 ) then -! -! where( im.gt.Mregular ) regumm = Mregular * half -! -! endif -! -!! if( myid.eq.0 ) write(ounit,'("global : " 10x " : "i3") im ="i3" , halfmm ="f5.1" , regum ="f5.1" ;")') ( ii, im(ii), halfmm(ii), regumm(ii), ii = 1, mn ) -! -!!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! -! -!!latex \subsubsection{\type{ime} and \type{ine} : extended resolution Fourier mode identification} -!!latex \begin{enumerate} -!!latex \item The ``extended'' Fourier resolution is defined by \internal{lMpol} $ = 4 $ \inputvar{Mpol}, \internal{lNtor} $ = 4 $\inputvar{Ntor}. -!!latex \end{enumerate} -! -!! lMpol = Mpol ; lNtor = Ntor ! no enhanced resolution for metrics; -!! lMpol = 2*Mpol ; lNtor = 2*Ntor ! enhanced resolution for metrics; -! lMpol = 4*Mpol ; lNtor = 4*Ntor ! extra-enhanced resolution for metrics; -! -! mne = 1 + lNtor + lMpol * ( 2 * lNtor + 1 ) ! resolution of metrics; enhanced resolution; see metrix; -! -! SALLOCATE( ime, (1:mne), 0 ) -! SALLOCATE( ine, (1:mne), 0 ) -! -! call gi00ab( lMpol, lNtor, Nfp, mne, ime(1:mne), ine(1:mne) ) -! -!!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! -! -!!latex \subsubsection{\type{mns}, \type{ims} and \type{ins} : Fourier mode identification for straight-fieldline angle} -! -! lMpol = iMpol ; lNtor = iNtor -! -! if( iMpol.le.0 ) lMpol = Mpol - iMpol -! if( iNtor.le.0 ) lNtor = Ntor - iNtor -! if( Ntor.eq.0 ) lNtor = 0 -! -! mns = 1 + lNtor + lMpol * ( 2 * lNtor + 1 ) ! resolution of straight-field line transformation on interfaces; see tr00ab; soon to be redundant; -! -! SALLOCATE( ims, (1:mns), 0 ) -! SALLOCATE( ins, (1:mns), 0 ) -! -! call gi00ab( lMpol, lNtor, Nfp, mns, ims(1:mns), ins(1:mns) ) ! note that the field periodicity factor is included in ins; -! -!!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! -! -!! set internal parameters that depend on numericlist; -! -!!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! -! -!! set internal parameters that depend on locallist; -! -!!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! -! -!! set internal parameters that depend on globallist; -! -!!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! -! -!! set internal parameters that depend on diagnosticslist; -! -! if( Lcheck.eq.5 ) then ; forcetol = 1.0e+12 ; nPpts = 0 ! will check Hessian using finite-differences; -! endif -! -!!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! -! -!!latex \subsubsection{\type{iRbc(1:mn,0:Mvol}, \type{iZbs(1:mn,0:Mvol}, \type{iRbs(1:mn,0:Mvol} and \type{iZbc(1:mn,0:Mvol} : geometry} -! -!!latex \begin{enumerate} -!!latex \item \type{iRbc}, \type{iZbs}, \type{iRbs} and \type{iZbc} : Fourier harmonics of interface geometry; -!!latex \item \type{iVns}, \type{iVnc}, \type{iBns} and \type{iBns} : Fourier harmonics of normal field at computational boundary; -!!latex \end{enumerate} -! -! SALLOCATE( iRbc, (1:mn,0:Mvol), zero ) ! interface Fourier harmonics; -! SALLOCATE( iZbs, (1:mn,0:Mvol), zero ) -! SALLOCATE( iRbs, (1:mn,0:Mvol), zero ) -! SALLOCATE( iZbc, (1:mn,0:Mvol), zero ) -! -! if( Lperturbed.eq.1 ) then -! SALLOCATE( dRbc, (1:mn,0:Mvol), zero ) ! interface Fourier harmonics; -! SALLOCATE( dZbs, (1:mn,0:Mvol), zero ) -! SALLOCATE( dRbs, (1:mn,0:Mvol), zero ) -! SALLOCATE( dZbc, (1:mn,0:Mvol), zero ) -! endif -! -!!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! -! -! SALLOCATE( iVns, (1:mn), zero ) -! SALLOCATE( iBns, (1:mn), zero ) -! SALLOCATE( iVnc, (1:mn), zero ) -! SALLOCATE( iBnc, (1:mn), zero ) -! -! !SALLOCATE( lRbc, (1:mn), zero ) ! not used; SRH: 27 Feb 18; -! !SALLOCATE( lZbs, (1:mn), zero ) -! !SALLOCATE( lRbs, (1:mn), zero ) -! !SALLOCATE( lZbc, (1:mn), zero ) -! -!!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! -! -!!latex \subsubsection{\type{ajk} : construction of coordinate axis} -! -!!latex \begin{enumerate} -!!latex \item This is only used in \link{rzaxis} to perform the poloidal integration and is defined quite simply: \newline -!!latex \internal{ajk[i]} $\equiv 2\pi$ if $m_i = 0$, and \newline -!!latex \internal{ajk[i]} $\equiv 0 $ if $m_i \ne 0$. -!!latex \end{enumerate} -! -! SALLOCATE( ajk, (1:mn), zero ) ! this must be allocated & assigned now, as it is used in readin; primarily used in packxi; 02 Jan 15; -! -! do kk = 1, mn ; mk = im(kk) ; nk = in(kk) -! -! if( mk.eq.0 ) ajk(kk) = pi2 -! -! enddo ! end of do kk; -! -!!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! -! -! if( myid.eq.0 ) then ! read plasma boundary & computational boundary; initialize interface geometry; -! -! if( Igeometry.eq.3 .and. Rbc(0,+1)+Rbc(0,-1).gt.zero .and. Zbs(0,+1)-Zbs(0,-1).gt.zero ) then ; Lchangeangle = .true. -! else ; Lchangeangle = .false. -! endif -! -! if( Lchangeangle ) write(ounit,'("readin : " 10x " : CHANGING ANGLE ;")') -! -! do ii = 1, mn ; mm = im(ii) ; nn = in(ii) / Nfp ! set plasma boundary, computational boundary; 29 Apr 15; -! -! if( Lchangeangle ) then ; jj = -1 ; kk = -nn ! change sign of poloidal angle; -! else ; jj = +1 ; kk = +nn -! endif -! -! if( mm.eq.0 .and. nn.eq.0 ) then -! -! ;iRbc(ii,Nvol) = Rbc( nn, mm) ! plasma boundary is ALWAYS given by namelist Rbc & Zbs; -! ;iZbs(ii,Nvol) = zero -! if( NOTstellsym ) then -! ;iRbs(ii,Nvol) = zero -! ;iZbc(ii,Nvol) = Zbc( nn, mm) -! else -! ;iRbs(ii,Nvol) = zero -! ;iZbc(ii,Nvol) = zero -! endif -! -! if( Lfreebound.eq.1 ) then -! -! iRbc(ii,Mvol) = Rwc( nn, mm) ! computational boundary is ALWAYS given by namelist Rbc & Zbs; -! iZbs(ii,Mvol) = zero -! if( NOTstellsym ) then -! iRbs(ii,Mvol) = zero -! iZbc(ii,Mvol) = Zwc( nn, mm) -! else -! iRbs(ii,Mvol) = zero -! iZbc(ii,Mvol) = zero -! endif -! -! iVns(ii ) = zero -! iBns(ii ) = zero -! if( NOTstellsym ) then -! iVnc(ii ) = Vnc( nn, mm) ! I guess that this must be zero, because \div B = 0 ; -! iBnc(ii ) = Bnc( nn, mm) ! I guess that this must be zero, because \div B = 0 ; -! else -! iVnc(ii ) = zero -! iBnc(ii ) = zero -! endif -! -! endif ! end of if( Lfreebound.eq.1 ) ; -! -! else ! if( mm.eq.0 .and. nn.eq.0 ) then ; matches -! -! ;iRbc(ii,Nvol) = Rbc( kk, mm) + Rbc(-kk,-mm) ! plasma boundary is ALWAYS given by namelist Rbc & Zbs; -! ;iZbs(ii,Nvol) = ( Zbs( kk, mm) - Zbs(-kk,-mm) ) * jj -! if( NOTstellsym ) then -! ;iRbs(ii,Nvol) = ( Rbs( kk, mm) - Rbs(-kk,-mm) ) * jj -! ;iZbc(ii,Nvol) = Zbc( kk, mm) + Zbc(-kk,-mm) -! else -! ;iRbs(ii,Nvol) = zero -! ;iZbc(ii,Nvol) = zero -! endif -! -! if( Lfreebound.eq.1 ) then -! -! iRbc(ii,Mvol) = Rwc( kk, mm) + Rwc(-kk,-mm) ! computational boundary is ALWAYS given by namelist Rbc & Zbs; -! iZbs(ii,Mvol) = ( Zws( kk, mm) - Zws(-kk,-mm) ) * jj -! if( NOTstellsym ) then -! iRbs(ii,Mvol) = ( Rws( kk, mm) - Rws(-kk,-mm) ) * jj -! iZbc(ii,Mvol) = Zwc( kk, mm) + Zwc(-kk,-mm) -! else -! iRbs(ii,Mvol) = zero -! iZbc(ii,Mvol) = zero -! endif -! -! endif ! matches if( Lfreebound.eq.1 ) ; -! -! iVns(ii ) = ( Vns( kk, mm) - Vns(-kk,-mm) ) * jj -! iBns(ii ) = ( Bns( kk, mm) - Bns(-kk,-mm) ) * jj -! if( NOTstellsym ) then -! iVnc(ii ) = Vnc( kk, mm) + Vnc(-kk,-mm) -! iBnc(ii ) = Bnc( kk, mm) + Bnc(-kk,-mm) -! else -! iVnc(ii ) = zero -! iBnc(ii ) = zero -! endif -! -! endif ! end of if( mm.eq.0 .and. nn.eq.0 ) ; -! -! enddo ! end of do ii = 1, mn; -! -! -! select case( Linitialize ) ! 24 Oct 12; -! -! case( :0 ) ! Linitialize=0 ; initial guess for geometry of the interior surfaces is given in the input file; -! -! SALLOCATE( RZRZ, (1:4,1:Nvol), zero ) ! temp array for reading input; -! -! if( Lchangeangle ) then ; jj = -1 ! change sign of poloidal angle; Loizu Nov 18; -! else ; jj = +1 -! endif -! -! do ! will read in Fourier harmonics until the end of file is reached; -! -! read(iunit,*,iostat=ios) mm, nn, RZRZ(1:4,1:Nvol) !if change of angle applies, transformation assumes m>=0 and for m=0 only n>=0; -! if( ios.ne.0 ) exit -! -! do ii = 1, mn ; mi = im(ii) ; ni = in(ii) ! loop over harmonics within range; -! if( mm.eq.0 .and. mi.eq.0 .and. nn*Nfp.eq.ni ) then -! iRbc(ii,1:Nvol-1) = RZRZ(1,1:Nvol-1) ! select relevant harmonics; -! iZbs(ii,1:Nvol-1) = RZRZ(2,1:Nvol-1) ! select relevant harmonics; -! if( NOTstellsym ) then -! iRbs(ii,1:Nvol-1) = RZRZ(3,1:Nvol-1) ! select relevant harmonics; -! iZbc(ii,1:Nvol-1) = RZRZ(4,1:Nvol-1) ! select relevant harmonics; -! else -! iRbs(ii,1:Nvol-1) = zero ! select relevant harmonics; -! iZbc(ii,1:Nvol-1) = zero ! select relevant harmonics; -! endif -! elseif( mm.eq.mi .and. nn*Nfp.eq.jj*ni ) then -! iRbc(ii,1:Nvol-1) = RZRZ(1,1:Nvol-1) ! select relevant harmonics; -! iZbs(ii,1:Nvol-1) = jj*RZRZ(2,1:Nvol-1) ! select relevant harmonics; -! if( NOTstellsym ) then -! iRbs(ii,1:Nvol-1) = jj*RZRZ(3,1:Nvol-1) ! select relevant harmonics; -! iZbc(ii,1:Nvol-1) = RZRZ(4,1:Nvol-1) ! select relevant harmonics; -! else -! iRbs(ii,1:Nvol-1) = zero ! select relevant harmonics; -! iZbc(ii,1:Nvol-1) = zero ! select relevant harmonics; -! endif -! endif -! enddo ! end of do ii; -! -! enddo ! end of do; -! -! DALLOCATE(RZRZ) -! -! end select ! end select case( Linitialize ); -! -! if( Igeometry.eq.3 ) then -! if( Rac(0).gt.zero ) then ! user has supplied logically possible coordinate axis; -! iRbc(1:Ntor+1,0) = Rac(0:Ntor) -! iZbs(1:Ntor+1,0) = Zas(0:Ntor) -! iRbs(1:Ntor+1,0) = Ras(0:Ntor) -! iZbc(1:Ntor+1,0) = Zac(0:Ntor) -! else ! see preset for poloidal-average specification of coordinate axis and geometrical initialization; -! endif ! end of if( Igeometry.eq.3 ) then ; -! endif -! -! endif ! end of if myid.eq.0 loop; only the master will read the input file; all variables need to be broadcast; -! -!!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! -!>>>>>>> origin/Lconstraint_-2 end subroutine ! broadcast_inputs diff --git a/preset.f90 b/preset.f90 index 3da41b4d..41dfe9d5 100644 --- a/preset.f90 +++ b/preset.f90 @@ -53,6 +53,7 @@ subroutine preset ! set internal parameters that depend on physicslist; + select case( Istellsym ) case( 0 ) ; YESstellsym = .false. ; NOTstellsym = .true. case( 1 ) ; YESstellsym = .true. ; NOTstellsym = .false. @@ -298,6 +299,8 @@ subroutine preset iZbc(ii,Mvol) = zero endif + endif ! matches if( Lfreebound.eq.1 ) ; + iVns(ii ) = ( Vns( kk, mm) - Vns(-kk,-mm) ) * jj iBns(ii ) = ( Bns( kk, mm) - Bns(-kk,-mm) ) * jj if( NOTstellsym ) then @@ -308,8 +311,6 @@ subroutine preset iBnc(ii ) = zero endif - endif ! matches if( Lfreebound.eq.1 ) ; - endif ! end of if( mm.eq.0 .and. nn.eq.0 ) ; enddo ! end of do ii = 1, mn; From abb700b1481aba4a376c2efbba6d6bb00986371f Mon Sep 17 00:00:00 2001 From: Richard Nies Date: Thu, 27 May 2021 15:45:44 -0400 Subject: [PATCH 04/20] NUMBA accelerated Zernike construction. --- Utilities/pythontools/py_spec/output/_processing.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/Utilities/pythontools/py_spec/output/_processing.py b/Utilities/pythontools/py_spec/output/_processing.py index 72c9998f..b8172cee 100644 --- a/Utilities/pythontools/py_spec/output/_processing.py +++ b/Utilities/pythontools/py_spec/output/_processing.py @@ -1,4 +1,5 @@ import numpy as np +from numba import jit def get_grid_and_jacobian_and_metric( @@ -322,6 +323,7 @@ def get_B_covariant(self, Bcontrav, g): return Bco +@jit(nopython=True) def _get_zernike(sarr, lrad, mpol): """ Get the value of the zernike polynomials, their first and second derivatives @@ -335,7 +337,7 @@ def _get_zernike(sarr, lrad, mpol): rm1 = np.zeros_like(r) # r to the power of m-1'th rm2 = np.zeros_like(r) # r to the power of m-2'th - zernike = np.zeros([ns, lrad + 1, mpol + 1], dtype=np.float64) + zernike = np.zeros(shape=(ns, lrad + 1, mpol + 1), dtype=np.float64) dzernike = np.zeros_like(zernike) d2zernike = np.zeros_like(zernike) From f01ecaad2c7e08c010fd38e55e457de3f23d1a25 Mon Sep 17 00:00:00 2001 From: Richard Nies Date: Sat, 17 Jul 2021 16:14:31 -0400 Subject: [PATCH 05/20] More effective reading in of magnetic field. --- .../pythontools/py_spec/output/_processing.py | 188 +++++++++++------- Utilities/pythontools/py_spec/output/spec.py | 4 +- 2 files changed, 122 insertions(+), 70 deletions(-) diff --git a/Utilities/pythontools/py_spec/output/_processing.py b/Utilities/pythontools/py_spec/output/_processing.py index b8172cee..c7c4b1fa 100644 --- a/Utilities/pythontools/py_spec/output/_processing.py +++ b/Utilities/pythontools/py_spec/output/_processing.py @@ -174,8 +174,8 @@ def get_B( tarr=np.linspace(0, 0, 1), zarr=np.linspace(0, 0, 1), ): - Bcontra_and_s_der = get_B_and_s_der(self, lvol, jacobian, sarr, tarr, zarr) - return Bcontra_and_s_der[0:3] + Bcontra_and_der = get_B_and_der(self, lvol, jacobian, sarr, tarr, zarr) + return Bcontra_and_der[0:3] def get_s_der_B( self, @@ -185,10 +185,33 @@ def get_s_der_B( tarr=np.linspace(0, 0, 1), zarr=np.linspace(0, 0, 1), ): - Bcontra_and_s_der = get_B_and_s_der(self, lvol, jacobian, sarr, tarr, zarr) - return Bcontra_and_s_der[3:6] + Bcontra_and_der = get_B_and_der(self, lvol, jacobian, sarr, tarr, zarr) + return Bcontra_and_der[3:6] -def get_B_and_s_der( +def get_t_der_B( + self, + lvol=0, + jacobian=None, + sarr=np.linspace(1, 1, 1), + tarr=np.linspace(0, 0, 1), + zarr=np.linspace(0, 0, 1), +): + Bcontra_and_der = get_B_and_der(self, lvol, jacobian, sarr, tarr, zarr) + return Bcontra_and_der[6:9] + +def get_z_der_B( + self, + lvol=0, + jacobian=None, + sarr=np.linspace(1, 1, 1), + tarr=np.linspace(0, 0, 1), + zarr=np.linspace(0, 0, 1), +): + Bcontra_and_der = get_B_and_der(self, lvol, jacobian, sarr, tarr, zarr) + return Bcontra_and_der[9:12] + + +def get_B_and_der( self, lvol=0, jacobian=None, @@ -226,87 +249,114 @@ def get_B_and_s_der( LZernike = self.input.physics.Igeometry > 1 and lvol == 0 + if not LZernike: + # make basis recombination for cheb + lcoeff = np.arange(0, Lrad + 1) + 1 + + Ate = Ate / lcoeff[None, :] + Aze = Aze / lcoeff[None, :] + Ato = Ato / lcoeff[None, :] + Azo = Azo / lcoeff[None, :] + + Ate[:, 0] = np.sum(Ate * (-1.0) ** lcoeff[None, :], 1) + Aze[:, 0] = np.sum(Aze * (-1.0) ** lcoeff[None, :], 1) + Ato[:, 0] = np.sum(Ato * (-1.0) ** lcoeff[None, :], 1) + Azo[:, 0] = np.sum(Azo * (-1.0) ** lcoeff[None, :], 1) + + + c = ( + im[nax, :, nax, nax] * Azo.T[:, :, nax, nax] + + in_[nax, :, nax, nax] * Ato.T[:, :, nax, nax] + ) * cosa[nax, :, :, :] - ( + im[nax, :, nax, nax] * Aze.T[:, :, nax, nax] + + in_[nax, :, nax, nax] * Ate.T[:, :, nax, nax] + ) * sina[nax, :, :, :] + + ct = im[nax, :, nax, nax] * ( - ( + im[nax, :, nax, nax] * Azo.T[:, :, nax, nax] + + in_[nax, :, nax, nax] * Ato.T[:, :, nax, nax] + ) * sina[nax, :, :, :] - ( + im[nax, :, nax, nax] * Aze.T[:, :, nax, nax] + + in_[nax, :, nax, nax] * Ate.T[:, :, nax, nax] + ) * cosa[nax, :, :, :] ) + + cz = in_[nax, :, nax, nax] * ( ( + im[nax, :, nax, nax] * Azo.T[:, :, nax, nax] + + in_[nax, :, nax, nax] * Ato.T[:, :, nax, nax] + ) * sina[nax, :, :, :] + ( + im[nax, :, nax, nax] * Aze.T[:, :, nax, nax] + + in_[nax, :, nax, nax] * Ate.T[:, :, nax, nax] + ) * cosa[nax, :, :, :] ) + + + c1 = ( + Aze.T[:, :, nax, nax] * cosa[nax, :, :, :] + + Azo.T[:, :, nax, nax] * sina[nax, :, :, :]) + + c1t = im[nax, :, nax, nax] * ( + - Aze.T[:, :, nax, nax] * sina[nax, :, :, :] + + Azo.T[:, :, nax, nax] * cosa[nax, :, :, :]) + + c1z = in_[nax, :, nax, nax] * ( + + Aze.T[:, :, nax, nax] * sina[nax, :, :, :] + - Azo.T[:, :, nax, nax] * cosa[nax, :, :, :]) + + + c2 = ( + Ate.T[:, :, nax, nax] * cosa[nax, :, :, :] + + Ato.T[:, :, nax, nax] * sina[nax, :, :, :]) + + c2t = im[nax, :, nax, nax] * ( + - Ate.T[:, :, nax, nax] * sina[nax, :, :, :] + + Ato.T[:, :, nax, nax] * cosa[nax, :, :, :]) + + c2z = in_[nax, :, nax, nax] * ( + + Ate.T[:, :, nax, nax] * sina[nax, :, :, :] + - Ato.T[:, :, nax, nax] * cosa[nax, :, :, :]) + + if LZernike: # Zernike polynomial being used from ._processing import _get_zernike zernike, dzernike, d2zernike = _get_zernike(sbar, Lrad, Mpol) - c = ( - im[nax, :, nax, nax] * Azo.T[:, :, nax, nax] - + in_[nax, :, nax, nax] * Ato.T[:, :, nax, nax] - ) * cosa[nax, :, :, :] - ( - im[nax, :, nax, nax] * Aze.T[:, :, nax, nax] - + in_[nax, :, nax, nax] * Ate.T[:, :, nax, nax] - ) * sina[ - nax, :, :, : - ] - - Bs = np.sum( zernike[:, :, im, None, None] * c[None, :, :, :, :], axis=(1, 2)) - dsBs = np.sum(dzernike[:, :, im, None, None] * c[None, :, :, :, :], axis=(1, 2)) - - c1 = ( - Aze.T[:, :, nax, nax] * cosa[nax, :, :, :] - + Azo.T[:, :, nax, nax] * sina[nax, :, :, :] - ) - Bt = -np.sum( dzernike[:, :, im, None, None] * c1[None, :, :, :, :], axis=(1, 2)) - dsBt = -np.sum(d2zernike[:, :, im, None, None] * c1[None, :, :, :, :], axis=(1, 2)) + Bs = np.sum( zernike[:, :, im, None, None] * c[None, :, :, :, :], axis=(1, 2)) + dsBs = np.sum(dzernike[:, :, im, None, None] * c[None, :, :, :, :], axis=(1, 2)) + dtBs = np.sum( zernike[:, :, im, None, None] * ct[None, :, :, :, :], axis=(1, 2)) + dzBs = np.sum( zernike[:, :, im, None, None] * cz[None, :, :, :, :], axis=(1, 2)) - c2 = ( - Ate.T[:, :, nax, nax] * cosa[nax, :, :, :] - + Ato.T[:, :, nax, nax] * sina[nax, :, :, :] - ) + Bt = -np.sum( dzernike[:, :, im, None, None] * c1[None, :, :, :, :], axis=(1, 2)) + dsBt = -np.sum(d2zernike[:, :, im, None, None] * c1[None, :, :, :, :], axis=(1, 2)) + dtBt = -np.sum( dzernike[:, :, im, None, None] * c1t[None, :, :, :, :], axis=(1, 2)) + dzBt = -np.sum( dzernike[:, :, im, None, None] * c1z[None, :, :, :, :], axis=(1, 2)) - Bz = np.sum( dzernike[:, :, im, None, None] * c2[None, :, :, :, :], axis=(1, 2)) - dsBz = np.sum(d2zernike[:, :, im, None, None] * c2[None, :, :, :, :], axis=(1, 2)) + Bz = np.sum( dzernike[:, :, im, None, None] * c2[None, :, :, :, :], axis=(1, 2)) + dsBz = np.sum(d2zernike[:, :, im, None, None] * c2[None, :, :, :, :], axis=(1, 2)) + dtBz = np.sum( dzernike[:, :, im, None, None] * c2t[None, :, :, :, :], axis=(1, 2)) + dzBz = np.sum( dzernike[:, :, im, None, None] * c2z[None, :, :, :, :], axis=(1, 2)) else: # Chebyshev being used import numpy.polynomial.chebyshev as Cheb - lcoeff = np.arange(0, Lrad + 1) + 1 - # make basis recombination for cheb - Ate = Ate / lcoeff[None, :] - Aze = Aze / lcoeff[None, :] - Ato = Ato / lcoeff[None, :] - Azo = Azo / lcoeff[None, :] - - Ate[:, 0] = np.sum(Ate * (-1.0) ** lcoeff[None, :], 1) - Aze[:, 0] = np.sum(Aze * (-1.0) ** lcoeff[None, :], 1) - Ato[:, 0] = np.sum(Ato * (-1.0) ** lcoeff[None, :], 1) - Azo[:, 0] = np.sum(Azo * (-1.0) ** lcoeff[None, :], 1) - - # Ch ,mn,t ,z - c = ( - im[nax, :, nax, nax] * Azo.T[:, :, nax, nax] - + in_[nax, :, nax, nax] * Ato.T[:, :, nax, nax] - ) * cosa[nax, :, :, :] - ( - im[nax, :, nax, nax] * Aze.T[:, :, nax, nax] - + in_[nax, :, nax, nax] * Ate.T[:, :, nax, nax] - ) * sina[ - nax, :, :, : - ] - Bs = np.rollaxis(np.sum(Cheb.chebval(sarr, c), axis=0), 2) dsBs = np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c)), axis=0), 2) + dtBs = np.rollaxis(np.sum(Cheb.chebval(sarr, ct), axis=0), 2) + dzBs = np.rollaxis(np.sum(Cheb.chebval(sarr, cz), axis=0), 2) - c1 = ( - Aze.T[:, :, nax, nax] * cosa[nax, :, :, :] - + Azo.T[:, :, nax, nax] * sina[nax, :, :, :] - ) - Bt = -np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c1)), axis=0),2) - dsBt = -np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c1,m=2)), axis=0),2) + Bt = -np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c1)), axis=0),2) + dsBt = -np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c1,m=2)), axis=0),2) + dtBt = -np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c1t)), axis=0),2) + dzBt = -np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c1z)), axis=0),2) - c2 = ( - Ate.T[:, :, nax, nax] * cosa[nax, :, :, :] - + Ato.T[:, :, nax, nax] * sina[nax, :, :, :] - ) - Bz = np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c2)), axis=0),2) - dsBz = np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c2,m=2)), axis=0),2) - + Bz = np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c2)), axis=0),2) + dsBz = np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c2,m=2)), axis=0),2) + dtBz = np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c2t)), axis=0),2) + dzBz = np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c2z)), axis=0),2) - Bcontrav_and_s_der = np.array([Bs, Bt, Bz, dsBs, dsBt, dsBz]) / jacobian - return Bcontrav_and_s_der + Bcontrav_and_der = np.array([Bs, Bt, Bz, dsBs, dsBt, dsBz, dtBs, dtBt, dtBz, dzBs, dzBt, dzBz]) / jacobian + return Bcontrav_and_der # Bcontrav = get_B(s,lvol=lvol,jacobian=jacobian,sarr=sarr,tarr=tarr,zarr=zarr) @@ -323,7 +373,7 @@ def get_B_covariant(self, Bcontrav, g): return Bco -@jit(nopython=True) +#@jit(nopython=True) def _get_zernike(sarr, lrad, mpol): """ Get the value of the zernike polynomials, their first and second derivatives diff --git a/Utilities/pythontools/py_spec/output/spec.py b/Utilities/pythontools/py_spec/output/spec.py index 3a2cb3ac..fd1ffa92 100644 --- a/Utilities/pythontools/py_spec/output/spec.py +++ b/Utilities/pythontools/py_spec/output/spec.py @@ -36,7 +36,9 @@ class SPECout: metric, get_B, get_s_der_B, - get_B_and_s_der, + get_t_der_B, + get_z_der_B, + get_B_and_der, get_modB, get_B_covariant ) From 4dea0f44c9dfffde37e3589e368ceeaf6849eecc Mon Sep 17 00:00:00 2001 From: Richard Nies Date: Thu, 9 Sep 2021 14:49:35 -0400 Subject: [PATCH 06/20] Clean-up of py_spec processing. --- Utilities/pythontools/py_spec/output/_processing.py | 13 ++++--------- 1 file changed, 4 insertions(+), 9 deletions(-) diff --git a/Utilities/pythontools/py_spec/output/_processing.py b/Utilities/pythontools/py_spec/output/_processing.py index c7c4b1fa..a30159b1 100644 --- a/Utilities/pythontools/py_spec/output/_processing.py +++ b/Utilities/pythontools/py_spec/output/_processing.py @@ -1,6 +1,4 @@ import numpy as np -from numba import jit - def get_grid_and_jacobian_and_metric( self, @@ -227,6 +225,7 @@ def get_B_and_der( # Lrad = s.input.physics.Lrad[lvol] + Ate = self.vector_potential.Ate[lvol] Aze = self.vector_potential.Aze[lvol] Ato = self.vector_potential.Ato[lvol] @@ -263,7 +262,6 @@ def get_B_and_der( Ato[:, 0] = np.sum(Ato * (-1.0) ** lcoeff[None, :], 1) Azo[:, 0] = np.sum(Azo * (-1.0) ** lcoeff[None, :], 1) - c = ( im[nax, :, nax, nax] * Azo.T[:, :, nax, nax] + in_[nax, :, nax, nax] * Ato.T[:, :, nax, nax] @@ -271,7 +269,7 @@ def get_B_and_der( im[nax, :, nax, nax] * Aze.T[:, :, nax, nax] + in_[nax, :, nax, nax] * Ate.T[:, :, nax, nax] ) * sina[nax, :, :, :] - + ct = im[nax, :, nax, nax] * ( - ( im[nax, :, nax, nax] * Azo.T[:, :, nax, nax] + in_[nax, :, nax, nax] * Ato.T[:, :, nax, nax] @@ -313,7 +311,7 @@ def get_B_and_der( c2z = in_[nax, :, nax, nax] * ( + Ate.T[:, :, nax, nax] * sina[nax, :, :, :] - Ato.T[:, :, nax, nax] * cosa[nax, :, :, :]) - + if LZernike: # Zernike polynomial being used @@ -355,13 +353,11 @@ def get_B_and_der( dtBz = np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c2t)), axis=0),2) dzBz = np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c2z)), axis=0),2) + Bcontrav_and_der = np.array([Bs, Bt, Bz, dsBs, dsBt, dsBz, dtBs, dtBt, dtBz, dzBs, dzBt, dzBz]) / jacobian return Bcontrav_and_der -# Bcontrav = get_B(s,lvol=lvol,jacobian=jacobian,sarr=sarr,tarr=tarr,zarr=zarr) - - def get_modB(self, Bcontrav, g): """Input - Bcontrav has to come from get_B function""" modB = np.sqrt(np.einsum("iabc,jiabc,jabc->abc", Bcontrav, g, Bcontrav)) @@ -373,7 +369,6 @@ def get_B_covariant(self, Bcontrav, g): return Bco -#@jit(nopython=True) def _get_zernike(sarr, lrad, mpol): """ Get the value of the zernike polynomials, their first and second derivatives From aefcc777a6b426d7ca7dee68094aac6dc894a543 Mon Sep 17 00:00:00 2001 From: Richard Nies Date: Thu, 9 Sep 2021 15:27:33 -0400 Subject: [PATCH 07/20] minor clean-ups --- brcast.f90 | 2 +- ma02aa.f90 | 6 +++--- preset.f90 | 1 - 3 files changed, 4 insertions(+), 5 deletions(-) diff --git a/brcast.f90 b/brcast.f90 index 0efd2ca3..9dfb9de6 100644 --- a/brcast.f90 +++ b/brcast.f90 @@ -134,7 +134,7 @@ subroutine brcast( lvol ) !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! ! if( lvol.gt.Nvol .and. Lconstraint.eq.-1 .and. Wcurent ) then ! 27 Feb 17; - if( (lvol.gt.Nvol .or. Lconstraint .eq.-2) .and. Wcurent ) then ! 27 Feb 17; + if( (lvol.gt.Nvol .or. Lconstraint .eq.-2) .and. Wcurent ) then ! 27 Feb 17; !write(ounit,'("brcast : " 10x " : myid="i3" ; broadcasting : curtor="es13.5" ; curpol="es13.5" ;")') myid, curtor, curpol RlBCAST( curtor, 1, llmodnp ) RlBCAST( curpol, 1, llmodnp ) diff --git a/ma02aa.f90 b/ma02aa.f90 index 0619a11a..87b32f18 100644 --- a/ma02aa.f90 +++ b/ma02aa.f90 @@ -418,9 +418,9 @@ subroutine ma02aa( lvol, NN ) if( LBlinear ) then ! assume Beltrami field is parameterized by helicity multiplier (and poloidal flux); - - lastcpu = GETTIME - + + lastcpu = GETTIME + if( Lplasmaregion ) then Xdof(1:2) = xoffset + (/ mu(lvol), dpflux(lvol) /) ! initial guess for degrees of freedom; offset from zero so that relative error is small; diff --git a/preset.f90 b/preset.f90 index a710e37c..d066478f 100644 --- a/preset.f90 +++ b/preset.f90 @@ -53,7 +53,6 @@ subroutine preset ! set internal parameters that depend on physicslist; - select case( Istellsym ) case( 0 ) ; YESstellsym = .false. ; NOTstellsym = .true. case( 1 ) ; YESstellsym = .true. ; NOTstellsym = .false. From ad9f86d569c9330a8084ce9f240e6703e1a3563e Mon Sep 17 00:00:00 2001 From: Richard Nies Date: Thu, 9 Sep 2021 15:29:24 -0400 Subject: [PATCH 08/20] Further minor clean-ups. --- Utilities/pythontools/py_spec/output/_processing.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/Utilities/pythontools/py_spec/output/_processing.py b/Utilities/pythontools/py_spec/output/_processing.py index a30159b1..43fe0ef2 100644 --- a/Utilities/pythontools/py_spec/output/_processing.py +++ b/Utilities/pythontools/py_spec/output/_processing.py @@ -1,5 +1,6 @@ import numpy as np + def get_grid_and_jacobian_and_metric( self, lvol=0, @@ -164,6 +165,7 @@ def metric( _, _, _, g = get_grid_and_jacobian_and_metric(self, lvol, sarr, tarr, zarr) return g + def get_B( self, lvol=0, @@ -175,6 +177,7 @@ def get_B( Bcontra_and_der = get_B_and_der(self, lvol, jacobian, sarr, tarr, zarr) return Bcontra_and_der[0:3] + def get_s_der_B( self, lvol=0, @@ -186,6 +189,7 @@ def get_s_der_B( Bcontra_and_der = get_B_and_der(self, lvol, jacobian, sarr, tarr, zarr) return Bcontra_and_der[3:6] + def get_t_der_B( self, lvol=0, From 9d033643bb206b446956255a9b69b0599370200a Mon Sep 17 00:00:00 2001 From: Richard Nies Date: Mon, 27 Sep 2021 19:50:10 -0400 Subject: [PATCH 09/20] Removed print statements, added checks of inputs. --- global.f90 | 1 + ma02aa.f90 | 3 --- mp00ac.f90 | 10 ---------- 3 files changed, 1 insertion(+), 13 deletions(-) diff --git a/global.f90 b/global.f90 index cc288925..8e564883 100644 --- a/global.f90 +++ b/global.f90 @@ -1076,6 +1076,7 @@ subroutine check_inputs() FATAL( readin, abs(one+gamma).lt.vsmall, 1+gamma appears in denominator in dforce ) !< \todo Please check this; SRH: 27 Feb 18; FATAL( readin, abs(one-gamma).lt.vsmall, 1-gamma appears in denominator in fu00aa ) !< \todo Please check this; SRH: 27 Feb 18; FATAL( readin, Lconstraint.lt.-2 .or. Lconstraint.gt.3, illegal Lconstraint ) + FATAL( readin, Lconstraint.eq.-2 .and. Nvol.ne.1, Lconstraint -2 only for single volume calculation ) FATAL( readin, Igeometry.eq.1 .and. rpol.lt.vsmall, poloidal extent of slab too small or negative ) FATAL( readin, Igeometry.eq.1 .and. rtor.lt.vsmall, toroidal extent of slab too small or negative ) diff --git a/ma02aa.f90 b/ma02aa.f90 index 87b32f18..b6e9c385 100644 --- a/ma02aa.f90 +++ b/ma02aa.f90 @@ -491,7 +491,6 @@ subroutine ma02aa( lvol, NN ) elseif (ihybrj .ge. 0) then dtflux(lvol) = Xdof(1) - xoffset end if - print *,"dtflux(lvol): ", dtflux(lvol) else select case( ihybrj ) case( 0: ) ; mu(lvol) = Xdof(1) - xoffset @@ -553,8 +552,6 @@ subroutine ma02aa( lvol, NN ) endif ! end of if( LBlinear ) then; - print *,"End of ma02aa.f90- dtflux(lvol): ", dtflux(lvol) - !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!- !latex \subsection{debugging: finite-difference confirmation of the derivatives of the rotational-transform} diff --git a/mp00ac.f90 b/mp00ac.f90 index dbf1ea03..b3d6f37f 100644 --- a/mp00ac.f90 +++ b/mp00ac.f90 @@ -223,9 +223,6 @@ subroutine mp00ac( Ndof, Xdof, Fdof, Ddof, Ldfjac, iflag ) ! argument list is fi lmu = mu(lvol) dpf = dpflux(lvol) dtf = Xdof(1) - xoffset - print *,"dtf: ", dtf - print *,"dpf: ", dpf - print *,"lmu: ", lmu end if else ! Lvacuumregion; @@ -617,12 +614,6 @@ subroutine mp00ac( Ndof, Xdof, Fdof, Ddof, Ldfjac, iflag ) ! argument list is fi if( iflag.eq.1 ) Fdof(1 ) = dItGpdxtp(1,0,lvol) - curpol ! Derivative w.r.t. toroidal flux only if( iflag.eq.2 ) Ddof(1,1) = dItGpdxtp(1,1,lvol) - if (iflag .eq. 1) then - print *,"curtor (computed): ", dItGpdxtp(0,0,lvol) - print *,"curpol (computed): ", dItGpdxtp(1,0,lvol) - print *,"curpol (requested): ", curpol - print *,"Fdof: ", Fdof(1) - end if case( -1 ) ! Lconstraint=-1; @@ -789,7 +780,6 @@ subroutine mp00ac( Ndof, Xdof, Fdof, Ddof, Ldfjac, iflag ) ! argument list is fi if ( Lconstraint .eq. -2) then !mu(lvol) = lmu dtflux(lvol) = dtf - print *,"mu(lvol): ", mu(lvol) end if iflag = -2 ! return "acceptance" flag through to ma02aa via ifail; early termination; From eccc08289702ad2e13c4d804002152588e2c968f Mon Sep 17 00:00:00 2001 From: Richard Nies Date: Mon, 27 Sep 2021 20:11:37 -0400 Subject: [PATCH 10/20] Added documentation. --- ma02aa.f90 | 3 ++- mp00ac.f90 | 1 + 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/ma02aa.f90 b/ma02aa.f90 index b6e9c385..e06af56d 100644 --- a/ma02aa.f90 +++ b/ma02aa.f90 @@ -399,6 +399,7 @@ subroutine ma02aa( lvol, NN ) !> to iteratively find the appropriately constrained solution, i.e. \f${\bf f}(\boldsymbol{\mu})=0\f$. !>
  • The function \f${\bf f}(\boldsymbol{\mu})\f$, which is computed by mp00ac(), is defined by the input parameter \c Lconstraint: !>
      +!>
    • If \c Lconstraint = -2, then \f$\boldsymbol{\mu}\f$ is *not* varied and \c Nxdof=0.
      • !>
    • If \c Lconstraint = -1, then \f$\boldsymbol{\mu}\f$ is *not* varied and \c Nxdof=0.
      • !>
    • If \c Lconstraint = 0,2, then \f$\boldsymbol{\mu}\f$ is varied to satisfy the enclosed current constraints, and \c Nxdof=2.
      • !>
    • If \c Lconstraint = 1, then \f$\boldsymbol{\mu}\f$ is varied to satisfy @@ -430,7 +431,7 @@ subroutine ma02aa( lvol, NN ) end if select case( Lconstraint ) - case( -2 ) ; ; Nxdof = 1 ! toroidal flux IS varied to match linking linking current; + case( -2 ) ; ; Nxdof = 1 ! toroidal flux IS varied to match linking current; case( -1 ) ; ; Nxdof = 0 ! multiplier & poloidal flux NOT varied ; ; ; iflag = 1 !(we don't need derivatives) case( 0 ) ; ; Nxdof = 0 ! multiplier & poloidal flux NOT varied diff --git a/mp00ac.f90 b/mp00ac.f90 index b3d6f37f..e8ebae59 100644 --- a/mp00ac.f90 +++ b/mp00ac.f90 @@ -62,6 +62,7 @@ !>
    • For \c Lcoordinatesingularity=T , the returned function is: !> \f{eqnarray}{ {\bf f}(\mu,\Delta\psi_p) \equiv !> \left\{ \begin{array}{cccccccr} +!> (& 0&,& ?&)^T, & \textrm{if }\texttt{Lconstraint} &=& -2 \\ !> (& 0&,& 0&)^T, & \textrm{if }\texttt{Lconstraint} &=& -1 \\ !> (& 0&,& 0&)^T, & \textrm{if }\texttt{Lconstraint} &=& 0 \\ !> (&{{\,\iota\!\!\!}-}(+1)-\texttt{iota (lvol )}&,& 0&)^T, & \textrm{if }\texttt{Lconstraint} &=& 1 \\ From 0030936b7236185168e392205d2dd0e27362807e Mon Sep 17 00:00:00 2001 From: Richard Nies Date: Wed, 29 Sep 2021 11:33:26 -0400 Subject: [PATCH 11/20] Corrected documentation. --- src/ma02aa.f90 | 4 ++-- src/mp00ac.f90 | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/ma02aa.f90 b/src/ma02aa.f90 index e06af56d..bca9fa2c 100644 --- a/src/ma02aa.f90 +++ b/src/ma02aa.f90 @@ -368,7 +368,7 @@ subroutine ma02aa( lvol, NN ) !> !>
        !> -!>
      • In addition to the enclosed toroidal flux, \f$\Delta \psi_t\f$, which is held constant in the plasma volumes, +!>
      • In addition to the enclosed toroidal flux, \f$\Delta \psi_t\f$, which is held constant in the plasma volumes (\c Lconstraint != -2), !> the Beltrami field in a given volume is assumed to be parameterized by \f$\mu\f$ and \f$\Delta \psi_p\f$. !> (Note that \f$\Delta \psi_p\f$ is not defined in a torus.)
        • !>
      • These are "packed" into an array, e.g. \f$\boldsymbol{\mu} \equiv (\mu, \Delta\psi_p)^T\f$, so that standard library routines , @@ -383,6 +383,7 @@ subroutine ma02aa( lvol, NN ) !>
      • \todo If \c Lconstraint = 2, then \f$\mu=\boldsymbol{\mu}_1\f$ is varied in order to satisfy the helicity constraint, !> and \f$\Delta\psi_p=\boldsymbol{\mu}_2\f$ is *not* varied, and \c Nxdof=1. !> (under re-construction) +!>
      • If \c Lconstraint = -2, then \f$\boldsymbol{\mu}\f$ is *not* varied, but \f$\Delta \psi_t\f$ is, and \c Nxdof=1. Only for single-volume calculations.
        • !> !> !>
      • @@ -399,7 +400,6 @@ subroutine ma02aa( lvol, NN ) !> to iteratively find the appropriately constrained solution, i.e. \f${\bf f}(\boldsymbol{\mu})=0\f$. !>
    • The function \f${\bf f}(\boldsymbol{\mu})\f$, which is computed by mp00ac(), is defined by the input parameter \c Lconstraint: !>
        -!>
      • If \c Lconstraint = -2, then \f$\boldsymbol{\mu}\f$ is *not* varied and \c Nxdof=0.
        • !>
      • If \c Lconstraint = -1, then \f$\boldsymbol{\mu}\f$ is *not* varied and \c Nxdof=0.
        • !>
      • If \c Lconstraint = 0,2, then \f$\boldsymbol{\mu}\f$ is varied to satisfy the enclosed current constraints, and \c Nxdof=2.
        • !>
      • If \c Lconstraint = 1, then \f$\boldsymbol{\mu}\f$ is varied to satisfy diff --git a/src/mp00ac.f90 b/src/mp00ac.f90 index e8ebae59..7935b98d 100644 --- a/src/mp00ac.f90 +++ b/src/mp00ac.f90 @@ -62,7 +62,7 @@ !>
      • For \c Lcoordinatesingularity=T , the returned function is: !> \f{eqnarray}{ {\bf f}(\mu,\Delta\psi_p) \equiv !> \left\{ \begin{array}{cccccccr} -!> (& 0&,& ?&)^T, & \textrm{if }\texttt{Lconstraint} &=& -2 \\ +!> (& 0&,& 0&)^T, & \textrm{if }\texttt{Lconstraint} &=& -2 \\ !> (& 0&,& 0&)^T, & \textrm{if }\texttt{Lconstraint} &=& -1 \\ !> (& 0&,& 0&)^T, & \textrm{if }\texttt{Lconstraint} &=& 0 \\ !> (&{{\,\iota\!\!\!}-}(+1)-\texttt{iota (lvol )}&,& 0&)^T, & \textrm{if }\texttt{Lconstraint} &=& 1 \\ @@ -219,7 +219,7 @@ subroutine mp00ac( Ndof, Xdof, Fdof, Ddof, Ldfjac, iflag ) ! argument list is fi else ; dpf = dpflux(lvol) endif - ! only poloidal flux is modified + ! only toroidal flux is modified if (Lconstraint .eq. -2) then lmu = mu(lvol) dpf = dpflux(lvol) From 101ae5b610fc21f0cc393af7e87c3fce3396d182 Mon Sep 17 00:00:00 2001 From: Richard Nies Date: Mon, 11 Oct 2021 16:33:34 +0100 Subject: [PATCH 12/20] Changes to documentation + added Bn!=zero option in fixed-boundary as new option inputlist (Lbdybnzero) --- src/curent.f90 | 29 +++++++++++++++++------------ src/inputlist.f90 | 4 ++++ src/ma02aa.f90 | 4 ++-- src/mp00ac.f90 | 1 + src/preset.f90 | 4 ++++ 5 files changed, 28 insertions(+), 14 deletions(-) diff --git a/src/curent.f90 b/src/curent.f90 index 05f29074..695c08f3 100644 --- a/src/curent.f90 +++ b/src/curent.f90 @@ -63,7 +63,7 @@ subroutine curent( lvol, mn, Nt, Nz, iflag, ldItGp ) use fileunits, only : ounit - use inputlist, only : Wmacros, Wcurent, Lrad, Lconstraint + use inputlist, only : Wmacros, Wcurent, Lrad, Lbdybnzero use cputiming, only : Tcurent @@ -105,7 +105,7 @@ subroutine curent( lvol, mn, Nt, Nz, iflag, ldItGp ) !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! - if (Lconstraint .eq. -2) then + if (.not. Lbdybnzero) then innout = 1.0 lss = 1.0 end if @@ -121,7 +121,9 @@ subroutine curent( lvol, mn, Nt, Nz, iflag, ldItGp ) call invfft( mn, im(1:mn), in(1:mn), efmn(1:mn), ofmn(1:mn), cfmn(1:mn), sfmn(1:mn), & Nt, Nz, Bsupz(1:Ntz,ideriv), Bsupt(1:Ntz,ideriv) ) ! map to real space; - if (Lconstraint .eq. -2) then + if (Lbdybnzero) then + Bsups = zero + else call build_vector_potential(lvol, innout, ideriv, 0) do ii = 1, mn ! loop over Fourier harmonics; 17 May 21; @@ -135,8 +137,6 @@ subroutine curent( lvol, mn, Nt, Nz, iflag, ldItGp ) Nt, Nz, Bsups(1:Ntz,ideriv), Bsups_2(1:Ntz,ideriv)) Bsups(1:Ntz,ideriv) = Bsups(1:Ntz,ideriv) + Bsups_2(1:ntz,ideriv) - else - Bsups = zero end if enddo ! end of do ideriv; 31 Jan 13; @@ -159,17 +159,22 @@ subroutine curent( lvol, mn, Nt, Nz, iflag, ldItGp ) if( iflag.eq. 2 .and. ideriv.lt.0 ) cycle ! derivatives of currents wrt geometry is not required; 20 Jun 14; if( iflag.eq.-1 .and. ideriv.gt.0 ) cycle ! derivatives of currents wrt enclosed toroidal and enclosed poloidal flux are not required; 20 Jun 14; - if (Lconstraint .eq. -2) then - ijreal(1:Ntz) = (Bsups(1:Ntz,ideriv) * guvij(1:Ntz,2,1,0) - Bsupt(1:Ntz,ideriv) * guvij(1:Ntz,2,2,0) + Bsupz(1:Ntz,ideriv) * guvij(1:Ntz,2,3,0) ) / sg(1:Ntz,0) - ijimag(1:Ntz) = (Bsups(1:Ntz,ideriv) * guvij(1:Ntz,1,3,0) - Bsupt(1:Ntz,ideriv) * guvij(1:Ntz,2,3,0) + Bsupz(1:Ntz,ideriv) * guvij(1:Ntz,3,3,0) ) / sg(1:Ntz,0) + if (Lbdybnzero) then + ijreal(1:Ntz) = ( - Bsupt(1:Ntz,ideriv) * guvij(1:Ntz,2,2,0) + Bsupz(1:Ntz,ideriv) * guvij(1:Ntz,2,3,0) ) / sg(1:Ntz,0) + ijimag(1:Ntz) = ( - Bsupt(1:Ntz,ideriv) * guvij(1:Ntz,2,3,0) + Bsupz(1:Ntz,ideriv) * guvij(1:Ntz,3,3,0) ) / sg(1:Ntz,0) else - ijreal(1:Ntz) = ( - Bsupt(1:Ntz,ideriv) * guvij(1:Ntz,2,2,0) + Bsupz(1:Ntz,ideriv) * guvij(1:Ntz,2,3,0) ) / sg(1:Ntz,0) - ijimag(1:Ntz) = ( - Bsupt(1:Ntz,ideriv) * guvij(1:Ntz,2,3,0) + Bsupz(1:Ntz,ideriv) * guvij(1:Ntz,3,3,0) ) / sg(1:Ntz,0) + ijreal(1:Ntz) = ( - Bsupt(1:Ntz,ideriv) * guvij(1:Ntz,2,2,0) + Bsupz(1:Ntz,ideriv) * guvij(1:Ntz,2,3,0) + Bsups(1:Ntz,ideriv) * guvij(1:Ntz,2,1,0)) / sg(1:Ntz,0) + ijimag(1:Ntz) = ( - Bsupt(1:Ntz,ideriv) * guvij(1:Ntz,2,3,0) + Bsupz(1:Ntz,ideriv) * guvij(1:Ntz,3,3,0) + Bsups(1:Ntz,ideriv) * guvij(1:Ntz,1,3,0)) / sg(1:Ntz,0) end if if( ideriv.eq.-1 ) then ! add derivatives of metrics with respect to interface geometry; 15 Sep 16; - ijreal(1:Ntz) = ijreal(1:Ntz) + ( - Bsupt(1:Ntz, 0) * guvij(1:Ntz,2,2,1) + Bsupz(1:Ntz, 0) * guvij(1:Ntz,2,3,1) ) / sg(1:Ntz,0) - ijimag(1:Ntz) = ijimag(1:Ntz) + ( - Bsupt(1:Ntz, 0) * guvij(1:Ntz,2,3,1) + Bsupz(1:Ntz, 0) * guvij(1:Ntz,3,3,1) ) / sg(1:Ntz,0) + if (Lbdybnzero) then + ijreal(1:Ntz) = ijreal(1:Ntz) + ( - Bsupt(1:Ntz, 0) * guvij(1:Ntz,2,2,1) + Bsupz(1:Ntz, 0) * guvij(1:Ntz,2,3,1) ) / sg(1:Ntz,0) + ijimag(1:Ntz) = ijimag(1:Ntz) + ( - Bsupt(1:Ntz, 0) * guvij(1:Ntz,2,3,1) + Bsupz(1:Ntz, 0) * guvij(1:Ntz,3,3,1) ) / sg(1:Ntz,0) + else + ijreal(1:Ntz) = ijreal(1:Ntz) + ( - Bsupt(1:Ntz, 0) * guvij(1:Ntz,2,2,1) + Bsupz(1:Ntz, 0) * guvij(1:Ntz,2,3,1) + Bsups(1:Ntz, 0) * guvij(1:Ntz,2,1,1)) / sg(1:Ntz,0) + ijimag(1:Ntz) = ijimag(1:Ntz) + ( - Bsupt(1:Ntz, 0) * guvij(1:Ntz,2,3,1) + Bsupz(1:Ntz, 0) * guvij(1:Ntz,3,3,1) + Bsups(1:Ntz, 0) * guvij(1:Ntz,1,3,1)) / sg(1:Ntz,0) + end if endif ifail = 0 diff --git a/src/inputlist.f90 b/src/inputlist.f90 index 04c9088c..b9fe19fe 100644 --- a/src/inputlist.f90 +++ b/src/inputlist.f90 @@ -64,6 +64,7 @@ module inputlist !<
      INTEGER :: Lconstraint = -1 !< selects constraints; primarily used in ma02aa() and mp00ac(). !<
        + !<
      • if \c Lconstraint==-2, then in the plasma region \f$\Delta \psi_t\f$ is varied to match the prescriped "linking" current, \c curpol.
        • !<
      • if \c Lconstraint==-1, then in the plasma regions \f$\Delta\psi_t\f$, \f$\mu\f$ and \f$\Delta \psi_p\f$ are *not* varied !< and in the vacuum region (only for free-boundary) \f$\Delta\psi_t\f$ and \f$\Delta \psi_p\f$ are *not* varied, and \f$\mu = 0\f$.
        • !<
      • if \c Lconstraint==0, then in the plasma regions \f$\Delta\psi_t\f$, \f$\mu\f$ and \f$\Delta \psi_p\f$ are *not* varied @@ -184,6 +185,7 @@ module inputlist !<
      • toroidal size is \f$L = 2\pi*\f$\c rtor
        • !<
      INTEGER :: Lreflect = 0 !< =1 reflect the upper and lower bound in slab, =0 do not reflect + LOGICAL :: Lbdybnzero = .true. !< =.true. assume Bsups=0 on boundary (single-volume calculation), =.false. obtain Bsups from Vns and Vnc input REAL :: Rac( 0:MNtor ) = 0.0 !< stellarator symmetric coordinate axis; REAL :: Zas( 0:MNtor ) = 0.0 !< stellarator symmetric coordinate axis; @@ -641,6 +643,7 @@ module inputlist rpol ,& rtor ,& Lreflect ,& + Lbdybnzero ,& Rac ,& Zas ,& Ras ,& @@ -820,6 +823,7 @@ subroutine initialize_inputs mupfits = 8 Lreflect = 0 + Lbdybnzero = .true. ! numericlist diff --git a/src/ma02aa.f90 b/src/ma02aa.f90 index bca9fa2c..08ff0f9c 100644 --- a/src/ma02aa.f90 +++ b/src/ma02aa.f90 @@ -375,6 +375,7 @@ subroutine ma02aa( lvol, NN ) !> e.g. \c C05PCF, can be used to (iteratively) find the appropriately-constrained Beltrami solution, i.e. \f${\bf f}(\boldsymbol{\mu})=0\f$.
      • !>
    • The function \f${\bf f}(\boldsymbol{\mu})\f$, which is computed by mp00ac(), is defined by the input parameter \c Lconstraint: !>
        +!>
      • If \c Lconstraint = -2, then \f$\boldsymbol{\mu}\f$ is *not* varied, but \f$\Delta \psi_t\f$ is, and \c Nxdof=1.
        • !>
      • If \c Lconstraint = -1, 0, then \f$\boldsymbol{\mu}\f$ is *not* varied and \c Nxdof=0.
        • !>
      • If \c Lconstraint = 1, then \f$\boldsymbol{\mu}\f$ is varied to satisfy the transform constraints; !> and \c Nxdof=1 in the simple torus and \c Nxdof=2 in the annular regions. @@ -382,8 +383,7 @@ subroutine ma02aa( lvol, NN ) !> and only \f$\mu=\boldsymbol{\mu}_1\f$ is varied in order to satisfy the transform constraint on the "outer" interface of that volume.)
        • !>
      • \todo If \c Lconstraint = 2, then \f$\mu=\boldsymbol{\mu}_1\f$ is varied in order to satisfy the helicity constraint, !> and \f$\Delta\psi_p=\boldsymbol{\mu}_2\f$ is *not* varied, and \c Nxdof=1. -!> (under re-construction) -!>
      • If \c Lconstraint = -2, then \f$\boldsymbol{\mu}\f$ is *not* varied, but \f$\Delta \psi_t\f$ is, and \c Nxdof=1. Only for single-volume calculations.
        • +!> (under re-construction) !> !> !>
      • diff --git a/src/mp00ac.f90 b/src/mp00ac.f90 index 7935b98d..9b9678e2 100644 --- a/src/mp00ac.f90 +++ b/src/mp00ac.f90 @@ -72,6 +72,7 @@ !>
    • For \c Lcoordinatesingularity=F , the returned function is: !> \f{eqnarray}{ {\bf f}(\mu,\Delta\psi_p) \equiv !> \left\{ \begin{array}{cccccccr} +!> (& 0&,& 0&)^T, & \textrm{if }\texttt{Lconstraint} &=& -2 \\ !> (& 0&,& 0&)^T, & \textrm{if }\texttt{Lconstraint} &=& -1 \\ !> (& 0&,& 0&)^T, & \textrm{if }\texttt{Lconstraint} &=& 0 \\ !> (&{{\,\iota\!\!\!}-}(-1)-\texttt{oita(lvol-1)}&,&{{\,\iota\!\!\!}-}(+1)-\texttt{iota(lvol)}&)^T, & \textrm{if }\texttt{Lconstraint} &=& 1 \\ diff --git a/src/preset.f90 b/src/preset.f90 index d066478f..aeeed3d0 100644 --- a/src/preset.f90 +++ b/src/preset.f90 @@ -300,6 +300,8 @@ subroutine preset endif ! matches if( Lfreebound.eq.1 ) ; + if( Lfreebound.eq.1 .or. .not. Lbdybnzero) then + iVns(ii ) = ( Vns( kk, mm) - Vns(-kk,-mm) ) * jj iBns(ii ) = ( Bns( kk, mm) - Bns(-kk,-mm) ) * jj if( NOTstellsym ) then @@ -310,6 +312,8 @@ subroutine preset iBnc(ii ) = zero endif + endif ! matches if( Lfreebound.eq.1 .or. .not. Lbdybnzero) ; + endif ! end of if( mm.eq.0 .and. nn.eq.0 ) ; enddo ! end of do ii = 1, mn; From 501a652ed9bc6af8e3e93d8601cc1966b115d105 Mon Sep 17 00:00:00 2001 From: Richard Nies Date: Mon, 11 Oct 2021 16:45:12 +0100 Subject: [PATCH 13/20] Removed fossil comment + corrected check of Nvol must be 1 for Lbdybnzero=false --- src/global.f90 | 6 +++--- src/ma02aa.f90 | 4 ---- 2 files changed, 3 insertions(+), 7 deletions(-) diff --git a/src/global.f90 b/src/global.f90 index 19239185..719040a0 100644 --- a/src/global.f90 +++ b/src/global.f90 @@ -1094,14 +1094,14 @@ subroutine check_inputs() cput = GETTIME - write(ounit,1010) cput-cpus, Igeometry, Istellsym, Lreflect + write(ounit,1010) cput-cpus, Igeometry, Istellsym, Lreflect, Lbdybnzero write(ounit,1011) Lfreebound, phiedge, curtor, curpol write(ounit,1012) gamma write(ounit,1013) Nfp, Nvol, Mvol, Mpol, Ntor write(ounit,1014) pscale, Ladiabatic, Lconstraint, mupftol, mupfits write(ounit,1015) Lrad(1:min(Mvol,32)) -1010 format("readin : ",f10.2," : Igeometry=",i3," ; Istellsym=",i3," ; Lreflect="i3" ;") +1010 format("readin : ",f10.2," : Igeometry=",i3," ; Istellsym=",i3," ; Lreflect="i3" ; Lbdybnzero="L2" ;") 1011 format("readin : ", 10x ," : Lfreebound=",i3," ; phiedge="es23.15" ; curtor="es23.15" ; curpol="es23.15" ;") 1012 format("readin : ", 10x ," : gamma="es23.15" ;") 1013 format("readin : ", 10x ," : Nfp=",i3," ; Nvol=",i3," ; Mvol=",i3," ; Mpol=",i3," ; Ntor=",i3," ;") @@ -1129,7 +1129,7 @@ subroutine check_inputs() FATAL( readin, abs(one+gamma).lt.vsmall, 1+gamma appears in denominator in dforce ) !< \todo Please check this; SRH: 27 Feb 18; FATAL( readin, abs(one-gamma).lt.vsmall, 1-gamma appears in denominator in fu00aa ) !< \todo Please check this; SRH: 27 Feb 18; FATAL( readin, Lconstraint.lt.-2 .or. Lconstraint.gt.3, illegal Lconstraint ) - FATAL( readin, Lconstraint.eq.-2 .and. Nvol.ne.1, Lconstraint -2 only for single volume calculation ) + FATAL( readin, .not. Lbdybnzero .and. Lfreebound.eq.1, Lbdybnzero=.false only for fixed-boundary calculation ) FATAL( readin, Igeometry.eq.1 .and. rpol.lt.vsmall, poloidal extent of slab too small or negative ) FATAL( readin, Igeometry.eq.1 .and. rtor.lt.vsmall, toroidal extent of slab too small or negative ) diff --git a/src/ma02aa.f90 b/src/ma02aa.f90 index 08ff0f9c..a22c3e95 100644 --- a/src/ma02aa.f90 +++ b/src/ma02aa.f90 @@ -553,10 +553,6 @@ subroutine ma02aa( lvol, NN ) endif ! end of if( LBlinear ) then; -!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!- - -!latex \subsection{debugging: finite-difference confirmation of the derivatives of the rotational-transform} - !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!- From 21fd29268188486bd89aaa9225722893f769661b Mon Sep 17 00:00:00 2001 From: Richard Nies Date: Mon, 11 Oct 2021 16:47:24 +0100 Subject: [PATCH 14/20] Two checks for Lbdybnzero=.false.: must be fixed-boundary and must be single-volume --- src/global.f90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/global.f90 b/src/global.f90 index 719040a0..593d257b 100644 --- a/src/global.f90 +++ b/src/global.f90 @@ -1130,6 +1130,7 @@ subroutine check_inputs() FATAL( readin, abs(one-gamma).lt.vsmall, 1-gamma appears in denominator in fu00aa ) !< \todo Please check this; SRH: 27 Feb 18; FATAL( readin, Lconstraint.lt.-2 .or. Lconstraint.gt.3, illegal Lconstraint ) FATAL( readin, .not. Lbdybnzero .and. Lfreebound.eq.1, Lbdybnzero=.false only for fixed-boundary calculation ) + FATAL( readin, .not. Lbdybnzero .and. Nvol.ne.1, Lbdybnzero=.false only for single-volume calculation ) FATAL( readin, Igeometry.eq.1 .and. rpol.lt.vsmall, poloidal extent of slab too small or negative ) FATAL( readin, Igeometry.eq.1 .and. rtor.lt.vsmall, toroidal extent of slab too small or negative ) From 4257f89e04b1432c9965dd7aaf6422a223239418 Mon Sep 17 00:00:00 2001 From: Richard Nies Date: Mon, 11 Oct 2021 21:39:46 +0100 Subject: [PATCH 15/20] Checks for testing purpose. --- src/curent.f90 | 4 ++-- src/global.f90 | 2 +- src/mp00ac.f90 | 16 ++++++++++++++-- 3 files changed, 17 insertions(+), 5 deletions(-) diff --git a/src/curent.f90 b/src/curent.f90 index 695c08f3..1f268dee 100644 --- a/src/curent.f90 +++ b/src/curent.f90 @@ -63,7 +63,7 @@ subroutine curent( lvol, mn, Nt, Nz, iflag, ldItGp ) use fileunits, only : ounit - use inputlist, only : Wmacros, Wcurent, Lrad, Lbdybnzero + use inputlist, only : Wmacros, Wcurent, Lrad, Lconstraint, Lbdybnzero use cputiming, only : Tcurent @@ -105,7 +105,7 @@ subroutine curent( lvol, mn, Nt, Nz, iflag, ldItGp ) !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! - if (.not. Lbdybnzero) then + if (Lconstraint .eq. -2) then innout = 1.0 lss = 1.0 end if diff --git a/src/global.f90 b/src/global.f90 index 593d257b..616b1a59 100644 --- a/src/global.f90 +++ b/src/global.f90 @@ -1130,7 +1130,7 @@ subroutine check_inputs() FATAL( readin, abs(one-gamma).lt.vsmall, 1-gamma appears in denominator in fu00aa ) !< \todo Please check this; SRH: 27 Feb 18; FATAL( readin, Lconstraint.lt.-2 .or. Lconstraint.gt.3, illegal Lconstraint ) FATAL( readin, .not. Lbdybnzero .and. Lfreebound.eq.1, Lbdybnzero=.false only for fixed-boundary calculation ) - FATAL( readin, .not. Lbdybnzero .and. Nvol.ne.1, Lbdybnzero=.false only for single-volume calculation ) +! FATAL( readin, .not. Lbdybnzero .and. Nvol.ne.1, Lbdybnzero=.false only for single-volume calculation ) FATAL( readin, Igeometry.eq.1 .and. rpol.lt.vsmall, poloidal extent of slab too small or negative ) FATAL( readin, Igeometry.eq.1 .and. rtor.lt.vsmall, toroidal extent of slab too small or negative ) diff --git a/src/mp00ac.f90 b/src/mp00ac.f90 index 9b9678e2..e09d37d7 100644 --- a/src/mp00ac.f90 +++ b/src/mp00ac.f90 @@ -224,7 +224,14 @@ subroutine mp00ac( Ndof, Xdof, Fdof, Ddof, Ldfjac, iflag ) ! argument list is fi if (Lconstraint .eq. -2) then lmu = mu(lvol) dpf = dpflux(lvol) - dtf = Xdof(1) - xoffset + if (lvol .eq. 1) then + dtf = Xdof(1) - xoffset + else + dtf = dtflux(lvol) + endif + print *,"dtf: ", dtf + print *,"dpf: ", dpf + print *,"lmu: ", lmu end if else ! Lvacuumregion; @@ -613,7 +620,12 @@ subroutine mp00ac( Ndof, Xdof, Fdof, Ddof, Ldfjac, iflag ) ! argument list is fi WCALL( mp00ac, curent,( lvol, mn, Nt, Nz, iflag, dItGpdxtp(0:1,-1:2,lvol)) ) ! Constraint on curpol only - if( iflag.eq.1 ) Fdof(1 ) = dItGpdxtp(1,0,lvol) - curpol +! if( iflag.eq.1 ) Fdof(1 ) = dItGpdxtp(1,0,lvol) - curpol +!TODO + if( iflag.eq.1 ) then + Fdof(1 ) = dItGpdxtp(1,0,lvol) - curpol + print *, "curpol (actual): ", dItGpdxtp(1,0,lvol), ", curpol (target): ", curpol + endif ! Derivative w.r.t. toroidal flux only if( iflag.eq.2 ) Ddof(1,1) = dItGpdxtp(1,1,lvol) From eed9ff77a63c7a90913d3400ec5c7f46ba9ac2a6 Mon Sep 17 00:00:00 2001 From: Richard Nies Date: Tue, 12 Oct 2021 17:42:59 +0100 Subject: [PATCH 16/20] Cleaned up comments. Corrected checks of consistent inputs (Lbdybnzero only allowed with Lconstraint=-2 and Lfreebound=0; and Lconstraint=-2 only allowed with Nvol=1). --- src/global.f90 | 4 ++-- src/inputlist.f90 | 2 +- src/ma02aa.f90 | 2 +- src/mp00ac.f90 | 12 ++---------- 4 files changed, 6 insertions(+), 14 deletions(-) diff --git a/src/global.f90 b/src/global.f90 index 616b1a59..d0c574ca 100644 --- a/src/global.f90 +++ b/src/global.f90 @@ -1129,8 +1129,8 @@ subroutine check_inputs() FATAL( readin, abs(one+gamma).lt.vsmall, 1+gamma appears in denominator in dforce ) !< \todo Please check this; SRH: 27 Feb 18; FATAL( readin, abs(one-gamma).lt.vsmall, 1-gamma appears in denominator in fu00aa ) !< \todo Please check this; SRH: 27 Feb 18; FATAL( readin, Lconstraint.lt.-2 .or. Lconstraint.gt.3, illegal Lconstraint ) - FATAL( readin, .not. Lbdybnzero .and. Lfreebound.eq.1, Lbdybnzero=.false only for fixed-boundary calculation ) -! FATAL( readin, .not. Lbdybnzero .and. Nvol.ne.1, Lbdybnzero=.false only for single-volume calculation ) + FATAL( readin, Lconstraint.eq.-2 .and. Nvol.ne.1, Lconstraint=-2 only for single-volume calculation ) + FATAL( readin, .not. Lbdybnzero .and. (Lconstraint.ne.-2 .or. Lfreebound.eq.1) , Lbdybnzero=.false only for fixed-boundary calculation with Lconstraint=-2 ) FATAL( readin, Igeometry.eq.1 .and. rpol.lt.vsmall, poloidal extent of slab too small or negative ) FATAL( readin, Igeometry.eq.1 .and. rtor.lt.vsmall, toroidal extent of slab too small or negative ) diff --git a/src/inputlist.f90 b/src/inputlist.f90 index b9fe19fe..2ba6595e 100644 --- a/src/inputlist.f90 +++ b/src/inputlist.f90 @@ -185,7 +185,7 @@ module inputlist !<
    • toroidal size is \f$L = 2\pi*\f$\c rtor
      • !<
    INTEGER :: Lreflect = 0 !< =1 reflect the upper and lower bound in slab, =0 do not reflect - LOGICAL :: Lbdybnzero = .true. !< =.true. assume Bsups=0 on boundary (single-volume calculation), =.false. obtain Bsups from Vns and Vnc input + LOGICAL :: Lbdybnzero = .true. !< for fixed-boundary, =.true. assume Bsups=0 on boundary, =.false. obtain Bsups from Vns and Vnc input (currently only for Lconstraint=-2) REAL :: Rac( 0:MNtor ) = 0.0 !< stellarator symmetric coordinate axis; REAL :: Zas( 0:MNtor ) = 0.0 !< stellarator symmetric coordinate axis; diff --git a/src/ma02aa.f90 b/src/ma02aa.f90 index a22c3e95..7b6d803b 100644 --- a/src/ma02aa.f90 +++ b/src/ma02aa.f90 @@ -368,7 +368,7 @@ subroutine ma02aa( lvol, NN ) !> !>
      !> -!>
    • In addition to the enclosed toroidal flux, \f$\Delta \psi_t\f$, which is held constant in the plasma volumes (\c Lconstraint != -2), +!>
    • In addition to the enclosed toroidal flux, \f$\Delta \psi_t\f$, which is held constant in the plasma volumes (for \c Lconstraint != -2), !> the Beltrami field in a given volume is assumed to be parameterized by \f$\mu\f$ and \f$\Delta \psi_p\f$. !> (Note that \f$\Delta \psi_p\f$ is not defined in a torus.)
      • !>
    • These are "packed" into an array, e.g. \f$\boldsymbol{\mu} \equiv (\mu, \Delta\psi_p)^T\f$, so that standard library routines , diff --git a/src/mp00ac.f90 b/src/mp00ac.f90 index e09d37d7..6e024c3a 100644 --- a/src/mp00ac.f90 +++ b/src/mp00ac.f90 @@ -224,14 +224,11 @@ subroutine mp00ac( Ndof, Xdof, Fdof, Ddof, Ldfjac, iflag ) ! argument list is fi if (Lconstraint .eq. -2) then lmu = mu(lvol) dpf = dpflux(lvol) - if (lvol .eq. 1) then + if (lvol .eq. mvol) then dtf = Xdof(1) - xoffset else dtf = dtflux(lvol) endif - print *,"dtf: ", dtf - print *,"dpf: ", dpf - print *,"lmu: ", lmu end if else ! Lvacuumregion; @@ -620,12 +617,7 @@ subroutine mp00ac( Ndof, Xdof, Fdof, Ddof, Ldfjac, iflag ) ! argument list is fi WCALL( mp00ac, curent,( lvol, mn, Nt, Nz, iflag, dItGpdxtp(0:1,-1:2,lvol)) ) ! Constraint on curpol only -! if( iflag.eq.1 ) Fdof(1 ) = dItGpdxtp(1,0,lvol) - curpol -!TODO - if( iflag.eq.1 ) then - Fdof(1 ) = dItGpdxtp(1,0,lvol) - curpol - print *, "curpol (actual): ", dItGpdxtp(1,0,lvol), ", curpol (target): ", curpol - endif + if( iflag.eq.1 ) Fdof(1 ) = dItGpdxtp(1,0,lvol) - curpol ! Derivative w.r.t. toroidal flux only if( iflag.eq.2 ) Ddof(1,1) = dItGpdxtp(1,1,lvol) From ef4aa39af75db078edc01687babbbcbea2fb733d Mon Sep 17 00:00:00 2001 From: Richard Nies Date: Thu, 21 Oct 2021 00:07:00 +0100 Subject: [PATCH 17/20] Added test case. --- InputFiles/TestCases/G3V01Lm2Fi.info | 1 + InputFiles/TestCases/G3V01Lm2Fi.sp | 566 +++++++++++++++++++++++++++ 2 files changed, 567 insertions(+) create mode 100644 InputFiles/TestCases/G3V01Lm2Fi.info create mode 100644 InputFiles/TestCases/G3V01Lm2Fi.sp diff --git a/InputFiles/TestCases/G3V01Lm2Fi.info b/InputFiles/TestCases/G3V01Lm2Fi.info new file mode 100644 index 00000000..defc67a1 --- /dev/null +++ b/InputFiles/TestCases/G3V01Lm2Fi.info @@ -0,0 +1 @@ +test case to verify Lconstraint=-2 (vary toroidal flux to satisfy poloidal linking current constraint) and Lbdybnzero=.false. (non-zero value of Bn on boundary of fixed-boundary calculation). diff --git a/InputFiles/TestCases/G3V01Lm2Fi.sp b/InputFiles/TestCases/G3V01Lm2Fi.sp new file mode 100644 index 00000000..9dcc5d6e --- /dev/null +++ b/InputFiles/TestCases/G3V01Lm2Fi.sp @@ -0,0 +1,566 @@ +! auto generated by SPECNamelist.py 2021-10-21 00:01:17 +&physicslist + igeometry = 3 + istellsym = 1 + lfreebound = 0 + phiedge = 1.0 + curtor = 0 + curpol = 0 + gamma = 0.0 + nfp = 2 + nvol = 1 + mpol = 8 + ntor = 8 + lrad = 10 + tflux = 2.326694073505053 + pflux = 0.0 + helicity = 3.477699658548601 + pscale = 0.0 + ladiabatic = 0 + pressure = 1.0 + adiabatic = 1.0 + mu = 0.0 + lconstraint = -2 + pl = 0, 0 + ql = 0, 0 + pr = 0, 0 + qr = 0, 0 + iota = 0.0, 0.280941793933848 + lp = 0, 0 + lq = 0, 0 + rp = 0, 0 + rq = 0, 0 + oita = 0.0, 0.280941793933848 + mupftol = 1e-16 + mupfits = 16 + rpol = 1.0 + rtor = 1.0 + rbc(-8:8,0) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.01299829262699682, + -0.0049241737074355385, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rbc(-8:8,1) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.004916498659397338, 0.12888068732779726, + 0.42319286094669983, -0.021257759852879716, 0.003995604635266758, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rbc(-8:8,2) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0008142021375225027, -0.16136229488354392, + 0.13259437283052, -0.0012247422810209936, 0.001859980407135166, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rbc(-8:8,3) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rbc(-8:8,4) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rbc(-8:8,5) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rbc(-8:8,6) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rbc(-8:8,7) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rbc(-8:8,8) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zbs(-8:8,0) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.044105780299325124, + 0.004897626703269465, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zbs(-8:8,1) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.008429218722148045, 0.11104611530330084, + -0.6189708354320282, 0.028721518892198167, 0.0012547220152077777, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zbs(-8:8,2) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.012746565240090657, -0.02209424370288384, + -0.04303116340918449, 0.011306578430716792, -0.0003626874352627375, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zbs(-8:8,3) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zbs(-8:8,4) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zbs(-8:8,5) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zbs(-8:8,6) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zbs(-8:8,7) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zbs(-8:8,8) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rbs(-8:8,0) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rbs(-8:8,1) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rbs(-8:8,2) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rbs(-8:8,3) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rbs(-8:8,4) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rbs(-8:8,5) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rbs(-8:8,6) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rbs(-8:8,7) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rbs(-8:8,8) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zbc(-8:8,0) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zbc(-8:8,1) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zbc(-8:8,2) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zbc(-8:8,3) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zbc(-8:8,4) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zbc(-8:8,5) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zbc(-8:8,6) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zbc(-8:8,7) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zbc(-8:8,8) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rwc(-8:8,0) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rwc(-8:8,1) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rwc(-8:8,2) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rwc(-8:8,3) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rwc(-8:8,4) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rwc(-8:8,5) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rwc(-8:8,6) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rwc(-8:8,7) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rwc(-8:8,8) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zws(-8:8,0) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zws(-8:8,1) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zws(-8:8,2) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zws(-8:8,3) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zws(-8:8,4) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zws(-8:8,5) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zws(-8:8,6) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zws(-8:8,7) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zws(-8:8,8) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rws(-8:8,0) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rws(-8:8,1) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rws(-8:8,2) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rws(-8:8,3) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rws(-8:8,4) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rws(-8:8,5) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rws(-8:8,6) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rws(-8:8,7) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rws(-8:8,8) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zwc(-8:8,0) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zwc(-8:8,1) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zwc(-8:8,2) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zwc(-8:8,3) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zwc(-8:8,4) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zwc(-8:8,5) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zwc(-8:8,6) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zwc(-8:8,7) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zwc(-8:8,8) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + vns(-8:8,0) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0032484357678936454, + 0.04759773531023998, 0.07484164088047297, 0.3491068966898682, + 0.07735720204569335, -0.04213990568445082, 0.0, 0.0 + vns(-8:8,1) = 0.0, 0.0, 0.0514342828823717, 0.24622358712752726, 0.001359929358603049, + -0.025834017904898112, -0.06434987281939539, 0.015416556506977486, + -0.1125546708357967, -0.05658841194313129, -0.03786455364343739, + 0.19124227012881556, -0.3645519514104967, -0.006179595178307682, + 0.057362672292951444, 0.0, 0.0 + vns(-8:8,2) = 0.0, 0.0, -0.14481554819901354, -0.42168301273828834, + 0.16233979647432484, 0.31330537429640026, 0.10321400286640497, + 0.05709562946522295, -0.1371688175533019, -0.007672840714461098, + -0.08546794825843257, 0.17723799697567313, 0.15340356518419954, + -0.035080425027089725, -0.02171370713627696, 0.0, 0.0 + vns(-8:8,3) = 0.0, 0.0, 0.15655040702185627, 0.4174834890894755, -0.2642683215937827, + 0.02223199880007334, 0.021430317480327095, 0.07416427163336486, + -0.19931478721188262, -0.06519672896461476, -0.03358927168429972, + -0.26428948720531764, -0.04440912116686778, 0.03456978237671945, + 0.013324200166517247, 0.0, 0.0 + vns(-8:8,4) = 0.0, 0.0, -0.08329033006640713, -0.25694437993929026, + 0.08212352196213508, -0.018307879744856273, -0.01828899126210851, + -0.03882713915994556, -0.0736253412037836, 0.0265586650645998, + 0.1499083616226613, 0.03881004406830072, -0.030764770976077865, + 0.0064636135049393455, 0.006005704806467853, 0.0, 0.0 + vns(-8:8,5) = 0.0, 0.0, 0.008434770759069172, 0.14359133925639994, 0.310005408232767, + -0.5121129424372788, 0.01576343325071893, 0.12870282106583014, + 0.037174669913650815, 0.0007705639677875958, -0.035265435500200705, + 0.01494259944518275, 0.04833979355061212, -0.006121826938287074, + -0.016188729223063716, 0.0, 0.0 + vns(-8:8,6) = 0.0, 0.0, -0.005842569954613564, -0.08871007887173403, + -0.47788806490147845, 0.25892135117550424, 0.045829079344694476, + 0.22170271928799812, 0.11974295868370594, 0.05897964449697521, + -0.015329979047629092, 0.020342062858547474, 0.007038146692507075, + -0.01176393507816326, -0.006930676016976327, 0.0, 0.0 + vns(-8:8,7) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + vns(-8:8,8) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + vnc(-8:8,0) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.0, -0.0, -0.0, + -0.0, -0.0, -0.0, -0.0, 0.0, 0.0 + vnc(-8:8,1) = 0.0, 0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, + -0.0, -0.0, -0.0, -0.0, -0.0, 0.0, 0.0 + vnc(-8:8,2) = 0.0, 0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, + -0.0, -0.0, -0.0, -0.0, -0.0, 0.0, 0.0 + vnc(-8:8,3) = 0.0, 0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, + -0.0, -0.0, -0.0, -0.0, -0.0, 0.0, 0.0 + vnc(-8:8,4) = 0.0, 0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, + -0.0, -0.0, -0.0, -0.0, -0.0, 0.0, 0.0 + vnc(-8:8,5) = 0.0, 0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, + -0.0, -0.0, -0.0, -0.0, -0.0, 0.0, 0.0 + vnc(-8:8,6) = 0.0, 0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, + -0.0, -0.0, -0.0, -0.0, -0.0, 0.0, 0.0 + vnc(-8:8,7) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + vnc(-8:8,8) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + bns(-8:8,0) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + bns(-8:8,1) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + bns(-8:8,2) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + bns(-8:8,3) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + bns(-8:8,4) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + bns(-8:8,5) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + bns(-8:8,6) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + bns(-8:8,7) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + bns(-8:8,8) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + bnc(-8:8,0) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + bnc(-8:8,1) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + bnc(-8:8,2) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + bnc(-8:8,3) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + bnc(-8:8,4) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + bnc(-8:8,5) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + bnc(-8:8,6) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + bnc(-8:8,7) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + bnc(-8:8,8) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rac(0:8) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zas(0:8) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + ras(0:8) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zac(0:8) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + lbdybnzero = .false. +/ + +&numericlist + linitialize = 0 + lautoinitbn = 0 + lzerovac = 0 + ndiscrete = 2 + nquad = -1 + impol = -4 + intor = -4 + lsparse = 0 + lsvdiota = 0 + imethod = 3 + iorder = 2 + iprecon = 1 + iotatol = -1.0 + lextrap = 0 + mregular = -1 + lrzaxis = 2 +/ + +&locallist + lbeltrami = 4 + linitgues = 1 +/ + +&globallist + lfindzero = 2 + escale = 0.0 + opsilon = 1.0 + pcondense = 4.0 + epsilon = 1.0 + wpoloidal = 1.0 + upsilon = 1.0 + forcetol = 1e-14 + c05xmax = 1e-06 + c05xtol = 1e-14 + c05factor = 0.0001 + lreadgf = .true. + mfreeits = 0 + gbntol = 1e-06 + gbnbld = 0.666 + vcasingeps = 1e-12 + vcasingtol = 1e-08 + vcasingits = 8 + vcasingper = 1 +/ + +&diagnosticslist + odetol = 1e-07 + nppts = 0 + nptrj = 8 + lhevalues = .false. + lhevectors = .false. + lhmatrix = .false. + lperturbed = 0 + dpp = -1 + dqq = -1 + lcheck = 0 + ltiming = .false. +/ + +&screenlist +/ + 0 0 5.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 0 1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 0 2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 0 3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 0 4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 0 5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 0 6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 0 7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 0 8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 0 9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 0 10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1 -10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1 -9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1 -8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1 -7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1 -6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1 -5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1 -4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1 -3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1 -2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1 -1 5.000000000000000e-01 5.000000000000000e-01 0.000000000000000e+00 0.000000000000000e+00 + 1 0 1.500000000000000e+00 -1.500000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1 1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1 2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1 3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1 4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1 5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1 6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1 7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1 8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1 9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1 10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 2 -10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 2 -9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 2 -8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 2 -7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 2 -6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 2 -5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 2 -4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 2 -3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 2 -2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 2 -1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 2 0 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 2 1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 2 2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 2 3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 2 4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 2 5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 2 6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 2 7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 2 8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 2 9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 2 10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 3 -10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 3 -9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 3 -8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 3 -7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 3 -6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 3 -5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 3 -4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 3 -3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 3 -2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 3 -1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 3 0 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 3 1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 3 2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 3 3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 3 4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 3 5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 3 6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 3 7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 3 8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 3 9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 3 10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 4 -10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 4 -9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 4 -8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 4 -7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 4 -6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 4 -5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 4 -4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 4 -3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 4 -2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 4 -1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 4 0 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 4 1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 4 2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 4 3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 4 4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 4 5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 4 6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 4 7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 4 8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 4 9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 4 10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 5 -10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 5 -9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 5 -8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 5 -7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 5 -6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 5 -5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 5 -4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 5 -3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 5 -2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 5 -1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 5 0 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 5 1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 5 2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 5 3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 5 4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 5 5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 5 6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 5 7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 5 8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 5 9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 5 10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 6 -10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 6 -9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 6 -8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 6 -7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 6 -6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 6 -5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 6 -4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 6 -3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 6 -2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 6 -1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 6 0 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 6 1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 6 2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 6 3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 6 4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 6 5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 6 6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 6 7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 6 8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 6 9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 6 10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 7 -10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 7 -9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 7 -8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 7 -7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 7 -6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 7 -5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 7 -4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 7 -3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 7 -2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 7 -1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 7 0 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 7 1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 7 2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 7 3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 7 4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 7 5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 7 6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 7 7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 7 8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 7 9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 7 10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 8 -10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 8 -9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 8 -8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 8 -7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 8 -6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 8 -5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 8 -4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 8 -3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 8 -2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 8 -1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 8 0 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 8 1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 8 2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 8 3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 8 4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 8 5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 8 6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 8 7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 8 8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 8 9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 8 10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 9 -10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 9 -9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 9 -8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 9 -7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 9 -6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 9 -5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 9 -4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 9 -3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 9 -2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 9 -1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 9 0 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 9 1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 9 2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 9 3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 9 4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 9 5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 9 6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 9 7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 9 8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 9 9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 9 10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 10 -10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 10 -9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 10 -8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 10 -7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 10 -6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 10 -5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 10 -4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 10 -3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 10 -2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 10 -1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 10 0 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 10 1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 10 2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 10 3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 10 4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 10 5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 10 6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 10 7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 10 8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 10 9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 10 10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 From 3a75635dad7dfc613c0520c40286471d2d694b14 Mon Sep 17 00:00:00 2001 From: Jonathan Schilling Date: Tue, 26 Oct 2021 14:42:50 +0200 Subject: [PATCH 18/20] SCREENLIST awk preprocessing removed from traditional Makefile --- Makefile | 30 ---------- src/global.f90 | 103 +++++++++++++++++++++++++++++++- src/inputlist.f90 | 149 ++++++++++++++++++++++++++++++---------------- 3 files changed, 199 insertions(+), 83 deletions(-) diff --git a/Makefile b/Makefile index ee50ab9c..8e48db22 100644 --- a/Makefile +++ b/Makefile @@ -314,24 +314,12 @@ dspec: $(addsuffix _d.o,$(ALLFILES)) $(MACROS) Makefile # inputlist needs special handling: expansion of DSCREENLIST and NSCREENLIST using awk inputlist_r.o: %_r.o: src/inputlist.f90 $(MACROS) - #@awk -v allfiles='$(ALLFILES)' 'BEGIN{nfiles=split(allfiles,files," ")} \ - #{if($$2=="DSCREENLIST") {for (i=1;i<=nfiles;i++) print " LOGICAL :: W"files[i]" = .false. "}}\ - #{if($$2=="NSCREENLIST") {for (i=1;i<=nfiles;i++) print " W"files[i]" , &"}}\ - #{print}' inputlist.f90 > mnputlist.f90 - #m4 -P $(MACROS) mnputlist.f90 > inputlist_m.F90 - #@rm -f mnputlist.f90 m4 -P $(MACROS) src/inputlist.f90 > inputlist_m.F90 $(FC) $(FLAGS) $(CFLAGS) $(RFLAGS) -o inputlist_r.o -c inputlist_m.F90 $(LIBS) @wc -l -L -w inputlist_m.F90 | awk '{print $$4" has "$$1" lines, "$$2" words, and the longest line is "$$3" characters ;"}' @echo '' inputlist_d.o: %_d.o: src/inputlist.f90 $(MACROS) - #@awk -v allfiles='$(ALLFILES)' 'BEGIN{nfiles=split(allfiles,files," ")} \ - #{if($$2=="DSCREENLIST") {for (i=1;i<=nfiles;i++) print " LOGICAL :: W"files[i]" = .false. "}}\ - #{if($$2=="NSCREENLIST") {for (i=1;i<=nfiles;i++) print " W"files[i]" , &"}}\ - #{print}' inputlist.f90 > mnputlist.f90 - #m4 -P $(MACROS) mnputlist.f90 > inputlist_m.F90 - #@rm -f mnputlist.f90 m4 -P $(MACROS) src/inputlist.f90 > inputlist_m.F90 $(FC) $(FLAGS) $(CFLAGS) $(DFLAGS) -o inputlist_d.o -c inputlist_m.F90 $(LIBS) @wc -l -L -w inputlist_m.F90 | awk '{print $$4" has "$$1" lines, "$$2" words, and the longest line is "$$3" characters ;"}' @@ -341,30 +329,12 @@ inputlist_d.o: %_d.o: src/inputlist.f90 $(MACROS) # global needs special handling: expansion of CPUVARIABLE, BSCREENLIST and WSCREENLIST using awk global_r.o: %_r.o: inputlist_r.o src/global.f90 $(MACROS) - #@awk -v allfiles='$(ALLFILES)' 'BEGIN{nfiles=split(allfiles,files," ")} \ - #{if($$2=="CPUVARIABLE") {for (i=1;i<=nfiles;i++) print " REAL :: T"files[i]" = 0.0, "files[i]"T = 0.0"}}\ - #{if($$2=="DSCREENLIST") {for (i=1;i<=nfiles;i++) print " LOGICAL :: W"files[i]" = .false. "}}\ - #{if($$2=="NSCREENLIST") {for (i=1;i<=nfiles;i++) print " W"files[i]" , &"}}\ - #{if($$2=="BSCREENLIST") {for (i=1;i<=nfiles;i++) print " LlBCAST(W"files[i]",1,0)"}}\ - #{if($$2=="WSCREENLIST") {s="'"'"'" ; d="'"\\\""'" ; for (i=1;i<=nfiles;i++) print " if( W"files[i]" ) write(iunit,"s"("d" W"files[i]" = "d"L1)"s")W"files[i]}}\ - #{print}' src/global.f90 > mlobal.f90 - #m4 -P $(MACROS) mlobal.f90 > global_m.F90 - #@rm -f mlobal.f90 m4 -P $(MACROS) src/global.f90 > global_m.F90 $(FC) $(FLAGS) $(CFLAGS) $(RFLAGS) -o global_r.o -c global_m.F90 $(LIBS) @wc -l -L -w global_m.F90 | awk '{print $$4" has "$$1" lines, "$$2" words, and the longest line is "$$3" characters ;"}' @echo '' global_d.o: %_d.o: inputlist_d.o src/global.f90 $(MACROS) - #@awk -v allfiles='$(ALLFILES)' 'BEGIN{nfiles=split(allfiles,files," ")} \ - #{if($$2=="CPUVARIABLE") {for (i=1;i<=nfiles;i++) print " REAL :: T"files[i]" = 0.0, "files[i]"T = 0.0"}}\ - #{if($$2=="DSCREENLIST") {for (i=1;i<=nfiles;i++) print " LOGICAL :: W"files[i]" = .false. "}}\ - #{if($$2=="NSCREENLIST") {for (i=1;i<=nfiles;i++) print " W"files[i]" , &"}}\ - #{if($$2=="BSCREENLIST") {for (i=1;i<=nfiles;i++) print " LlBCAST(W"files[i]",1,0)"}}\ - #{if($$2=="WSCREENLIST") {s="'"'"'" ; d="'"\\\""'" ; for (i=1;i<=nfiles;i++) print " if( W"files[i]" ) write(iunit,"s"("d" W"files[i]" = "d"L1)"s")W"files[i]}}\ - #{print}' src/global.f90 > mlobal.f90 - #m4 -P $(MACROS) mlobal.f90 > global_m.F90 - #@rm -f mlobal.f90 m4 -P $(MACROS) src/global.f90 > global_m.F90 $(FC) $(FLAGS) $(CFLAGS) $(DFLAGS) -o global_d.o -c global_m.F90 $(LIBS) @wc -l -L -w global_m.F90 | awk '{print $$4" has "$$1" lines, "$$2" words, and the longest line is "$$3" characters ;"}' diff --git a/src/global.f90 b/src/global.f90 index 3980627c..c67b2070 100644 --- a/src/global.f90 +++ b/src/global.f90 @@ -1526,7 +1526,56 @@ subroutine broadcast_inputs if( Wreadin ) then ; cput = GETTIME ; write(ounit,'("readin : ",f10.2," : broadcasting screenlist from ext.sp ;")') cput-cpus endif -! BSCREENLIST ! broadcast screenlist; this is expanded by Makefile; do not remove; + LlBCAST( Wmanual ,1,0) + LlBCAST( Wrzaxis ,1,0) + LlBCAST( Wpackxi ,1,0) + LlBCAST( Wvolume ,1,0) + LlBCAST( Wcoords ,1,0) + LlBCAST( Wbasefn ,1,0) + LlBCAST( Wmemory ,1,0) + LlBCAST( Wmetrix ,1,0) + LlBCAST( Wma00aa ,1,0) + LlBCAST( Wmatrix ,1,0) + LlBCAST( Wspsmat ,1,0) + LlBCAST( Wspsint ,1,0) + LlBCAST( Wmp00ac ,1,0) + LlBCAST( Wma02aa ,1,0) + LlBCAST( Wpackab ,1,0) + LlBCAST( Wtr00ab ,1,0) + LlBCAST( Wcurent ,1,0) + LlBCAST( Wdf00ab ,1,0) + LlBCAST( Wlforce ,1,0) + LlBCAST( Wintghs ,1,0) + LlBCAST( Wmtrxhs ,1,0) + LlBCAST( Wlbpol ,1,0) + LlBCAST( Wbrcast ,1,0) + LlBCAST( Wdfp100 ,1,0) + LlBCAST( Wdfp200 ,1,0) + LlBCAST( Wdforce ,1,0) + LlBCAST( Wnewton ,1,0) + LlBCAST( Wcasing ,1,0) + LlBCAST( Wbnorml ,1,0) + LlBCAST( Wjo00aa ,1,0) + LlBCAST( Wpp00aa ,1,0) + LlBCAST( Wpp00ab ,1,0) + LlBCAST( Wbfield ,1,0) + LlBCAST( Wstzxyz ,1,0) + LlBCAST( Whesian ,1,0) + LlBCAST( Wra00aa ,1,0) + LlBCAST( Wnumrec ,1,0) + LlBCAST( Wdcuhre ,1,0) + LlBCAST( Wminpack,1,0) + LlBCAST( Wiqpack ,1,0) + LlBCAST( Wrksuite,1,0) + LlBCAST( Wi1mach ,1,0) + LlBCAST( Wd1mach ,1,0) + LlBCAST( Wilut ,1,0) + LlBCAST( Witers ,1,0) + LlBCAST( Wsphdf5 ,1,0) + LlBCAST( Wpreset ,1,0) + LlBCAST( Wglobal ,1,0) + LlBCAST( Wxspech ,1,0) + LlBCAST( Wbuild_vector_potential, 1, 0 ) LlBCAST( Wreadin, 1, 0 ) LlBCAST( Wwrtend, 1, 0 ) LlBCAST( Wmacros, 1, 0 ) @@ -1826,7 +1875,57 @@ subroutine wrtend endif write(iunit,'("&screenlist")') -! WSCREENLIST ! write screenlist; this is expanded by Makefile ; do not remove; + if( Wmanual ) write(iunit,'(" Wmanual = ",L1 )') Wmanual + if( Wrzaxis ) write(iunit,'(" Wrzaxis = ",L1 )') Wrzaxis + if( Wpackxi ) write(iunit,'(" Wpackxi = ",L1 )') Wpackxi + if( Wvolume ) write(iunit,'(" Wvolume = ",L1 )') Wvolume + if( Wcoords ) write(iunit,'(" Wcoords = ",L1 )') Wcoords + if( Wbasefn ) write(iunit,'(" Wbasefn = ",L1 )') Wbasefn + if( Wmemory ) write(iunit,'(" Wmemory = ",L1 )') Wmemory + if( Wmetrix ) write(iunit,'(" Wmetrix = ",L1 )') Wmetrix + if( Wma00aa ) write(iunit,'(" Wma00aa = ",L1 )') Wma00aa + if( Wmatrix ) write(iunit,'(" Wmatrix = ",L1 )') Wmatrix + if( Wspsmat ) write(iunit,'(" Wspsmat = ",L1 )') Wspsmat + if( Wspsint ) write(iunit,'(" Wspsint = ",L1 )') Wspsint + if( Wmp00ac ) write(iunit,'(" Wmp00ac = ",L1 )') Wmp00ac + if( Wma02aa ) write(iunit,'(" Wma02aa = ",L1 )') Wma02aa + if( Wpackab ) write(iunit,'(" Wpackab = ",L1 )') Wpackab + if( Wtr00ab ) write(iunit,'(" Wtr00ab = ",L1 )') Wtr00ab + if( Wcurent ) write(iunit,'(" Wcurent = ",L1 )') Wcurent + if( Wdf00ab ) write(iunit,'(" Wdf00ab = ",L1 )') Wdf00ab + if( Wlforce ) write(iunit,'(" Wlforce = ",L1 )') Wlforce + if( Wintghs ) write(iunit,'(" Wintghs = ",L1 )') Wintghs + if( Wmtrxhs ) write(iunit,'(" Wmtrxhs = ",L1 )') Wmtrxhs + if( Wlbpol ) write(iunit,'(" Wlbpol = ",L1 )') Wlbpol + if( Wbrcast ) write(iunit,'(" Wbrcast = ",L1 )') Wbrcast + if( Wdfp100 ) write(iunit,'(" Wdfp100 = ",L1 )') Wdfp100 + if( Wdfp200 ) write(iunit,'(" Wdfp200 = ",L1 )') Wdfp200 + if( Wdforce ) write(iunit,'(" Wdforce = ",L1 )') Wdforce + if( Wnewton ) write(iunit,'(" Wnewton = ",L1 )') Wnewton + if( Wcasing ) write(iunit,'(" Wcasing = ",L1 )') Wcasing + if( Wbnorml ) write(iunit,'(" Wbnorml = ",L1 )') Wbnorml + if( Wjo00aa ) write(iunit,'(" Wjo00aa = ",L1 )') Wjo00aa + if( Wpp00aa ) write(iunit,'(" Wpp00aa = ",L1 )') Wpp00aa + if( Wpp00ab ) write(iunit,'(" Wpp00ab = ",L1 )') Wpp00ab + if( Wbfield ) write(iunit,'(" Wbfield = ",L1 )') Wbfield + if( Wstzxyz ) write(iunit,'(" Wstzxyz = ",L1 )') Wstzxyz + if( Whesian ) write(iunit,'(" Whesian = ",L1 )') Whesian + if( Wra00aa ) write(iunit,'(" Wra00aa = ",L1 )') Wra00aa + if( Wnumrec ) write(iunit,'(" Wnumrec = ",L1 )') Wnumrec + if( Wdcuhre ) write(iunit,'(" Wdcuhre = ",L1 )') Wdcuhre + if( Wminpack ) write(iunit,'(" Wminpack= ",L1 )') Wminpack + if( Wiqpack ) write(iunit,'(" Wiqpack = ",L1 )') Wiqpack + if( Wrksuite ) write(iunit,'(" Wrksuite= ",L1 )') Wrksuite + if( Wi1mach ) write(iunit,'(" Wi1mach = ",L1 )') Wi1mach + if( Wd1mach ) write(iunit,'(" Wd1mach = ",L1 )') Wd1mach + if( Wilut ) write(iunit,'(" Wilut = ",L1 )') Wilut + if( Witers ) write(iunit,'(" Witers = ",L1 )') Witers + if( Wsphdf5 ) write(iunit,'(" Wsphdf5 = ",L1 )') Wsphdf5 + if( Wpreset ) write(iunit,'(" Wpreset = ",L1 )') Wpreset + if( Wglobal ) write(iunit,'(" Wglobal = ",L1 )') Wglobal + if( Wxspech ) write(iunit,'(" Wxspech = ",L1 )') Wxspech + if( Wbuild_vector_potential) write(iunit,'(" Wbuild_vector_potential = ",L1 )') Wbuild_vector_potential + if( Wreadin ) write(iunit,'(" Wreadin = ",L1 )') Wreadin if( Wwrtend ) write(iunit,'(" Wwrtend = ",L1 )') Wwrtend if( Wmacros ) write(iunit,'(" Wmacros = ",L1 )') Wmacros diff --git a/src/inputlist.f90 b/src/inputlist.f90 index 2ba6595e..0406fbbe 100644 --- a/src/inputlist.f90 +++ b/src/inputlist.f90 @@ -548,56 +548,55 @@ module inputlist !> \brief The namelist \c screenlist controls screen output. !> Every subroutine, e.g. \c xy00aa.h, has its own write flag, \c Wxy00aa. !> @{ - LOGICAL :: Wmanual = .false. - LOGICAL :: Wrzaxis = .false. - LOGICAL :: Wpackxi = .false. - LOGICAL :: Wvolume = .false. - LOGICAL :: Wcoords = .false. - LOGICAL :: Wbasefn = .false. - LOGICAL :: Wmemory = .false. - LOGICAL :: Wmetrix = .false. - LOGICAL :: Wma00aa = .false. - LOGICAL :: Wmatrix = .false. - LOGICAL :: Wspsmat = .false. - LOGICAL :: Wspsint = .false. - LOGICAL :: Wmp00ac = .false. - LOGICAL :: Wma02aa = .false. - LOGICAL :: Wpackab = .false. - LOGICAL :: Wtr00ab = .false. - LOGICAL :: Wcurent = .false. - LOGICAL :: Wdf00ab = .false. - LOGICAL :: Wlforce = .false. - LOGICAL :: Wintghs = .false. - LOGICAL :: Wmtrxhs = .false. - LOGICAL :: Wlbpol = .false. - LOGICAL :: Wbrcast = .false. - LOGICAL :: Wdfp100 = .false. - LOGICAL :: Wdfp200 = .false. - LOGICAL :: Wdforce = .false. - LOGICAL :: Wnewton = .false. - LOGICAL :: Wcasing = .false. - LOGICAL :: Wbnorml = .false. - LOGICAL :: Wjo00aa = .false. - LOGICAL :: Wpp00aa = .false. - LOGICAL :: Wpp00ab = .false. - LOGICAL :: Wbfield = .false. - LOGICAL :: Wstzxyz = .false. - LOGICAL :: Whesian = .false. - LOGICAL :: Wra00aa = .false. - LOGICAL :: Wnumrec = .false. - LOGICAL :: Wdcuhre = .false. - LOGICAL :: Wminpack = .false. - LOGICAL :: Wiqpack = .false. - LOGICAL :: Wrksuite = .false. - LOGICAL :: Wi1mach = .false. - LOGICAL :: Wd1mach = .false. - LOGICAL :: Wilut = .false. - LOGICAL :: Witers = .false. - LOGICAL :: Wsphdf5 = .false. - LOGICAL :: Wpreset = .false. - LOGICAL :: Wglobal = .false. - LOGICAL :: Wxspech = .false. -! DSCREENLIST !< define screenlist; this is expanded by Makefile; DO NOT REMOVE; each file compiled by Makefile has its own write flag; + LOGICAL :: Wmanual = .false. + LOGICAL :: Wrzaxis = .false. + LOGICAL :: Wpackxi = .false. + LOGICAL :: Wvolume = .false. + LOGICAL :: Wcoords = .false. + LOGICAL :: Wbasefn = .false. + LOGICAL :: Wmemory = .false. + LOGICAL :: Wmetrix = .false. + LOGICAL :: Wma00aa = .false. + LOGICAL :: Wmatrix = .false. + LOGICAL :: Wspsmat = .false. + LOGICAL :: Wspsint = .false. + LOGICAL :: Wmp00ac = .false. + LOGICAL :: Wma02aa = .false. + LOGICAL :: Wpackab = .false. + LOGICAL :: Wtr00ab = .false. + LOGICAL :: Wcurent = .false. + LOGICAL :: Wdf00ab = .false. + LOGICAL :: Wlforce = .false. + LOGICAL :: Wintghs = .false. + LOGICAL :: Wmtrxhs = .false. + LOGICAL :: Wlbpol = .false. + LOGICAL :: Wbrcast = .false. + LOGICAL :: Wdfp100 = .false. + LOGICAL :: Wdfp200 = .false. + LOGICAL :: Wdforce = .false. + LOGICAL :: Wnewton = .false. + LOGICAL :: Wcasing = .false. + LOGICAL :: Wbnorml = .false. + LOGICAL :: Wjo00aa = .false. + LOGICAL :: Wpp00aa = .false. + LOGICAL :: Wpp00ab = .false. + LOGICAL :: Wbfield = .false. + LOGICAL :: Wstzxyz = .false. + LOGICAL :: Whesian = .false. + LOGICAL :: Wra00aa = .false. + LOGICAL :: Wnumrec = .false. + LOGICAL :: Wdcuhre = .false. + LOGICAL :: Wminpack = .false. + LOGICAL :: Wiqpack = .false. + LOGICAL :: Wrksuite = .false. + LOGICAL :: Wi1mach = .false. + LOGICAL :: Wd1mach = .false. + LOGICAL :: Wilut = .false. + LOGICAL :: Witers = .false. + LOGICAL :: Wsphdf5 = .false. + LOGICAL :: Wpreset = .false. + LOGICAL :: Wglobal = .false. + LOGICAL :: Wxspech = .false. LOGICAL :: Wbuild_vector_potential = .false. !< \todo: what is this? LOGICAL :: Wreadin = .false. !< write screen output of readin() LOGICAL :: Wwrtend = .false. !< write screen output of wrtend() @@ -740,7 +739,55 @@ module inputlist scaling namelist/screenlist/& -! NSCREENLIST ! namelist screenlist; this is expanded by Makefile; DO NOT REMOVE; + Wmanual , & + Wrzaxis , & + Wpackxi , & + Wvolume , & + Wcoords , & + Wbasefn , & + Wmemory , & + Wmetrix , & + Wma00aa , & + Wmatrix , & + Wspsmat , & + Wspsint , & + Wmp00ac , & + Wma02aa , & + Wpackab , & + Wtr00ab , & + Wcurent , & + Wdf00ab , & + Wlforce , & + Wintghs , & + Wmtrxhs , & + Wlbpol , & + Wbrcast , & + Wdfp100 , & + Wdfp200 , & + Wdforce , & + Wnewton , & + Wcasing , & + Wbnorml , & + Wjo00aa , & + Wpp00aa , & + Wpp00ab , & + Wbfield , & + Wstzxyz , & + Whesian , & + Wra00aa , & + Wnumrec , & + Wdcuhre , & + Wminpack , & + Wiqpack , & + Wrksuite , & + Wi1mach , & + Wd1mach , & + Wilut , & + Witers , & + Wsphdf5 , & + Wpreset , & + Wglobal , & + Wxspech , & Wbuild_vector_potential , & Wreadin , & Wwrtend , & From 8184075d849c3fe42bba23098dd99c722e2ecd9f Mon Sep 17 00:00:00 2001 From: Richard Nies Date: Fri, 17 Dec 2021 11:36:26 -0500 Subject: [PATCH 19/20] FFT to read B field and derivatives. --- .../pythontools/py_spec/output/_processing.py | 593 ++++++++++++++++-- Utilities/pythontools/py_spec/output/spec.py | 5 + 2 files changed, 531 insertions(+), 67 deletions(-) diff --git a/Utilities/pythontools/py_spec/output/_processing.py b/Utilities/pythontools/py_spec/output/_processing.py index 43fe0ef2..5a48a23f 100644 --- a/Utilities/pythontools/py_spec/output/_processing.py +++ b/Utilities/pythontools/py_spec/output/_processing.py @@ -1,5 +1,7 @@ import numpy as np +import gc +from numba import jit def get_grid_and_jacobian_and_metric( self, @@ -166,54 +168,220 @@ def metric( return g -def get_B( +def get_B_FFT( self, lvol=0, jacobian=None, sarr=np.linspace(1, 1, 1), - tarr=np.linspace(0, 0, 1), - zarr=np.linspace(0, 0, 1), + Nt=1, + Nz=1 ): - Bcontra_and_der = get_B_and_der(self, lvol, jacobian, sarr, tarr, zarr) - return Bcontra_and_der[0:3] + if jacobian is None: + R, Z, jacobian, g = get_grid_and_jacobian_and_metric( + self, lvol=lvol, sarr=sarr, tarr=tarr, zarr=zarr + ) + # Lrad = s.input.physics.Lrad[lvol] + Ate = self.vector_potential.Ate[lvol] + Aze = self.vector_potential.Aze[lvol] + Ato = self.vector_potential.Ato[lvol] + Azo = self.vector_potential.Azo[lvol] -def get_s_der_B( + mn = Ate.shape[0] + im = self.output.im + in_ = self.output.in_ + Lrad = self.input.physics.Lrad[lvol] + Mpol = self.input.physics.Mpol + Nfp = self.input.physics.Nfp + Ns = len(sarr) + + fac = [] + sbar = (sarr + 1) / 2 + + nax = np.newaxis + + # Zernike polynomial being used + from ._processing import _get_zernike + zernike, dzernike, _ = _get_zernike(sbar, Lrad, Mpol) + + # Obtain Bs + Bs_mn_cos = np.sum( zernike[:, :, im] * (im[nax,:]*Azo.T[nax, :, :] + in_[nax,:]*Ato.T[nax, :, :]) , axis=(1)) + Bs_mn_sin =-np.sum( zernike[:, :, im] * (im[nax,:]*Aze.T[nax, :, :] + in_[nax,:]*Ate.T[nax, :, :]) , axis=(1)) + Bs = invfft_B(Bs_mn_cos, Bs_mn_sin, mn, im, in_, Nfp, Ns, Nt, Nz) + + # Obtain Bt + Bt_mn_cos = -np.sum( dzernike[:, :, im] * Aze.T[nax, :, :], axis=(1)) + Bt_mn_sin = -np.sum( dzernike[:, :, im] * Azo.T[nax, :, :], axis=(1)) + Bt = invfft_B(Bt_mn_cos, Bt_mn_sin, mn, im, in_, Nfp, Ns, Nt, Nz) + + # Obtain Bz + Bz_mn_cos = np.sum( dzernike[:, :, im] * Ate.T[nax, :, :], axis=(1)) + Bz_mn_sin = np.sum( dzernike[:, :, im] * Ato.T[nax, :, :], axis=(1)) + Bz = invfft_B(Bz_mn_cos, Bz_mn_sin, mn, im, in_, Nfp, Ns, Nt, Nz) + + return np.array([Bs, Bt, Bz]) / jacobian + +def get_s_der_B_FFT( self, lvol=0, jacobian=None, sarr=np.linspace(1, 1, 1), - tarr=np.linspace(0, 0, 1), - zarr=np.linspace(0, 0, 1), + Nt=1, + Nz=1 ): - Bcontra_and_der = get_B_and_der(self, lvol, jacobian, sarr, tarr, zarr) - return Bcontra_and_der[3:6] + if jacobian is None: + R, Z, jacobian, g = get_grid_and_jacobian_and_metric( + self, lvol=lvol, sarr=sarr, tarr=tarr, zarr=zarr + ) + # Lrad = s.input.physics.Lrad[lvol] + Ate = self.vector_potential.Ate[lvol] + Aze = self.vector_potential.Aze[lvol] + Ato = self.vector_potential.Ato[lvol] + Azo = self.vector_potential.Azo[lvol] -def get_t_der_B( + mn = Ate.shape[0] + im = self.output.im + in_ = self.output.in_ + Lrad = self.input.physics.Lrad[lvol] + Mpol = self.input.physics.Mpol + Nfp = self.input.physics.Nfp + Ns = len(sarr) + + fac = [] + sbar = (sarr + 1) / 2 + + nax = np.newaxis + + # Zernike polynomial being used + from ._processing import _get_zernike + _, dzernike, d2zernike = _get_zernike(sbar, Lrad, Mpol) + + # Obtain dsBs + dsBs_mn_cos = np.sum( dzernike[:, :, im] * (im[nax,:]*Azo.T[nax, :, :] + in_[nax,:]*Ato.T[nax, :, :]) , axis=(1)) + dsBs_mn_sin =-np.sum( dzernike[:, :, im] * (im[nax,:]*Aze.T[nax, :, :] + in_[nax,:]*Ate.T[nax, :, :]) , axis=(1)) + dsBs = invfft_B(dsBs_mn_cos, dsBs_mn_sin, mn, im, in_, Nfp, Ns, Nt, Nz) + + # Obtain dsBt + dsBt_mn_cos = -np.sum( d2zernike[:, :, im] * Aze.T[nax, :, :], axis=(1)) + dsBt_mn_sin = -np.sum( d2zernike[:, :, im] * Azo.T[nax, :, :], axis=(1)) + dsBt = invfft_B(dsBt_mn_cos, dsBt_mn_sin, mn, im, in_, Nfp, Ns, Nt, Nz) + + # Obtain dsBz + dsBz_mn_cos = np.sum( d2zernike[:, :, im] * Ate.T[nax, :, :], axis=(1)) + dsBz_mn_sin = np.sum( d2zernike[:, :, im] * Ato.T[nax, :, :], axis=(1)) + dsBz = invfft_B(dsBz_mn_cos, dsBz_mn_sin, mn, im, in_, Nfp, Ns, Nt, Nz) + + return np.array([dsBs, dsBt, dsBz]) / jacobian + +def get_t_der_B_FFT( self, lvol=0, jacobian=None, sarr=np.linspace(1, 1, 1), - tarr=np.linspace(0, 0, 1), - zarr=np.linspace(0, 0, 1), + Nt=1, + Nz=1 ): - Bcontra_and_der = get_B_and_der(self, lvol, jacobian, sarr, tarr, zarr) - return Bcontra_and_der[6:9] + if jacobian is None: + R, Z, jacobian, g = get_grid_and_jacobian_and_metric( + self, lvol=lvol, sarr=sarr, tarr=tarr, zarr=zarr + ) -def get_z_der_B( + # Lrad = s.input.physics.Lrad[lvol] + Ate = self.vector_potential.Ate[lvol] + Aze = self.vector_potential.Aze[lvol] + Ato = self.vector_potential.Ato[lvol] + Azo = self.vector_potential.Azo[lvol] + + mn = Ate.shape[0] + im = self.output.im + in_ = self.output.in_ + Lrad = self.input.physics.Lrad[lvol] + Mpol = self.input.physics.Mpol + Nfp = self.input.physics.Nfp + Ns = len(sarr) + + fac = [] + sbar = (sarr + 1) / 2 + + nax = np.newaxis + + # Zernike polynomial being used + from ._processing import _get_zernike + zernike, dzernike, _ = _get_zernike(sbar, Lrad, Mpol) + + # Obtain dtBs + dtBs_mn_sin = -np.sum( zernike[:, :, im] * im[nax,:] * (im[nax,:]*Azo.T[nax, :, :] + in_[nax,:]*Ato.T[nax, :, :]) , axis=(1)) + dtBs_mn_cos = -np.sum( zernike[:, :, im] * im[nax,:] * (im[nax,:]*Aze.T[nax, :, :] + in_[nax,:]*Ate.T[nax, :, :]) , axis=(1)) + dtBs = invfft_B(dtBs_mn_cos, dtBs_mn_sin, mn, im, in_, Nfp, Ns, Nt, Nz) + + # Obtain dtBt + dtBt_mn_sin = np.sum( dzernike[:, :, im] * im[nax,:] * Aze.T[nax, :, :], axis=(1)) + dtBt_mn_cos = -np.sum( dzernike[:, :, im] * im[nax,:] * Azo.T[nax, :, :], axis=(1)) + dtBt = invfft_B(dtBt_mn_cos, dtBt_mn_sin, mn, im, in_, Nfp, Ns, Nt, Nz) + + # Obtain dtBz + dtBz_mn_sin = -np.sum( dzernike[:, :, im] * im[nax,:] * Ate.T[nax, :, :], axis=(1)) + dtBz_mn_cos = np.sum( dzernike[:, :, im] * im[nax,:] * Ato.T[nax, :, :], axis=(1)) + dtBz = invfft_B(dtBz_mn_cos, dtBz_mn_sin, mn, im, in_, Nfp, Ns, Nt, Nz) + + return np.array([dtBs, dtBt, dtBz]) / jacobian + + +def get_z_der_B_FFT( self, lvol=0, jacobian=None, sarr=np.linspace(1, 1, 1), - tarr=np.linspace(0, 0, 1), - zarr=np.linspace(0, 0, 1), + Nt=1, + Nz=1 ): - Bcontra_and_der = get_B_and_der(self, lvol, jacobian, sarr, tarr, zarr) - return Bcontra_and_der[9:12] + if jacobian is None: + R, Z, jacobian, g = get_grid_and_jacobian_and_metric( + self, lvol=lvol, sarr=sarr, tarr=tarr, zarr=zarr + ) + # Lrad = s.input.physics.Lrad[lvol] + Ate = self.vector_potential.Ate[lvol] + Aze = self.vector_potential.Aze[lvol] + Ato = self.vector_potential.Ato[lvol] + Azo = self.vector_potential.Azo[lvol] -def get_B_and_der( + mn = Ate.shape[0] + im = self.output.im + in_ = self.output.in_ + Lrad = self.input.physics.Lrad[lvol] + Mpol = self.input.physics.Mpol + Nfp = self.input.physics.Nfp + Ns = len(sarr) + + fac = [] + sbar = (sarr + 1) / 2 + + nax = np.newaxis + + # Zernike polynomial being used + from ._processing import _get_zernike + zernike, dzernike, _ = _get_zernike(sbar, Lrad, Mpol) + + # Obtain dzBs + dzBs_mn_sin = np.sum( zernike[:, :, im] * in_[nax,:] * (im[nax,:]*Azo.T[nax, :, :] + in_[nax,:]*Ato.T[nax, :, :]) , axis=(1)) + dzBs_mn_cos = np.sum( zernike[:, :, im] * in_[nax,:] * (im[nax,:]*Aze.T[nax, :, :] + in_[nax,:]*Ate.T[nax, :, :]) , axis=(1)) + dzBs = invfft_B(dzBs_mn_cos, dzBs_mn_sin, mn, im, in_, Nfp, Ns, Nt, Nz) + + # Obtain dzBt + dzBt_mn_sin = -np.sum( dzernike[:, :, im] * in_[nax,:] * Aze.T[nax, :, :], axis=(1)) + dzBt_mn_cos = np.sum( dzernike[:, :, im] * in_[nax,:] * Azo.T[nax, :, :], axis=(1)) + dzBt = invfft_B(dzBt_mn_cos, dzBt_mn_sin, mn, im, in_, Nfp, Ns, Nt, Nz) + + # Obtain dzBz + dzBz_mn_sin = np.sum( dzernike[:, :, im] * in_[nax,:] * Ate.T[nax, :, :], axis=(1)) + dzBz_mn_cos = -np.sum( dzernike[:, :, im] * in_[nax,:] * Ato.T[nax, :, :], axis=(1)) + dzBz = invfft_B(dzBz_mn_cos, dzBz_mn_sin, mn, im, in_, Nfp, Ns, Nt, Nz) + + return np.array([dzBs, dzBt, dzBz]) / jacobian + +def get_B( self, lvol=0, jacobian=None, @@ -221,7 +389,6 @@ def get_B_and_der( tarr=np.linspace(0, 0, 1), zarr=np.linspace(0, 0, 1), ): - if jacobian is None: R, Z, jacobian, g = get_grid_and_jacobian_and_metric( self, lvol=lvol, sarr=sarr, tarr=tarr, zarr=zarr @@ -252,7 +419,13 @@ def get_B_and_der( LZernike = self.input.physics.Igeometry > 1 and lvol == 0 - if not LZernike: + if LZernike: + # Zernike polynomial being used + from ._processing import _get_zernike + zernike, dzernike, _ = _get_zernike(sbar, Lrad, Mpol) + else: + # Chebyshev being used + import numpy.polynomial.chebyshev as Cheb # make basis recombination for cheb lcoeff = np.arange(0, Lrad + 1) + 1 @@ -266,6 +439,8 @@ def get_B_and_der( Ato[:, 0] = np.sum(Ato * (-1.0) ** lcoeff[None, :], 1) Azo[:, 0] = np.sum(Azo * (-1.0) ** lcoeff[None, :], 1) + + # Obtain Bs c = ( im[nax, :, nax, nax] * Azo.T[:, :, nax, nax] + in_[nax, :, nax, nax] * Ato.T[:, :, nax, nax] @@ -274,91 +449,361 @@ def get_B_and_der( + in_[nax, :, nax, nax] * Ate.T[:, :, nax, nax] ) * sina[nax, :, :, :] - ct = im[nax, :, nax, nax] * ( - ( - im[nax, :, nax, nax] * Azo.T[:, :, nax, nax] - + in_[nax, :, nax, nax] * Ato.T[:, :, nax, nax] - ) * sina[nax, :, :, :] - ( - im[nax, :, nax, nax] * Aze.T[:, :, nax, nax] - + in_[nax, :, nax, nax] * Ate.T[:, :, nax, nax] - ) * cosa[nax, :, :, :] ) + if LZernike: + Bs = np.sum( zernike[:, :, im, None, None] * c[None, :, :, :, :], axis=(1, 2)) + else: + Bs = np.rollaxis(np.sum(Cheb.chebval(sarr, c), axis=0), 2) - cz = in_[nax, :, nax, nax] * ( ( + # Obtain Bt + c1 = ( + Aze.T[:, :, nax, nax] * cosa[nax, :, :, :] + + Azo.T[:, :, nax, nax] * sina[nax, :, :, :]) + + if LZernike: + Bt = -np.sum( dzernike[:, :, im, None, None] * c1[None, :, :, :, :], axis=(1, 2)) + else: + Bt = -np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c1)), axis=0),2) + + # Obtain Bz + c2 = ( + Ate.T[:, :, nax, nax] * cosa[nax, :, :, :] + + Ato.T[:, :, nax, nax] * sina[nax, :, :, :]) + + if LZernike: + Bz = np.sum( dzernike[:, :, im, None, None] * c2[None, :, :, :, :], axis=(1, 2)) + else: + Bz = np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c2)), axis=0),2) + + return np.array([Bs, Bt, Bz]) / jacobian + + +def get_s_der_B( + self, + lvol=0, + jacobian=None, + sarr=np.linspace(1, 1, 1), + tarr=np.linspace(0, 0, 1), + zarr=np.linspace(0, 0, 1), +): + + if jacobian is None: + R, Z, jacobian, g = get_grid_and_jacobian_and_metric( + self, lvol=lvol, sarr=sarr, tarr=tarr, zarr=zarr + ) + + # Lrad = s.input.physics.Lrad[lvol] + + + Ate = self.vector_potential.Ate[lvol] + Aze = self.vector_potential.Aze[lvol] + Ato = self.vector_potential.Ato[lvol] + Azo = self.vector_potential.Azo[lvol] + + mn = Ate.shape[0] + im = self.output.im + in_ = self.output.in_ + Lrad = self.input.physics.Lrad[lvol] + Mpol = self.input.physics.Mpol + + fac = [] + sbar = (sarr + 1) / 2 + + nax = np.newaxis + # [mn,it,iz] + ang_arg = im[:, nax, nax] * tarr[nax, :, nax] - in_[:, nax, nax] * zarr[nax, nax, :] + cosa = np.cos(ang_arg) + sina = np.sin(ang_arg) + + LZernike = self.input.physics.Igeometry > 1 and lvol == 0 + + if LZernike: + # Zernike polynomial being used + from ._processing import _get_zernike + _, dzernike, d2zernike = _get_zernike(sbar, Lrad, Mpol) + else: + # Chebyshev being used + import numpy.polynomial.chebyshev as Cheb + # make basis recombination for cheb + lcoeff = np.arange(0, Lrad + 1) + 1 + + Ate = Ate / lcoeff[None, :] + Aze = Aze / lcoeff[None, :] + Ato = Ato / lcoeff[None, :] + Azo = Azo / lcoeff[None, :] + + Ate[:, 0] = np.sum(Ate * (-1.0) ** lcoeff[None, :], 1) + Aze[:, 0] = np.sum(Aze * (-1.0) ** lcoeff[None, :], 1) + Ato[:, 0] = np.sum(Ato * (-1.0) ** lcoeff[None, :], 1) + Azo[:, 0] = np.sum(Azo * (-1.0) ** lcoeff[None, :], 1) + + # Obtain dsBs + c = ( im[nax, :, nax, nax] * Azo.T[:, :, nax, nax] + in_[nax, :, nax, nax] * Ato.T[:, :, nax, nax] - ) * sina[nax, :, :, :] + ( + ) * cosa[nax, :, :, :] - ( im[nax, :, nax, nax] * Aze.T[:, :, nax, nax] + in_[nax, :, nax, nax] * Ate.T[:, :, nax, nax] - ) * cosa[nax, :, :, :] ) + ) * sina[nax, :, :, :] + if LZernike: + dsBs = np.sum(dzernike[:, :, im, None, None] * c[None, :, :, :, :], axis=(1, 2)) + else: + dsBs = np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c)), axis=0), 2) + # Obtain dsBt c1 = ( Aze.T[:, :, nax, nax] * cosa[nax, :, :, :] + Azo.T[:, :, nax, nax] * sina[nax, :, :, :]) - - c1t = im[nax, :, nax, nax] * ( - - Aze.T[:, :, nax, nax] * sina[nax, :, :, :] - + Azo.T[:, :, nax, nax] * cosa[nax, :, :, :]) - - c1z = in_[nax, :, nax, nax] * ( - + Aze.T[:, :, nax, nax] * sina[nax, :, :, :] - - Azo.T[:, :, nax, nax] * cosa[nax, :, :, :]) + if LZernike: + dsBt = -np.sum(d2zernike[:, :, im, None, None] * c1[None, :, :, :, :], axis=(1, 2)) + else: + dsBt = -np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c1,m=2)), axis=0),2) + # Obtain dsBz c2 = ( Ate.T[:, :, nax, nax] * cosa[nax, :, :, :] + Ato.T[:, :, nax, nax] * sina[nax, :, :, :]) - c2t = im[nax, :, nax, nax] * ( - - Ate.T[:, :, nax, nax] * sina[nax, :, :, :] - + Ato.T[:, :, nax, nax] * cosa[nax, :, :, :]) + if LZernike: + dsBz = np.sum(d2zernike[:, :, im, None, None] * c2[None, :, :, :, :], axis=(1, 2)) + else: + dsBz = np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c2,m=2)), axis=0),2) - c2z = in_[nax, :, nax, nax] * ( - + Ate.T[:, :, nax, nax] * sina[nax, :, :, :] - - Ato.T[:, :, nax, nax] * cosa[nax, :, :, :]) + return np.array([dsBs, dsBt, dsBz]) / jacobian + + +def get_t_der_B( + self, + lvol=0, + jacobian=None, + sarr=np.linspace(1, 1, 1), + tarr=np.linspace(0, 0, 1), + zarr=np.linspace(0, 0, 1), +): + if jacobian is None: + R, Z, jacobian, g = get_grid_and_jacobian_and_metric( + self, lvol=lvol, sarr=sarr, tarr=tarr, zarr=zarr + ) + + # Lrad = s.input.physics.Lrad[lvol] + + Ate = self.vector_potential.Ate[lvol] + Aze = self.vector_potential.Aze[lvol] + Ato = self.vector_potential.Ato[lvol] + Azo = self.vector_potential.Azo[lvol] + + mn = Ate.shape[0] + im = self.output.im + in_ = self.output.in_ + Lrad = self.input.physics.Lrad[lvol] + Mpol = self.input.physics.Mpol + + fac = [] + sbar = (sarr + 1) / 2 + + nax = np.newaxis + # [mn,it,iz] + ang_arg = im[:, nax, nax] * tarr[nax, :, nax] - in_[:, nax, nax] * zarr[nax, nax, :] + cosa = np.cos(ang_arg) + sina = np.sin(ang_arg) + + LZernike = self.input.physics.Igeometry > 1 and lvol == 0 if LZernike: # Zernike polynomial being used from ._processing import _get_zernike + zernike, dzernike, _ = _get_zernike(sbar, Lrad, Mpol) + else: + # Chebyshev being used + import numpy.polynomial.chebyshev as Cheb + # make basis recombination for cheb + lcoeff = np.arange(0, Lrad + 1) + 1 - zernike, dzernike, d2zernike = _get_zernike(sbar, Lrad, Mpol) + Ate = Ate / lcoeff[None, :] + Aze = Aze / lcoeff[None, :] + Ato = Ato / lcoeff[None, :] + Azo = Azo / lcoeff[None, :] - Bs = np.sum( zernike[:, :, im, None, None] * c[None, :, :, :, :], axis=(1, 2)) - dsBs = np.sum(dzernike[:, :, im, None, None] * c[None, :, :, :, :], axis=(1, 2)) + Ate[:, 0] = np.sum(Ate * (-1.0) ** lcoeff[None, :], 1) + Aze[:, 0] = np.sum(Aze * (-1.0) ** lcoeff[None, :], 1) + Ato[:, 0] = np.sum(Ato * (-1.0) ** lcoeff[None, :], 1) + Azo[:, 0] = np.sum(Azo * (-1.0) ** lcoeff[None, :], 1) + + # Obtain dtBs + ct = im[nax, :, nax, nax] * ( - ( + im[nax, :, nax, nax] * Azo.T[:, :, nax, nax] + + in_[nax, :, nax, nax] * Ato.T[:, :, nax, nax] + ) * sina[nax, :, :, :] - ( + im[nax, :, nax, nax] * Aze.T[:, :, nax, nax] + + in_[nax, :, nax, nax] * Ate.T[:, :, nax, nax] + ) * cosa[nax, :, :, :] ) + + if LZernike: dtBs = np.sum( zernike[:, :, im, None, None] * ct[None, :, :, :, :], axis=(1, 2)) - dzBs = np.sum( zernike[:, :, im, None, None] * cz[None, :, :, :, :], axis=(1, 2)) + else: + dtBs = np.rollaxis(np.sum(Cheb.chebval(sarr, ct), axis=0), 2) - Bt = -np.sum( dzernike[:, :, im, None, None] * c1[None, :, :, :, :], axis=(1, 2)) - dsBt = -np.sum(d2zernike[:, :, im, None, None] * c1[None, :, :, :, :], axis=(1, 2)) + # Obtain dtBt + c1t = im[nax, :, nax, nax] * ( + - Aze.T[:, :, nax, nax] * sina[nax, :, :, :] + + Azo.T[:, :, nax, nax] * cosa[nax, :, :, :]) + + if LZernike: dtBt = -np.sum( dzernike[:, :, im, None, None] * c1t[None, :, :, :, :], axis=(1, 2)) - dzBt = -np.sum( dzernike[:, :, im, None, None] * c1z[None, :, :, :, :], axis=(1, 2)) + else: + dtBt = -np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c1t)), axis=0),2) - Bz = np.sum( dzernike[:, :, im, None, None] * c2[None, :, :, :, :], axis=(1, 2)) - dsBz = np.sum(d2zernike[:, :, im, None, None] * c2[None, :, :, :, :], axis=(1, 2)) + # Obtain dtBz + c2t = im[nax, :, nax, nax] * ( + - Ate.T[:, :, nax, nax] * sina[nax, :, :, :] + + Ato.T[:, :, nax, nax] * cosa[nax, :, :, :]) + + if LZernike: dtBz = np.sum( dzernike[:, :, im, None, None] * c2t[None, :, :, :, :], axis=(1, 2)) - dzBz = np.sum( dzernike[:, :, im, None, None] * c2z[None, :, :, :, :], axis=(1, 2)) + else: + dtBz = np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c2t)), axis=0),2) + + return np.array([dtBs, dtBt, dtBz]) / jacobian + +def get_z_der_B( + self, + lvol=0, + jacobian=None, + sarr=np.linspace(1, 1, 1), + tarr=np.linspace(0, 0, 1), + zarr=np.linspace(0, 0, 1), +): + + if jacobian is None: + R, Z, jacobian, g = get_grid_and_jacobian_and_metric( + self, lvol=lvol, sarr=sarr, tarr=tarr, zarr=zarr + ) + + # Lrad = s.input.physics.Lrad[lvol] + + + Ate = self.vector_potential.Ate[lvol] + Aze = self.vector_potential.Aze[lvol] + Ato = self.vector_potential.Ato[lvol] + Azo = self.vector_potential.Azo[lvol] + + mn = Ate.shape[0] + im = self.output.im + in_ = self.output.in_ + Lrad = self.input.physics.Lrad[lvol] + Mpol = self.input.physics.Mpol + + fac = [] + sbar = (sarr + 1) / 2 + + nax = np.newaxis + # [mn,it,iz] + ang_arg = im[:, nax, nax] * tarr[nax, :, nax] - in_[:, nax, nax] * zarr[nax, nax, :] + cosa = np.cos(ang_arg) + sina = np.sin(ang_arg) + + LZernike = self.input.physics.Igeometry > 1 and lvol == 0 + + if LZernike: + # Zernike polynomial being used + from ._processing import _get_zernike + zernike, dzernike, _ = _get_zernike(sbar, Lrad, Mpol) else: # Chebyshev being used import numpy.polynomial.chebyshev as Cheb + # make basis recombination for cheb + lcoeff = np.arange(0, Lrad + 1) + 1 - Bs = np.rollaxis(np.sum(Cheb.chebval(sarr, c), axis=0), 2) - dsBs = np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c)), axis=0), 2) - dtBs = np.rollaxis(np.sum(Cheb.chebval(sarr, ct), axis=0), 2) + Ate = Ate / lcoeff[None, :] + Aze = Aze / lcoeff[None, :] + Ato = Ato / lcoeff[None, :] + Azo = Azo / lcoeff[None, :] + + Ate[:, 0] = np.sum(Ate * (-1.0) ** lcoeff[None, :], 1) + Aze[:, 0] = np.sum(Aze * (-1.0) ** lcoeff[None, :], 1) + Ato[:, 0] = np.sum(Ato * (-1.0) ** lcoeff[None, :], 1) + Azo[:, 0] = np.sum(Azo * (-1.0) ** lcoeff[None, :], 1) + + # Obtain dzBs + cz = in_[nax, :, nax, nax] * ( ( + im[nax, :, nax, nax] * Azo.T[:, :, nax, nax] + + in_[nax, :, nax, nax] * Ato.T[:, :, nax, nax] + ) * sina[nax, :, :, :] + ( + im[nax, :, nax, nax] * Aze.T[:, :, nax, nax] + + in_[nax, :, nax, nax] * Ate.T[:, :, nax, nax] + ) * cosa[nax, :, :, :] ) + + if LZernike: + dzBs = np.sum( zernike[:, :, im, None, None] * cz[None, :, :, :, :], axis=(1, 2)) + else: dzBs = np.rollaxis(np.sum(Cheb.chebval(sarr, cz), axis=0), 2) - Bt = -np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c1)), axis=0),2) - dsBt = -np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c1,m=2)), axis=0),2) - dtBt = -np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c1t)), axis=0),2) + # Obtain dzBt + c1z = in_[nax, :, nax, nax] * ( + + Aze.T[:, :, nax, nax] * sina[nax, :, :, :] + - Azo.T[:, :, nax, nax] * cosa[nax, :, :, :]) + + if LZernike: + dzBt = -np.sum( dzernike[:, :, im, None, None] * c1z[None, :, :, :, :], axis=(1, 2)) + else: dzBt = -np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c1z)), axis=0),2) - Bz = np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c2)), axis=0),2) - dsBz = np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c2,m=2)), axis=0),2) - dtBz = np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c2t)), axis=0),2) + # Obtain dzBz + c2z = in_[nax, :, nax, nax] * ( + + Ate.T[:, :, nax, nax] * sina[nax, :, :, :] + - Ato.T[:, :, nax, nax] * cosa[nax, :, :, :]) + + if LZernike: + dzBz = np.sum( dzernike[:, :, im, None, None] * c2z[None, :, :, :, :], axis=(1, 2)) + else: dzBz = np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c2z)), axis=0),2) + return np.array([dzBs, dzBt, dzBz]) / jacobian + + +def get_B_and_der_FFT( + self, + lvol=0, + jacobian=None, + sarr=np.linspace(1, 1, 1), + Nt=1, + Nz=1, +): + + if jacobian is None: + R, Z, jacobian, g = get_grid_and_jacobian_and_metric( + self, lvol=lvol, sarr=sarr, tarr=tarr, zarr=zarr) + + Bs, Bt, Bz = get_B_FFT( self, lvol, jacobian, sarr, Nt, Nz) + dsBs, dsBt, dsBz = get_s_der_B_FFT(self, lvol, jacobian, sarr, Nt, Nz) + dtBs, dtBt, dtBz = get_t_der_B_FFT(self, lvol, jacobian, sarr, Nt, Nz) + dzBs, dzBt, dzBz = get_z_der_B_FFT(self, lvol, jacobian, sarr, Nt, Nz) + + Bcontrav_and_der = np.array([Bs, Bt, Bz, dsBs, dsBt, dsBz, dtBs, dtBt, dtBz, dzBs, dzBt, dzBz]) + return Bcontrav_and_der + +def get_B_and_der( + self, + lvol=0, + jacobian=None, + sarr=np.linspace(1, 1, 1), + tarr=np.linspace(0, 0, 1), + zarr=np.linspace(0, 0, 1), +): + + if jacobian is None: + R, Z, jacobian, g = get_grid_and_jacobian_and_metric( + self, lvol=lvol, sarr=sarr, tarr=tarr, zarr=zarr) - Bcontrav_and_der = np.array([Bs, Bt, Bz, dsBs, dsBt, dsBz, dtBs, dtBt, dtBz, dzBs, dzBt, dzBz]) / jacobian + Bs, Bt, Bz = get_B( self, lvol, jacobian, sarr, tarr, zarr) + dsBs, dsBt, dsBz = get_s_der_B(self, lvol, jacobian, sarr, tarr, zarr) + dtBs, dtBt, dtBz = get_t_der_B(self, lvol, jacobian, sarr, tarr, zarr) + dzBs, dzBt, dzBz = get_z_der_B(self, lvol, jacobian, sarr, tarr, zarr) + + Bcontrav_and_der = np.array([Bs, Bt, Bz, dsBs, dsBt, dsBz, dtBs, dtBt, dtBz, dzBs, dzBt, dzBz]) return Bcontrav_and_der @@ -450,3 +895,17 @@ def _get_zernike(sarr, lrad, mpol): d2zernike = d2zernike * 0.25 # to account for the factor of half in sbar = (1+s)/2 return zernike, dzernike, d2zernike + + +def invfft_B(B_mn_cos, B_mn_sin, mn, im, in_, Nfp, Ns, Nt, Nz): + B_mn = np.zeros(shape=(Ns,Nt,Nz), dtype=complex) + + for imn in range(mn): + mm = im[imn] + nn = int(in_[imn]/Nfp) + B_mn[:, mm,-nn] = 0.5*(B_mn_cos[:,imn] - B_mn_sin[:,imn]*1j) + B_mn[:,-mm, nn] = 0.5*(B_mn_cos[:,imn] + B_mn_sin[:,imn]*1j) + B_mn[:,0,0] = 2*B_mn[:,0,0] + B = np.fft.irfft2(B_mn, s=(Nt,Nz))*Nt*Nz + + return B diff --git a/Utilities/pythontools/py_spec/output/spec.py b/Utilities/pythontools/py_spec/output/spec.py index fd1ffa92..c5689e63 100644 --- a/Utilities/pythontools/py_spec/output/spec.py +++ b/Utilities/pythontools/py_spec/output/spec.py @@ -34,6 +34,11 @@ class SPECout: grid, jacobian, metric, + get_B_FFT, + get_s_der_B_FFT, + get_t_der_B_FFT, + get_z_der_B_FFT, + get_B_and_der_FFT, get_B, get_s_der_B, get_t_der_B, From 44a88bd75f96b8ef6516d57284cfb7149aea10a9 Mon Sep 17 00:00:00 2001 From: Richard Nies Date: Fri, 17 Dec 2021 11:48:51 -0500 Subject: [PATCH 20/20] Removed numba package. --- Utilities/pythontools/py_spec/output/_processing.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/Utilities/pythontools/py_spec/output/_processing.py b/Utilities/pythontools/py_spec/output/_processing.py index 5a48a23f..277c83b5 100644 --- a/Utilities/pythontools/py_spec/output/_processing.py +++ b/Utilities/pythontools/py_spec/output/_processing.py @@ -1,8 +1,6 @@ import numpy as np import gc -from numba import jit - def get_grid_and_jacobian_and_metric( self, lvol=0,