find_boozer_angles Subroutine

public subroutine find_boozer_angles(nrho, iota, fbase_in, LA_in, nu_in, tz_dim, tz_boozer, thetzeta_out)

Uses

  • proc~~find_boozer_angles~~UsesGraph proc~find_boozer_angles find_boozer_angles module~modgvec_fbase MODgvec_fBase proc~find_boozer_angles->module~modgvec_fbase module~modgvec_globals MODgvec_Globals proc~find_boozer_angles->module~modgvec_globals module~modgvec_fbase->module~modgvec_globals iso_fortran_env iso_fortran_env module~modgvec_globals->iso_fortran_env

on one flux surface, find for an given list of (thet_j,zeta_j), the corresponding (thet_j,zeta_j) positions, given Here, new boozer angles are theta=theta+Gt(theta,zeta) zeta=zeta+nu(theta,zeta), with Gt=lambda+iotanu and nu periodic functions and zero average and same base Note that in this routine, we will use a 2d root search with a newton method, setting [f1,f2]^T = [thet+A(thet,zeta)-thet=0, zeta+B(thet,zeta)-zeta*=0]^T that includes the derivatives (Jacobian), so that the newton step needs to the solved: -[f1] [ 1+dA/dthet dA/dzeta] [dthet] | | = | | | | -[f2] [ dB/dthet 1+dB/dzeta] [dzeta]

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nrho

number of surfaces, (second dimension of LA_in and nu_in modes)
real(kind=wp), intent(in) :: iota(nrho)

iota at the rho positions.
type(t_fBase), intent(in) :: fbase_in
real(kind=wp), intent(in) :: LA_in(1:fbase_in%modes,nrho)
real(kind=wp), intent(in) :: nu_in(1:fbase_in%modes,nrho)
integer, intent(in) :: tz_dim
real(kind=wp), intent(in) :: tz_boozer(2,tz_dim)
real(kind=wp), intent(out) :: thetzeta_out(2,tz_dim,nrho)

theta,zeta position in original angles, for given boozer angles

Calls

proc~~find_boozer_angles~~CallsGraph proc~find_boozer_angles find_boozer_angles interface~progressbar ProgressBar proc~find_boozer_angles->interface~progressbar proc~fbase_evaldof_xn t_fBase%fBase_evalDOF_xn proc~find_boozer_angles->proc~fbase_evaldof_xn proc~get_booz_newton get_booz_newton proc~find_boozer_angles->proc~get_booz_newton interface~progressbar->interface~progressbar dgemv dgemv proc~fbase_evaldof_xn->dgemv proc~fbase_eval_xn t_fBase%fBase_eval_xn proc~fbase_evaldof_xn->proc~fbase_eval_xn interface~newtonroot2d NewtonRoot2D proc~get_booz_newton->interface~newtonroot2d interface~newtonroot2d->interface~newtonroot2d

Called by

proc~~find_boozer_angles~~CalledByGraph proc~find_boozer_angles find_boozer_angles proc~self_find_boozer_angles t_sfl_boozer%self_find_boozer_angles proc~self_find_boozer_angles->proc~find_boozer_angles proc~self_find_boozer_angles_irho t_sfl_boozer%self_find_boozer_angles_irho proc~self_find_boozer_angles_irho->proc~find_boozer_angles proc~buildtransform_sfl t_transform_sfl%BuildTransform_SFL proc~buildtransform_sfl->proc~self_find_boozer_angles

Source Code

