initialize the type fBase maximum mode numbers, number of integration points, type of basis (sin/cos or sin and cos)
!$OMP PARALLEL DO SCHEDULE(STATIC) DEFAULT(SHARED) PRIVATE(i) DO i=1,sf%mn_IP sf%base_IP( i,:)=sf%eval( 0,sf%x_IP(:,i)) sf%base_dthet_IP(i,:)=sf%eval(DERIV_THET,sf%x_IP(:,i)) sf%base_dzeta_IP(i,:)=sf%eval(DERIV_ZETA,sf%x_IP(:,i)) END DO !$OMP END PARALLEL DO
!! 1D BASE
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(t_fBase), | intent(inout) | :: | sf |
self |
||
| integer, | intent(in) | :: | mn_max_in(2) |
maximum mode in m and n |
||
| integer, | intent(in) | :: | mn_nyq_in(2) |
number of integration points |
||
| integer, | intent(in) | :: | nfp_in |
number of field periods |
||
| character(len=8), | intent(in) | :: | sin_cos_in |
can be either only sine: " sin" only cosine: " cos" or full: "sin_cos" |
||
| logical, | intent(in) | :: | exclude_mn_zero_in |
=true: exclude m=n=0 mode in the basis (only important if cos is in basis) |
SUBROUTINE fBase_init( sf, mn_max_in,mn_nyq_in,nfp_in,sin_cos_in,exclude_mn_zero_in) ! MODULES USE MODgvec_Globals, ONLY: myRank, nRanks IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- ! INPUT VARIABLES INTEGER , INTENT(IN ) :: mn_max_in(2) !! maximum mode in m and n INTEGER , INTENT(IN ) :: mn_nyq_in(2) !! number of integration points INTEGER , INTENT(IN ) :: nfp_in !! number of field periods CHARACTER(LEN=8),INTENT(IN ) :: sin_cos_in !! can be either only sine: " _sin_" only cosine: " _cos_" or full: "_sin_cos_" LOGICAL ,INTENT(IN ) :: exclude_mn_zero_in !! =true: exclude m=n=0 mode in the basis (only important if cos is in basis) !----------------------------------------------------------------------------------------------------------------------------------- ! OUTPUT VARIABLES CLASS(t_fBase), INTENT(INOUT) :: sf !! self !----------------------------------------------------------------------------------------------------------------------------------- ! LOCAL VARIABLES INTEGER :: i,iMode,m,n,mIP,nIP,mn_excl,iRank INTEGER :: modes_sin,modes_cos REAL(wp):: mm,nn !=================================================================================================================================== IF(.NOT.test_called)THEN SWRITE(UNIT_stdOut,'(4X,A,2(A,I6," , ",I6),A,I4,A,L2,A)')'INIT fBase type:', & ' mn_max= (',mn_max_in, & ' ), mn_nyq = ',mn_nyq_in, & ' )\n nfp = ',nfp_in, & ' exclude_mn_zero = ',exclude_mn_zero_in, & ' , sin/cos : '//TRIM(sin_cos_in)//' ...' END IF IF(sf%initialized) THEN CALL abort(__STAMP__, & "Trying to reinit fBase!") END IF IF(mn_nyq_in(1).LT.(2*mn_max_in(1)+1)) & CALL abort(__STAMP__, & "error in fBase: mn_nyq in theta should be >= 2*mn_max(1)+1!",mn_nyq_in(1),REAL(mn_max_in(1))) IF(mn_nyq_in(2).LT.(2*mn_max_in(2)+1)) & CALL abort(__STAMP__, & "error in fBase: mn_nyq in zeta should be >= 2*mn_max(2)+1!",mn_nyq_in(2),REAL(mn_max_in(2))) sf%mn_max(1:2) = mn_max_in(1:2) sf%mn_nyq(1:2) = mn_nyq_in(1:2) sf%mn_IP = sf%mn_nyq(1)*sf%mn_nyq(2) sf%nfp = nfp_in sf%exclude_mn_zero = exclude_mn_zero_in IF(INDEX(TRIM(sin_cos_in),TRIM(sin_cos_map(_SIN_))).NE.0) THEN sf%sin_cos = _SIN_ ELSEIF(INDEX(TRIM(sin_cos_in),TRIM(sin_cos_map(_COS_))).NE.0) THEN sf%sin_cos = _COS_ ELSEIF(INDEX(TRIM(sin_cos_in),TRIM(sin_cos_map(_SINCOS_))).NE.0) THEN sf%sin_cos = _SINCOS_ ELSE CALL abort(__STAMP__, & "error in fBase: sin/cos not correctly specified, should be either _SIN_, _COS_ or _SINCOS_") END IF ASSOCIATE(& m_max => sf%mn_max(1) & , n_max => sf%mn_max(2) & , m_nyq => sf%mn_nyq(1) & , n_nyq => sf%mn_nyq(2) & , nfp => sf%nfp & , sin_cos => sf%sin_cos & , modes => sf%modes & , sin_range => sf%sin_range & , cos_range => sf%cos_range & ) mn_excl=MERGE(1,0,sf%exclude_mn_zero) !=1 if exclude=TRUE, =0 if exclude=FALSE ! modes_sin :: m=0: n=1...n_max , m=1...m_max: n=-n_max...n_max. REMARK: for sine, m=0,n=0 is automatically excluded ! modes_cos :: m=0: n=0...n_max , m=1...m_max: n=-n_max...n_max. mn_excl=True will exclude m=n=0 modes_sin= (n_max ) + m_max*(2*n_max+1) modes_cos= (n_max+1)-mn_excl + m_max*(2*n_max+1) SELECT CASE(sin_cos) CASE(_SIN_) modes = modes_sin sin_range(:) = (/0,modes_sin/) cos_range(:) = (/0,0/) !off sf%mTotal1D=m_max+1 CASE(_COS_) modes = modes_cos sin_range(:) = (/0,0/) !off cos_range(:) = (/0,modes_cos/) sf%mTotal1D=m_max+1 CASE(_SINCOS_) modes = modes_sin+modes_cos sin_range(:) = (/0,modes_sin/) cos_range(:) = (/modes_sin,modes/) sf%mTotal1D=2*(m_max+1) END SELECT !MPI DECOMPOSITION OF MODES INTO EQUAL MODE GROUPS, FOR PARALLELIZING COMPUTATIONS DONE SEPARATELY ON EACH MODE ALLOCATE(sf%offset_modes(0:nRanks)) sf%offset_modes(0)=0 DO iRank=0,nRanks-1 sf%offset_modes(iRank+1)=(modes*(iRank+1))/nRanks END DO sf%modes_str = sf%offset_modes(myRank )+1 sf%modes_end = sf%offset_modes(myRank+1) IF(MPIroot)THEN iRank=COUNT((sf%offset_modes(1:nRanks)-sf%offset_modes(0:nRanks-1)) .EQ. 0) IF (iRank.GT.0) THEN WRITE(UNIT_stdout,'(5X,A,I4,A,I4,A)') & 'WARNING: more MPI ranks than number of modes! ', iRank , ' ranks of ' ,nRanks, & ' have no modes associated. This only affects MPI scaling.' END IF END IF CALL fbase_alloc(sf) iMode=0 !first sine then cosine !SINE IF((sin_range(2)-sin_range(1)).EQ.modes_sin)THEN m=0 DO n=1,n_max iMode=iMode+1 sf%Xmn(:,iMode)=(/m,n*nfp/) !include nfp here END DO !n DO m=1,m_max; DO n=-n_max,n_max iMode=iMode+1 sf%Xmn(:,iMode)=(/m,n*nfp/) !include nfp here END DO; END DO !m,n END IF !sin_range>0 sf%mn_zero_mode=-1 !default !COSINE (for _SINCOS_, it comes after sine) IF((cos_range(2)-cos_range(1)).EQ.modes_cos)THEN m=0 IF(mn_excl.EQ.0) sf%mn_zero_mode=iMode+1 DO n=mn_excl,n_max iMode=iMode+1 sf%Xmn(:,iMode)=(/m,n*nfp/) !include nfp here END DO !n DO m=1,m_max; DO n=-n_max,n_max iMode=iMode+1 sf%Xmn(:,iMode)=(/m,n*nfp/) !include nfp here END DO; END DO !m,n END IF !cos_range>0 IF(iMode.NE.modes) CALL abort(__STAMP__,& ' Problem in Xmn ') DO iMode=1,modes m=sf%Xmn(1,iMode) n=sf%Xmn(2,iMode) ! set odd/even/zero for m-modes IF((m.EQ.0))THEN !m=0 IF((n.EQ.0))THEN !n=0 sf%zero_odd_even(iMode)=MN_ZERO ELSE !n /=0 sf%zero_odd_even(iMode)=M_ZERO END IF ELSE !m /=0 IF(MOD(m,2).EQ.0)THEN sf%zero_odd_even(iMode)=M_EVEN ELSE IF(m.EQ.1)THEN sf%zero_odd_even(iMode)=M_ODD_FIRST ELSE sf%zero_odd_even(iMode)=M_ODD END IF END IF END IF ! compute 1/norm, with norm=int_0^2pi int_0^2pi (base_mn(thet,zeta))^2 dthet dzeta IF((m.EQ.0).AND.(n.EQ.0))THEN !m=n=0 (only needed for cos) sf%snorm_base(iMode)=1.0_wp/(TWOPI*TWOPI) !norm=4pi^2 ELSE sf%snorm_base(iMode)=2.0_wp/(TWOPI*TWOPI) !norm=2pi^2 END IF END DO !iMode=1,modes sf%d_thet = TWOPI/REAL(m_nyq,wp) sf%d_zeta = TWOPI/REAL(n_nyq*nfp,wp) DO mIP=1,m_nyq sf%thet_IP(mIP)=(REAL(mIP,wp)-0.5_wp)*sf%d_thet END DO DO nIP=1,n_nyq sf%zeta_IP(nIP)=(REAL(nIP,wp)-0.5_wp)*sf%d_zeta END DO i=0 DO nIP=1,n_nyq DO mIP=1,m_nyq i=i+1 sf%x_IP(1,i)=sf%thet_IP(mIP) sf%x_IP(2,i)=sf%zeta_IP(nIP) END DO !m END DO !n sf%d_zeta = sf%d_zeta*REAL(nfp,wp) ! to get full integral [0,2pi) !! !$OMP PARALLEL DO SCHEDULE(STATIC) DEFAULT(SHARED) PRIVATE(i) !! DO i=1,sf%mn_IP !! sf%base_IP( i,:)=sf%eval( 0,sf%x_IP(:,i)) !! sf%base_dthet_IP(i,:)=sf%eval(DERIV_THET,sf%x_IP(:,i)) !! sf%base_dzeta_IP(i,:)=sf%eval(DERIV_ZETA,sf%x_IP(:,i)) !! END DO !! !$OMP END PARALLEL DO !SIN(m*theta-(n*nfp)*zeta) DO iMode=sin_range(1)+1,sin_range(2) mm=REAL(sf%Xmn(1,iMode),wp) nn=REAL(sf%Xmn(2,iMode),wp) sf%base_IP(:,iMode) = SIN(mm*sf%x_IP(1,:)-nn*sf%x_IP(2,:)) sf%base_dthet_IP(:,iMode)= mm*COS(mm*sf%x_IP(1,:)-nn*sf%x_IP(2,:)) sf%base_dzeta_IP(:,iMode)=-nn*COS(mm*sf%x_IP(1,:)-nn*sf%x_IP(2,:)) END DO !iMode !COS(m*theta-(n*nfp)*zeta) DO iMode=cos_range(1)+1,cos_range(2) mm=REAL(sf%Xmn(1,iMode),wp) nn=REAL(sf%Xmn(2,iMode),wp) sf%base_IP(:,iMode) = COS(mm*sf%x_IP(1,:)-nn*sf%x_IP(2,:)) sf%base_dthet_IP(:,iMode)=-mm*SIN(mm*sf%x_IP(1,:)-nn*sf%x_IP(2,:)) sf%base_dzeta_IP(:,iMode)= nn*SIN(mm*sf%x_IP(1,:)-nn*sf%x_IP(2,:)) END DO !iMode !!!! 1D BASE !sin(mt-nz) = sin(mt)*cos(nz)-cos(mt)*sin(nz) =as(1,mt)*b(1,nz) + as(2,mt)*b(2,nz) !cos(mt-nz) = cos(mt)*cos(nz)+sin(mt)*sin(nz) =ac(1,mt)*b(1,nz) + ac(2,mt)*b(2,nz) IF((sf%sin_cos.EQ._SIN_).OR.(sf%sin_cos.EQ._SINCOS_))THEN !a1s DO m=0,m_max mm=REAL(m,wp) DO mIP=1,m_nyq ASSOCIATE(xm=>sf%thet_IP(mIP)) sf%base1D_IPthet( mIP,1:2,1+m) =(/ SIN(mm*xm), -COS(mm*xm)/) sf%base1D_dthet_IPthet(mIP,1:2,1+m) =(/ mm*COS(mm*xm), mm*SIN(mm*xm)/) END ASSOCIATE END DO END DO END IF IF((sf%sin_cos.EQ._COS_).OR.(sf%sin_cos.EQ._SINCOS_))THEN i=sf%mTotal1D-(sf%mn_max(1)+1) !=offset, =0 if cos, =m_max+1 if sincos DO m=0,m_max mm=REAL(m,wp) DO mIP=1,m_nyq ASSOCIATE(xm=>sf%thet_IP(mIP)) sf%base1D_IPthet( mIP,1:2,i+1+m) =(/ COS(mm*xm), SIN(mm*xm)/) sf%base1D_dthet_IPthet(mIP,1:2,i+1+m) =(/-mm*SIN(mm*xm), mm*COS(mm*xm)/) END ASSOCIATE END DO END DO END IF DO n=-n_max,n_max nn=REAL(n*nfp,wp) DO nIP=1,n_nyq ASSOCIATE(xn=>sf%zeta_IP(nIP)) sf%base1D_IPzeta( 1:2,n,nIP) =(/ COS(nn*xn), SIN(nn*xn)/) sf%base1D_dzeta_IPzeta(1:2,n,nIP) =(/-nn*SIN(nn*xn), nn*COS(nn*xn)/) END ASSOCIATE END DO END DO END ASSOCIATE !sf sf%initialized=.TRUE. IF(.NOT.test_called) THEN SWRITE(UNIT_stdOut,'(4X,A)')'... DONE' CALL fBase_test(sf) END IF END SUBROUTINE fBase_init