Get_Boozer_sinterp

private subroutine Get_Boozer_sinterp(sf, X1_base_in, X2_base_in, LA_base_in, X1_in, X2_in, LA_in)

Builds the boozer transform coordinate theta^B = theta + lambda + iota(s)*nu(s,theta,zeta) zeta^B = zeta +nu(s,theta,zeta)

since in Boozer, the covariant magnetic field components are the current profiles, B = Itor(s) grad theta^B + Ipol(s) grad zeta^B + X grad s = Itor(s) grad (theta+lambda+iotanu) + Ipol(s) grad (zeta + nu) + X grad s = (Itor(1+dlambda/dtheta) + (Itoriota+Ipol)dnu/dtheta) grad theta + (Itor(dlambda/dzeta)+(Itoriota+Ipol)dnu/dzeta) => dnu/dtheta = (B_theta - Itor - Itordlambda/dtheta ) / (Itoriota+Ipol) => dnu/dzeta = (B_zeta - Ipol - Itordlambda/dzeta ) / (Itor*iota+Ipol) There is a integrability condition for nu: d/dzeta(dnu/dtheta)-d/dthet(dnu/dzeta)=d/dzeta(dB_theta/dtheta)-d/dthet(dB_theta/dzeta)=0 which is equivalent to impose J^s=0. now if lambda is recomputed via a projection of J^s=0 onto the same fourier series as nu, the compatibility condition is EXACTLY(!) fullfilled.

!ALLOCATE(LA_s(1:LA_fbase_nyq%modes))

Type Bound

t_sfl_boozer

Arguments

Type	Intent		Name
class(t_sfl_boozer),	intent(inout)	::	sf
class(t_base),	intent(in)	::	X1_base_in
class(t_base),	intent(in)	::	X2_base_in
class(t_base),	intent(in)	::	LA_base_in
real(kind=wp),	intent(in)	::	X1_in(1:X1_base_in%s%nbase,1:X1_base_in%f%modes)
real(kind=wp),	intent(in)	::	X2_in(1:X2_base_in%s%nbase,1:X2_base_in%f%modes)
real(kind=wp),	intent(in)	::	LA_in(1:LA_base_in%s%nbase,1:LA_base_in%f%modes)

Calls

Help

Called by

Help

Source Code