SUBROUTINE find_boozer_angles(nrho,iota,fbase_in,LA_in,nu_in,tz_dim,tz_boozer,thetzeta_out)
! MODULES
USE MODgvec_Globals,ONLY: UNIT_stdOut,PI,ProgressBar,testlevel
USE MODgvec_fbase  ,ONLY: t_fbase
IMPLICIT NONE
  ! INPUT VARIABLES -------------------------!
  INTEGER      ,INTENT(IN) :: nrho   !! number of surfaces, (second dimension  of LA_in and nu_in modes)
  REAL(wp)     ,INTENT(IN) :: iota(nrho)  !! iota at the rho positions.
  TYPE(t_fbase),INTENT(IN) ::fbase_in     !< same basis of lambda and nu
  REAL(wp)     ,INTENT(IN) :: LA_in(1:fbase_in%modes,nrho) !< fourier coefficients of thet*=thet+LA(theta,zeta)+iota*nu(theta,zeta)
  REAL(wp)     ,INTENT(IN) :: nu_in(1:fbase_in%modes,nrho) !< coefficients of zeta*=zeta+nu(theta,zeta)
  INTEGER      ,INTENT(IN) :: tz_dim                !< size of the list of in thetstar,zetastar
  REAL(wp)     ,INTENT(IN) :: tz_boozer(2,tz_dim) !< theta,zeta positions in boozer angle (same for all rho)

  ! OUTPUT VARIABLES -------------------------!
  REAL(wp) ,INTENT(OUT)   :: thetzeta_out(2,tz_dim,nrho)  !! theta,zeta position in original angles, for given boozer angles
  ! LOCAL VARIABLES --------------------------!
  INTEGER               :: irho,j
  REAL(wp)              :: bounds(2),x0(2)
  REAL(wp)              :: check(tz_dim),maxerr(2,nrho)
  REAL(wp)              :: Gt(1:fbase_in%modes)  !! transform in theta: lambda + iota*nu
  LOGICAL               :: docheck
  ! CODE --------------------------------------------------------------------------------------------------------------------------!
  __PERFON('find_boozer_angles')
  SWRITE(UNIT_StdOut,'(A,2(I8,A))')'Find boozer angles via 2D Newton on  nrho=',nrho,' times ntheta_zeta= ',tz_dim, " points"
  docheck=(testlevel.GT.0)
  bounds=(/PI, PI/fbase_in%nfp/)
  CALL ProgressBar(0,nrho)!init
  DO irho=1,nrho
    Gt(:)=LA_in(:,irho)+iota(irho)*nu_in(:,irho)
!$OMP PARALLEL DO SCHEDULE(STATIC) DEFAULT(NONE) &
!$OMP   PRIVATE(j,x0) FIRSTPRIVATE(irho,bounds) &
!$OMP   SHARED(tz_dim,tz_boozer,thetzeta_out,fbase_in,Gt,nu_in)
    DO j=1,tz_dim
        x0=tz_boozer(:,j)
        thetzeta_out(:,j,irho)=get_booz_newton(x0,bounds,fbase_in,Gt,nu_in(:,irho))
    END DO !j
!$OMP END PARALLEL DO
    CALL ProgressBar(irho,nrho)
  END DO !irho

  IF(docheck)THEN
    DO irho=1,nrho
      check=fbase_in%evalDOF_xn(tz_dim,thetzeta_out(:,:,irho),0,(LA_in(:,irho)+iota(irho)*nu_in(:,irho)))
      maxerr(1,irho)=maxval(abs(check+(thetzeta_out(1,:,irho)-tz_boozer(1,:))))
      check=fbase_in%evalDOF_xn(tz_dim,thetzeta_out(:,:,irho),0,nu_in(:,irho))
      maxerr(2,irho)=maxval(abs(check+(thetzeta_out(2,:,irho)-tz_boozer(2,:))))
    END DO

    IF(ANY(maxerr(:,:).GT.1.0e-12)) THEN
        WRITE(UNIT_stdout,*)'CHECK BOOZER THETA*',MAXVAL(maxerr(1,:))
        WRITE(UNIT_stdout,*)'CHECK BOOZER ZETA*', maxerr
        CALL abort(__STAMP__, &
        "find_boozer_angles: Error in transform")
    END IF
  END IF  !docheck
  SWRITE(UNIT_StdOut,'(A)') '...DONE.'
__PERFOFF('find_boozer_angles')

END SUBROUTINE find_boozer_angles