fBase_init – GVEC

private subroutine fBase_init(sf, mn_max_in, mn_nyq_in, nfp_in, sin_cos_in, exclude_mn_zero_in)

Uses

- MODgvec_Globals
Help

Graph Key

Nodes of different colours represent the following:

Solid arrows point from a submodule to the (sub)module which it is descended from. Dashed arrows point from a module or program unit to modules which it uses.

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(2,sf%x_IP(:,i)) sf%base_dzeta_IP(i,:)=sf%eval(3,sf%x_IP(:,i)) END DO !$OMP END PARALLEL DO

!! 1D BASE

Type Bound

t_fBase

Arguments

Type	Intent		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)

Calls

Help

Called by

Help

Source Code

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
  ! CODE --------------------------------------------------------------------------------------------------------------------------!
  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)),&
        TypeInfo="InvalidParameterError")
  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)),&
        TypeInfo="InvalidParameterError")

  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_",&
        TypeInfo="InvalidParameterError")
  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

fBase_init Subroutine

Contents

Source Code

private subroutine fBase_init(sf, mn_max_in, mn_nyq_in, nfp_in, sin_cos_in, exclude_mn_zero_in)

Uses

Type Bound

Arguments

Calls

Called by

Source Code