SUBROUTINE Get_Boozer_sinterp(sf,X1_base_in,X2_base_in,LA_base_in,X1_in,X2_in,LA_in)
  ! MODULES
  USE MODgvec_Globals,ONLY: UNIT_stdOut,ProgressBar
  USE MODgvec_base,ONLY: t_base
  USE MODgvec_fbase,ONLY: t_fbase, sin_cos_map
  USE MODgvec_LinAlg
  USE MODgvec_lambda_solve, ONLY:  Lambda_setup_and_solve
  IMPLICIT NONE
    ! INPUT VARIABLES -------------------------!
    CLASS(t_base),INTENT(IN) :: X1_base_in,X2_base_in,LA_base_in   !< base classes for input U_in
    REAL(wp),INTENT(IN):: X1_in(1:X1_base_in%s%nbase,1:X1_base_in%f%modes)
    REAL(wp),INTENT(IN):: X2_in(1:X2_base_in%s%nbase,1:X2_base_in%f%modes)
    REAL(wp),INTENT(IN):: LA_in(1:LA_base_in%s%nbase,1:LA_base_in%f%modes)
    ! OUTPUT VARIABLES -------------------------!
    CLASS(t_sfl_boozer), INTENT(INOUT) :: sf
    ! LOCAL VARIABLES --------------------------!
    INTEGER               :: mn_max(2),mn_nyq(2),irho,iMode,modes,i_mn,mn_IP
    INTEGER               :: nfp
    REAL(wp)              :: spos,dthet_dzeta,dPhids_int,iota_int,dChids_int
    REAL(wp)              :: b_thet,b_zeta
    REAL(wp)              :: detJ,Itor,Ipol,stmp

    REAL(wp)                          ::  X1_s(  1:X1_base_in%f%modes)
    REAL(wp)                          :: dX1ds_s(1:X1_base_in%f%modes)
    REAL(wp)                          ::  X2_s(  1:X2_base_in%f%modes)
    REAL(wp)                          :: dX2ds_s(1:X2_base_in%f%modes)
    REAL(wp),ALLOCATABLE              :: LA_s(:,:)
    REAL(wp),DIMENSION(sf%nu_fbase%modes) :: nu_m,nu_n
    REAL(wp),DIMENSION(sf%nu_fbase%mn_IP) :: Bcov_thet_IP,Bcov_zeta_IP
    REAL(wp),DIMENSION(sf%nu_fbase%mn_IP) :: dLAdthet_IP,dLAdzeta_IP
    REAL(wp),DIMENSION(sf%nu_fbase%mn_IP) :: LA_IP,fm_IP,fn_IP,gam_tt,gam_tz,gam_zz
    REAL(wp),DIMENSION(sf%nu_fbase%mn_IP) :: X1_IP,dX1ds_IP,dX1dthet,dX1dzeta
    REAL(wp),DIMENSION(sf%nu_fbase%mn_IP) :: X2_IP,dX2ds_IP,dX2dthet,dX2dzeta
    TYPE(t_fbase),ALLOCATABLE             :: X1_fbase_nyq
    TYPE(t_fbase),ALLOCATABLE             :: X2_fbase_nyq
    TYPE(t_fbase),ALLOCATABLE             :: LA_fbase_nyq
    ! CODE --------------------------------------------------------------------------------------------------------------------------!
    nfp = X1_base_in%f%nfp
    IF(nfp.NE.sf%nu_fbase%nfp) CALL abort(__STAMP__, &
                   'GET BOOZER ANGLE TRANSFORM,  sf%nu_fbase%nfp not the same as in X1')
    mn_max(1:2)=sf%nu_fbase%mn_max
    mn_nyq(1:2)=sf%nu_fbase%mn_nyq
    SWRITE(UNIT_StdOut,'(A,A,A,I4,3(A,2I6))')'GET BOOZER ANGLE TRANSFORM (', &
                                             TRIM(MERGE('RECOMPUTE      ','USE EQUILIBRIUM',sf%relambda)),' LAMBDA!), nfp=',nfp, &
                                            ', mn_max_in=',LA_base_in%f%mn_max,', mn_max_out=',mn_max,', mn_int=',mn_nyq
    __PERFON('get_boozer')
    __PERFON('init')

    mn_IP        = sf%nu_fbase%mn_IP  !total number of integration points
    modes        = sf%nu_fbase%modes  !number of modes in output
    dthet_dzeta  = sf%nu_fbase%d_thet*sf%nu_fbase%d_zeta !integration weights

    !same base for X1, but with new mn_nyq (for pre-evaluation of basis functions)
    X1_fbase_nyq =  t_fBase(X1_base_in%f%mn_max,  mn_nyq, &
                                  X1_base_in%f%nfp, &
                      sin_cos_map(X1_base_in%f%sin_cos), &
                                  X1_base_in%f%exclude_mn_zero)
    SWRITE(UNIT_StdOut,*)'        ...Init X1_nyq Base Done'

    X2_fbase_nyq =  t_fBase(X2_base_in%f%mn_max,  mn_nyq, &
                                  X2_base_in%f%nfp, &
                      sin_cos_map(X2_base_in%f%sin_cos), &
                                  X2_base_in%f%exclude_mn_zero)
    SWRITE(UNIT_StdOut,*)'        ...Init X2_nyq Base Done'
     IF(.NOT.sf%relambda)THEN
      !same base for lambda, but with new mn_nyq (for pre-evaluation of basis functions)
      LA_fbase_nyq =  t_fBase(LA_base_in%f%mn_max,  mn_nyq, &
                                  LA_base_in%f%nfp, &
                      sin_cos_map(LA_base_in%f%sin_cos), &
                                  LA_base_in%f%exclude_mn_zero)
      ALLOCATE(LA_s(1:LA_base_in%f%modes,sf%nrho))
      SWRITE(UNIT_StdOut,*)'        ...Init LA_nyq Base Done'
    END IF

    !!!ALLOCATE(LA_s(1:LA_fbase_nyq%modes))
    nu_m=0.0_wp; nu_n=0.0_wp

    __PERFOFF('init')


    CALL ProgressBar(0,sf%nrho) !INIT
    DO irho=1,sf%nrho
      __PERFON('eval_data')
      spos=sf%rho_pos(irho)

      dPhids_int  = sf%phiPrime(irho)
      iota_int    = sf%iota(irho)
      dChids_int  = dPhids_int*iota_int

      !interpolate radially
      X1_s(:)    = X1_base_in%s%evalDOF2D_s(spos,X1_base_in%f%modes,      0,X1_in(:,:))
      dX1ds_s(:) = X1_base_in%s%evalDOF2D_s(spos,X1_base_in%f%modes,DERIV_S,X1_in(:,:))

      X2_s(:)    = X2_base_in%s%evalDOF2D_s(spos,X2_base_in%f%modes,      0,X2_in(:,:))
      dX2ds_s(:) = X2_base_in%s%evalDOF2D_s(spos,X2_base_in%f%modes,DERIV_S,X2_in(:,:))

      IF(.NOT.sf%relambda) THEN
        LA_s(:,irho) = LA_base_in%s%evalDOF2D_s(spos,LA_base_in%f%modes,      0,LA_in(:,:))
      END IF

      !evaluate at integration points
      X1_IP    = X1_fbase_nyq%evalDOF_IP(         0, X1_s(  :))
      dX1ds_IP = X1_fbase_nyq%evalDOF_IP(         0,dX1ds_s(:))
      dX1dthet = X1_fbase_nyq%evalDOF_IP(DERIV_THET, X1_s(  :))
      dX1dzeta = X1_fbase_nyq%evalDOF_IP(DERIV_ZETA, X1_s(  :))

      X2_IP    = X2_fbase_nyq%evalDOF_IP(         0, X2_s(  :))
      dX2ds_IP = X2_fbase_nyq%evalDOF_IP(         0,dX2ds_s(:))
      dX2dthet = X2_fbase_nyq%evalDOF_IP(DERIV_THET, X2_s(  :))
      dX2dzeta = X2_fbase_nyq%evalDOF_IP(DERIV_ZETA, X2_s(  :))


      __PERFOFF('eval_data')
      __PERFON('eval_bsub')
      __PERFON('eval_metrics')

  !$OMP PARALLEL DO &
  !$OMP   SCHEDULE(STATIC) DEFAULT(NONE)  &
  !$OMP   PRIVATE(i_mn,detJ)  &
  !$OMP   SHARED(sf,mn_IP,dX1ds_IP,dX2ds_IP,dX1dthet,dX2dthet,dX1dzeta,dX2dzeta,X1_IP,X2_IP,gam_tt,gam_tz,gam_zz)
      !evaluate metrics on (theta,zeta)
      DO i_mn=1,mn_IP

        detJ        =  ( dX1ds_IP(i_mn)*dX2dthet(i_mn) -dX2ds_IP(i_mn)*dX1dthet(i_mn) ) &
                     * sf%hmap%eval_Jh_aux(X1_IP(i_mn),X2_IP(i_mn),sf%hmap_xv(i_mn)) !Jp*Jh
        gam_tt(i_mn)  = sf%hmap%eval_gij_aux(dX1dthet(i_mn),dX2dthet(i_mn),0.0_wp, &
                                              X1_IP  (i_mn), X2_IP(  i_mn),        &
                                             dX1dthet(i_mn),dX2dthet(i_mn),0.0_wp, &
                                             sf%hmap_xv(i_mn) )/detJ   !g_theta,theta
        gam_tz(i_mn)  = sf%hmap%eval_gij_aux(dX1dthet(i_mn),dX2dthet(i_mn),0.0_wp, &
                                              X1_IP  (i_mn), X2_IP(  i_mn),        &
                                             dX1dzeta(i_mn),dX2dzeta(i_mn),1.0_wp, &
                                             sf%hmap_xv(i_mn) )/detJ   !g_zeta,theta
        gam_zz(i_mn)  = sf%hmap%eval_gij_aux(dX1dzeta(i_mn),dX2dzeta(i_mn),1.0_wp, &
                                              X1_IP  (i_mn), X2_IP(  i_mn),        &
                                             dX1dzeta(i_mn),dX2dzeta(i_mn),1.0_wp, &
                                             sf%hmap_xv(i_mn) )/detJ   !g_zeta,zeta
      END DO !i_mn
  !$OMP END PARALLEL DO

      __PERFOFF('eval_metrics')

      IF(sf%relambda)THEN
      __PERFON('new_lambda')
        CALL Lambda_setup_and_solve(sf%nu_fbase,dPhids_int,dchids_int,gam_tt,gam_tz,gam_zz,sf%lambda(:,irho))
        !CALL Lambda_setup_and_solve(LA_fbase_nyq,dPhids_int,dchids_int,gam_tt,gam_tz,gam_zz,LA_s)
        LA_IP(:)    = sf%nu_fbase%evalDOF_IP(         0,sf%lambda(:,irho))
        dLAdthet_IP = sf%nu_fbase%evalDOF_IP(DERIV_THET,sf%lambda(:,irho))
        dLAdzeta_IP = sf%nu_fbase%evalDOF_IP(DERIV_ZETA,sf%lambda(:,irho))
        __PERFOFF('new_lambda')
      ELSE
        LA_IP(:)    = LA_fbase_nyq%evalDOF_IP(         0,LA_s(:,irho))
        dLAdthet_IP = LA_fbase_nyq%evalDOF_IP(DERIV_THET,LA_s(:,irho))
        dLAdzeta_IP = LA_fbase_nyq%evalDOF_IP(DERIV_ZETA,LA_s(:,irho))
      END IF


      Itor=0.0_wp;Ipol=0.0_wp
  !$OMP PARALLEL DO &
  !$OMP   SCHEDULE(STATIC) DEFAULT(NONE)  &
  !$OMP   PRIVATE(i_mn,b_thet,b_zeta)  &
  !$OMP   REDUCTION(+:Itor,Ipol) &
  !$OMP   SHARED(mn_IP,dchids_int,dPhids_int,dLAdzeta_IP,dLAdthet_IP,gam_tt,gam_tz,gam_zz,Bcov_thet_IP,Bcov_zeta_IP)
      !evaluate B_theta,B_zeta (and integrate for currents)
      DO i_mn=1,mn_IP
        b_thet = dchids_int- dPhids_int*dLAdzeta_IP(i_mn)    !b_theta
        b_zeta = dPhids_int*(1.0_wp   + dLAdthet_IP(i_mn))    !b_zeta

        Bcov_thet_IP(i_mn) = (gam_tt(i_mn)*b_thet + gam_tz(i_mn)*b_zeta)
        Bcov_zeta_IP(i_mn) = (gam_tz(i_mn)*b_thet + gam_zz(i_mn)*b_zeta)
        Itor=Itor+Bcov_thet_IP(i_mn)
        Ipol=Ipol+Bcov_zeta_IP(i_mn)
      END DO !i_mn
  !$OMP END PARALLEL DO
      Itor=(Itor/REAL(mn_IP,wp)) !Itor=zero mode of Bcov_thet
      Ipol=(Ipol/REAL(mn_IP,wp)) !Ipol=zero mode of Bcov_thet

  !    Itor=(1.0_wp/REAL(mn_IP,wp))*SUM(Bcov_thet_IP(:))
  !    Ipol=(1.0_wp/REAL(mn_IP,wp))*SUM(Bcov_zeta_IP(:))

      __PERFOFF('eval_bsub')
      __PERFON('project')

      stmp=1.0_wp/(Itor*iota_int+Ipol)
  !$OMP PARALLEL DO        &
  !$OMP   SCHEDULE(STATIC) DEFAULT(NONE) PRIVATE(i_mn)        &
  !$OMP   SHARED(mn_IP,Itor,Ipol,stmp,dLAdthet_IP,Bcov_thet_IP,fm_IP)
      DO i_mn=1,mn_IP
        fm_IP(i_mn)  = (Bcov_thet_IP(i_mn)-Itor-Itor*dLAdthet_IP(i_mn))*stmp
      END DO
  !$OMP END PARALLEL DO

      !projection: only onto base_dthet
      CALL sf%nu_fbase%projectIPtoDOF(.FALSE.,1.0_wp,DERIV_THET,fm_IP(:),nu_m(:))

      IF(sf%nu_fbase%mn_max(2).GT.0) THEN !3D case
  !$OMP PARALLEL DO        &
  !$OMP   SCHEDULE(STATIC) DEFAULT(NONE) PRIVATE(i_mn)        &
  !$OMP   SHARED(mn_IP,Itor,Ipol,stmp,dLAdzeta_IP,Bcov_zeta_IP,fn_IP)
        DO i_mn=1,mn_IP
          fn_IP(i_mn)= (Bcov_zeta_IP(i_mn)-Ipol-Itor*dLAdzeta_IP(i_mn))*stmp
        END DO
  !$OMP END PARALLEL DO

        !projection onto base_dzeta
        CALL sf%nu_fbase%projectIPtoDOF(.FALSE.,1.0_wp,DERIV_ZETA,fn_IP(:),nu_n(:))
      END IF !3D case (n_max >0)

      ! only if n=0, use formula from base_dthet projected G, else use base_dzeta projected G
      DO iMode=1,modes
        ASSOCIATE(m=>sf%nu_fbase%Xmn(1,iMode),n=>sf%nu_fbase%Xmn(2,iMode))
        IF(m.NE.0) nu_m(iMode)=nu_m(iMode)*(dthet_dzeta*sf%nu_fbase%snorm_base(iMode))/REAL(m*m,wp)
        IF(n.NE.0) nu_n(iMode)=nu_n(iMode)*(dthet_dzeta*sf%nu_fbase%snorm_base(iMode))/REAL(n*n,wp)
        IF(n.EQ.0)THEN
          sf%nu(iMode,irho)=nu_m(iMode)
        ELSE
          sf%nu(iMode,irho)=nu_n(iMode)
        END IF
        IF((m.NE.0).AND.(n.NE.0))THEN
          !compare G_mn results:,
          !WRITE(*,"(A,I3,X,A,I3,X,2(A,X,E11.4,X),A,E11.4)")'DEBUG m=',m,'n=',n,'nu_m',nu_m(iMode),'nu_n',nu_n(iMode),'nu_m - nu_n=',nu_m(iMode)-nu_n(iMode)
        END IF
        END ASSOCIATE !m,n
      END DO
      !write(*,*)'DEBUG ===',is



      __PERFOFF('project')
      CALL ProgressBar(irho,sf%nrho)
    END DO !is
    DEALLOCATE(X1_fbase_nyq)
    DEALLOCATE(X2_fbase_nyq)
    IF(.NOT. sf%relambda) THEN
      CALL sf%nu_fbase%change_base(LA_fbase_nyq,2,LA_s,sf%lambda) !save lambda to sfl_boozer structure!
      DEALLOCATE(LA_fbase_nyq) ; DEALLOCATE(LA_s)
    END IF
    SWRITE(UNIT_StdOut,'(A)') '...DONE.'
    __PERFOFF('get_boozer')
  END SUBROUTINE Get_Boozer_sinterp

Get_Boozer_sinterp Subroutine

Contents

Source Code

private subroutine Get_Boozer_sinterp(sf, X1_base_in, X2_base_in, LA_base_in, X1_in, X2_in, LA_in)

Uses

Type Bound

Arguments

Calls

Called by

Source Code