InitVMEC Subroutine

public subroutine InitVMEC()

Uses

  • proc~~initvmec~~UsesGraph proc~initvmec InitVMEC module~modgvec_cubic_spline MODgvec_cubic_spline proc~initvmec->module~modgvec_cubic_spline module~modgvec_globals MODgvec_Globals proc~initvmec->module~modgvec_globals module~modgvec_readintools MODgvec_ReadInTools proc~initvmec->module~modgvec_readintools module~modgvec_rprofile_bspl MODgvec_rProfile_bspl proc~initvmec->module~modgvec_rprofile_bspl module~modgvec_vmec_readin MODgvec_VMEC_Readin proc~initvmec->module~modgvec_vmec_readin module~modgvec_vmec_vars MODgvec_VMEC_Vars proc~initvmec->module~modgvec_vmec_vars module~modgvec_cubic_spline->module~modgvec_globals module~sll_m_bsplines sll_m_bsplines module~modgvec_cubic_spline->module~sll_m_bsplines iso_fortran_env iso_fortran_env module~modgvec_globals->iso_fortran_env module~modgvec_readintools->module~modgvec_globals module~modgvec_rprofile_bspl->module~modgvec_globals module~modgvec_rprofile_base MODgvec_rProfile_base module~modgvec_rprofile_bspl->module~modgvec_rprofile_base module~modgvec_rprofile_bspl->module~sll_m_bsplines module~modgvec_vmec_readin->module~modgvec_globals module~modgvec_vmec_vars->module~modgvec_cubic_spline module~modgvec_vmec_vars->module~modgvec_globals module~modgvec_vmec_vars->module~modgvec_rprofile_base module~modgvec_rprofile_base->module~modgvec_globals module~sll_m_assert sll_m_assert module~sll_m_bsplines->module~sll_m_assert module~sll_m_bsplines_base sll_m_bsplines_base module~sll_m_bsplines->module~sll_m_bsplines_base module~sll_m_bsplines_non_uniform sll_m_bsplines_non_uniform module~sll_m_bsplines->module~sll_m_bsplines_non_uniform module~sll_m_bsplines_uniform sll_m_bsplines_uniform module~sll_m_bsplines->module~sll_m_bsplines_uniform module~sll_m_errors sll_m_errors module~sll_m_bsplines->module~sll_m_errors module~sll_m_working_precision sll_m_working_precision module~sll_m_bsplines->module~sll_m_working_precision module~sll_m_bsplines_base->module~sll_m_assert module~sll_m_bsplines_base->module~sll_m_working_precision module~sll_m_bsplines_non_uniform->module~sll_m_assert module~sll_m_bsplines_non_uniform->module~sll_m_bsplines_base module~sll_m_bsplines_non_uniform->module~sll_m_working_precision module~sll_m_bsplines_uniform->module~sll_m_assert module~sll_m_bsplines_uniform->module~sll_m_bsplines_base module~sll_m_bsplines_uniform->module~sll_m_errors module~sll_m_bsplines_uniform->module~sll_m_working_precision module~sll_m_errors->iso_fortran_env

Initialize VMEC module

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! transform VMEC data (R,phi=zeta,Z) to GVEC right hand side system (R,Z,phi), swap sign of zeta !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

Arguments

None

Calls

proc~~initvmec~~CallsGraph proc~initvmec InitVMEC getint getint proc~initvmec->getint getstr getstr proc~initvmec->getstr interface~getfreeunit GETFREEUNIT proc~initvmec->interface~getfreeunit interface~interpolate_cubic_spline interpolate_cubic_spline proc~initvmec->interface~interpolate_cubic_spline proc~fitspline FitSpline proc~initvmec->proc~fitspline proc~fitsplinehalf FitSplineHalf proc~initvmec->proc~fitsplinehalf proc~readvmec ReadVMEC proc~initvmec->proc~readvmec proc~rprofile_eval_at_rho c_rProfile%rProfile_eval_at_rho proc~initvmec->proc~rprofile_eval_at_rho interface~getfreeunit->interface~getfreeunit interface~interpolate_cubic_spline->interface~interpolate_cubic_spline proc~cubspl_eval t_cubspl%cubspl_eval proc~fitsplinehalf->proc~cubspl_eval proc~readnemec ReadNEMEC proc~readvmec->proc~readnemec eval_at_rho2 eval_at_rho2 proc~rprofile_eval_at_rho->eval_at_rho2 proc~rho2_derivative rho2_derivative proc~rprofile_eval_at_rho->proc~rho2_derivative proc~rprofile_drho2 c_rProfile%rProfile_drho2 proc~rprofile_eval_at_rho->proc~rprofile_drho2 proc~rprofile_drho3 c_rProfile%rProfile_drho3 proc~rprofile_eval_at_rho->proc~rprofile_drho3 proc~rprofile_drho4 c_rProfile%rProfile_drho4 proc~rprofile_eval_at_rho->proc~rprofile_drho4 eval_basis eval_basis proc~cubspl_eval->eval_basis eval_basis_and_n_derivs eval_basis_and_n_derivs proc~cubspl_eval->eval_basis_and_n_derivs proc~readnemec->interface~getfreeunit proc~alloc_all alloc_all proc~readnemec->proc~alloc_all proc~halftofull HalfToFull proc~readnemec->proc~halftofull proc~poly_derivative_prefactor poly_derivative_prefactor proc~rho2_derivative->proc~poly_derivative_prefactor proc~rprofile_drho2->eval_at_rho2 proc~rprofile_drho2->proc~rho2_derivative proc~rprofile_drho3->eval_at_rho2 proc~rprofile_drho3->proc~rho2_derivative proc~rprofile_drho4->eval_at_rho2 proc~rprofile_drho4->proc~rho2_derivative proc~halftofull->proc~cubspl_eval

Source Code

SUBROUTINE InitVMEC
! MODULES
USE MODgvec_Globals,ONLY:UNIT_stdOut,abort,fmt_sep,GETFREEUNIT
USE MODgvec_rProfile_bspl, ONLY: t_rProfile_bspl
USE MODgvec_cubic_spline, ONLY: interpolate_cubic_spline
USE MODgvec_ReadInTools
USE MODgvec_VMEC_Vars
USE MODgvec_VMEC_Readin
IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT/OUTPUT VARIABLES
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
INTEGER              :: iMode,bc_unit
LOGICAL              :: useFilter
REAL(wp),ALLOCATABLE :: c_iota(:),knots_iota(:), c_pres(:), knots_pres(:)
REAL(wp),ALLOCATABLE :: c_phi(:),knots_phi(:), c_chi(:), knots_chi(:)
!===================================================================================================================================
IF(.NOT.MPIroot) RETURN
WRITE(UNIT_stdOut,'(A)')'  INIT VMEC INPUT ...'

!VMEC "wout*.nc"  file
VMECdataFile   = GETSTR("VMECwoutfile")
VMECFile_Format= GETINT("VMECwoutfile_format",Proposal=0)

CALL ReadVmec(VMECdataFile,VMECfile_format)

switchzeta=.TRUE. !always switch zeta=-phi_vmec
switchtheta=(signgs>0) !

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! transform VMEC data (R,phi=zeta,Z) to GVEC right hand side system (R,Z,phi), swap sign of zeta
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
WRITE(UNIT_stdOut,'(A)')'  ... switching from VMECs (R,phi,Z) to gvecs (R,Z,phi) coordinate system ...'
IF(.NOT.switchtheta)THEN
  WRITE(UNIT_stdOut,'(A)')'  ... VMEC signgs<0, so zeta direction is switched.'
  DO iMode=1,mn_mode
    IF(NINT(xm(iMode)).EQ.0)THEN
      !xn for m=0 are only positive, so we need to swap the sign of sinus coefficients
      ! -coef_{0n} sin(-n*(-zeta))= coef_{0n} sin(-n*zeta)
      !( cosine cos(-n*(-zeta))=cos(-n*zeta), symmetric in zeta for m=0)
      zmns(iMode,:)=-zmns(iMode,:)
      lmns(iMode,:)=-lmns(iMode,:)
      IF(lasym) THEN
        rmns(iMode,:)=-rmns(iMode,:)
      END IF
    ELSE
      !for m>0 , xn are always pairs of negative and positive modes,
      !so here we can simply swap the sign of the mode number
      xn(iMode)=-xn(iMode)
    END IF
  END DO !iMode=1,mn_mode
  !additional nyq data (bmnc,gmnc
  DO iMode=1,mn_mode_nyq
    ! nyq data is only cosine (bmnc,gmnc)
    IF(NINT(xm_nyq(iMode)).NE.0)THEN
      xn_nyq(iMode)=-xn_nyq(iMode)
    END If
  END DO !iMode=1,mn_mode_nyq
  !since sign of zeta has changed, swap sign of Jacobian, too.
  gmnc=-gmnc
  IF(lasym) gmns=-gmns

  ! also iota must change sign, since its sign depend on the coordinate system
  iotaf=-iotaf
ELSE
  WRITE(UNIT_stdOut,'(A)')'  ... VMEC signgs>0, so zeta and theta direction are switched.'
  !cos(-m*theta - (-n*zeta))=cos(m*theta-n*zeta)  !same sign for coefficients
  !only sin(-m*theta - (-n*zeta))=-sin(m*theta-n*zeta)  ! switch sign of coefficients
  zmns=-zmns
  lmns=-lmns
  IF(lasym) THEN
      rmns=-rmns
  END IF
END IF !not switchtheta

!WRITE VMEC_TO_GVEC BOUNDARY AND AXIS file
bc_Unit=GETFREEUNIT()
OPEN(UNIT     = bc_Unit       ,&
     FILE     = "vmec_to_gvec_boundary_and_axis.txt" ,&
     STATUS   = 'REPLACE'   ,&
     ACCESS   = 'SEQUENTIAL' )
  WRITE(bc_unit,'(A)')'! vmec boundary and axis written in the gvec coordinates (right-handed in rho,theta,zeta)'
  WRITE(bc_unit,'(A)')'!'
  WRITE(bc_unit,'(A)')'! non-zero boundary values of rmnc vmec -> X1_b_cos'
  DO iMode=1,mn_mode
    IF(ABS(rmnc(iMode,nFluxVMEC)).GT.1e-12)THEN
      WRITE(bc_unit,'("X1_b_cos(",I5,";",I5,")",X,"=",X,E22.15)')NINT(xm(iMode)),NINT(xn(iMode)/nfp),rmnc(iMode,nFluxVMEC)
    END IF
  END DO !imode
  IF(lasym)THEN
    WRITE(bc_unit,'(A)')'! non-zero boundary values of rmns vmec -> X1_b_sin'
    DO iMode=1,mn_mode
      IF(ABS(rmns(iMode,nFluxVMEC)).GT.1e-12)THEN
        WRITE(bc_unit,'("X1_b_sin(",I5,";",I5,")",X,"=",X,E22.15)')NINT(xm(iMode)),NINT(xn(iMode)/nfp),rmns(iMode,nFluxVMEC)
      END IF
    END DO !imode
  END IF
  IF(lasym)THEN
    WRITE(bc_unit,'(A)')'! non-zero boundary values of zmnc vmec -> X2_b_cos'
    DO iMode=1,mn_mode
      IF(ABS(zmnc(iMode,nFluxVMEC)).GT.1e-12)THEN
        WRITE(bc_unit,'("X2_b_cos(",I5,";",I5,")",X,"=",X,E22.15)')NINT(xm(iMode)),NINT(xn(iMode)/nfp),zmnc(iMode,nFluxVMEC)
      END IF
    END DO !imode
  END IF
  WRITE(bc_unit,'(A)')'! non-zero boundary values of zmns vmec -> X2_b_sin'
  DO iMode=1,mn_mode
    IF(ABS(zmns(iMode,nFluxVMEC)).GT.1e-12)THEN
      WRITE(bc_unit,'("X2_b_sin(",I5,";",I5,")",X,"=",X,E22.15)') NINT(xm(iMode)),NINT(xn(iMode)/nfp),zmns(iMode,nFluxVMEC)
    END IF
  END DO !imode
  !axis
  WRITE(bc_unit,'(A)')'! non-zero axis values of rmnc vmec -> X1_a_cos'
  DO iMode=1,mn_mode
    IF((NINT(xm(imode)).EQ.0).AND.(ABS(rmnc(iMode,nFluxVMEC)).GT.1e-12))THEN
      WRITE(bc_unit,'("X1_a_cos(",I1,";",I5,")",X,"=",X,E22.15)')0,NINT(xn(iMode)/nfp),rmnc(iMode,1)
    END IF
  END DO !imode
  IF(lasym)THEN
    WRITE(bc_unit,'(A)')'! non-zero axis values of rmns vmec -> X1_a_sin'
    DO iMode=1,mn_mode
      IF((NINT(xm(imode)).EQ.0).AND.(ABS(rmns(iMode,nFluxVMEC)).GT.1e-12))THEN
        WRITE(bc_unit,'("X1_a_sin(",I1,";",I5,")",X,"=",X,E22.15)')0,NINT(xn(iMode)/nfp),rmns(iMode,1)
      END IF
    END DO !imode
  END IF
  IF(lasym)THEN
    WRITE(bc_unit,'(A)')'! non-zero axis values of zmnc vmec -> X2_a_cos'
    DO iMode=1,mn_mode
      IF((NINT(xm(imode)).EQ.0).AND.(ABS(zmnc(iMode,nFluxVMEC)).GT.1e-12))THEN
        WRITE(bc_unit,'("X2_a_cos(",I1,";",I5,")",X,"=",X,E22.15)')0,NINT(xn(iMode)/nfp),zmnc(iMode,1)
      END IF
    END DO !imode
  END IF
  WRITE(bc_unit,'(A)')'! non-zero axis values of zmns vmec -> X2_a_sin'
  DO iMode=1,mn_mode
    IF((NINT(xm(imode)).EQ.0).AND.(ABS(zmns(iMode,nFluxVMEC)).GT.1e-12))THEN
      WRITE(bc_unit,'("X2_a_sin(",I1,";",I5,")",X,"=",X,E22.15)')0,NINT(xn(iMode)/nfp),zmns(iMode,1)
    END IF
  END DO !imode
CLOSE(bc_unit)


!toroidal flux from VMEC, now called PHI!
ALLOCATE(Phi_prof(nFluxVMEC))
Phi_prof = phi

!normalized toroidal flux (=flux variable s [0;1] in VMEC)
ALLOCATE(NormFlux_prof(nFluxVMEC))
NormFlux_prof=(Phi_prof-Phi_prof(1))/(Phi_prof(nFluxVMEC)-Phi_prof(1))
WRITE(UNIT_stdOut,'(4X,A,3F10.4)')'normalized toroidal flux of first three flux surfaces',NormFlux_prof(2:4)
!poloidal flux from VMEC
ALLOCATE(chi_prof(nFluxVMEC))
chi_prof=chi
WRITE(UNIT_stdOut,'(4X,A,3F10.4)')'min/max toroidal flux',MINVAL(phi)*TwoPi,MAXVAL(phi)*TwoPi
WRITE(UNIT_stdOut,'(4X,A,3F10.4)')'min/max poloidal flux',MINVAL(chi)*TwoPi,MAXVAL(chi)*TwoPi

WRITE(UNIT_stdOut,'(4X,A, I6)')'Total Number of flux surfaces: ',nFluxVMEC
WRITE(UNIT_stdOut,'(4X,A, I6)')'Total Number of mn-modes     : ',mn_mode
WRITE(UNIT_stdOut,'(4X,A,3I6)')'Max Mode m,n,nfp             : ',NINT(MAXVAL(xm)),NINT(MAXVAL(xn)),nfp
IF(lasym)THEN
  WRITE(UNIT_stdOut,'(6X,A,I6)')   '  lasym=T... ASYMMETRIC'
ELSE
  WRITE(UNIT_stdOut,'(6X,A,I6)')   '  lasym=F... SYMMETRIC (half fourier modes)'
END IF

useFilter=.TRUE. !GETLOGICAL('VMECuseFilter',Proposal=.TRUE.) !SHOULD BE ALWAYS TRUE...

ALLOCATE(xmabs(mn_mode))
DO iMode=1,mn_mode
  xmabs(iMode)=ABS(NINT(xm(iMode)))
  IF(useFilter)THEN
    IF(xmabs(iMode) > 3) THEN !Filtering for |m| > 3
      IF(MOD(xmabs(iMode),2) == 0) THEN
        xmabs(iMode)=2 !Even mode, remove rho**2
      ELSE
        xmabs(iMode)=3 !Odd mode, remove rho**3
      END IF
    END IF
  END IF !usefilter
END DO !iMode=1,mn_mode

!prepare Spline interpolation
ALLOCATE(rho(1:nFluxVMEC))
rho(:)=SQRT(NormFlux_prof(:))


ALLOCATE(Rmnc_Spl(mn_mode))
CALL FitSpline(mn_mode,nFluxVMEC,xmAbs,Rmnc,Rmnc_Spl)

ALLOCATE(Zmns_Spl(mn_mode))
CALL FitSpline(mn_mode,nFluxVMEC,xmAbs,Zmns,Zmns_Spl)

IF(lasym)THEN
  WRITE(Unit_stdOut,'(4X,A)')'LASYM=TRUE : R,Z,lambda in cos and sin!'
  ALLOCATE(Rmns_Spl(mn_mode))
  CALL FitSpline(mn_mode,nFluxVMEC,xmAbs,Rmns,Rmns_Spl)

  ALLOCATE(Zmnc_Spl(mn_mode))
  CALL FitSpline(mn_mode,nFluxVMEC,xmAbs,Zmnc,Zmnc_Spl)

END IF


ALLOCATE(lmns_Spl(mn_mode))
IF(lasym) ALLOCATE(lmnc_Spl(mn_mode))
!WRITE(*,*)'DEBUG:lambda_grid:',lambda_grid
IF(lambda_grid.EQ."half")THEN
  !lambda given on half grid
  CALL           FitSplineHalf(mn_mode,nFluxVMEC,xmAbs,lmns,lmns_Spl)
  IF(lasym) CALL FitSplineHalf(mn_mode,nFluxVMEC,xmAbs,lmnc,lmnc_Spl)
ELSEIF(lambda_grid.EQ."full")THEN
  CALL           FitSpline(mn_mode,nFluxVMEC,xmAbs,lmns,lmns_Spl)
  IF(lasym) CALL FitSpline(mn_mode,nFluxVMEC,xmAbs,lmnc,lmnc_Spl)
ELSE
  CALL abort(__STAMP__, &
             'no lambda_grid found!!!! lambda_grid='//TRIM(lambda_grid) )
END IF


CALL interpolate_cubic_spline(rho**2,presf,c_pres,knots_pres, (/0,0/))
vmec_pres_profile = t_rProfile_bspl(knots_pres, c_pres)
SDEALLOCATE(knots_pres)
SDEALLOCATE(c_pres)

CALL interpolate_cubic_spline(rho**2,Phi_prof,c_phi,knots_phi,(/0,0/))
vmec_phi_profile = t_rProfile_bspl(knots_phi, c_phi)
SDEALLOCATE(knots_phi)
SDEALLOCATE(c_phi)

CALL interpolate_cubic_spline(rho**2,chi_prof,c_chi,knots_chi,(/0,0/))
vmec_chi_profile = t_rProfile_bspl(knots_chi, c_chi)
SDEALLOCATE(knots_chi)
SDEALLOCATE(c_chi)

CALL interpolate_cubic_spline(rho**2,iotaf,c_iota,knots_iota,(/0,0/))
vmec_iota_profile = t_rProfile_bspl(knots_iota, c_iota)
SDEALLOCATE(knots_iota)
SDEALLOCATE(c_iota)

WRITE(Unit_stdOut,'(4X,A,3F12.4)')'tor. flux Phi  axis/middle/edge',Phi(1)*TwoPi,Phi(nFluxVMEC/2)*TwoPi,Phi(nFluxVMEC)*TwoPi
WRITE(Unit_stdOut,'(4X,A,3F12.4)')'pol. flux chi  axis/middle/edge',chi(1)*TwoPi,chi(nFluxVMEC/2)*TwoPi,chi(nFluxVMEC)*TwoPi
WRITE(Unit_stdOut,'(4X,A,3F12.4)')'   iota        axis/middle/edge',vmec_iota_profile%eval_at_rho(0.0_wp),&
vmec_iota_profile%eval_at_rho(rho(nFluxVMEC/2)),vmec_iota_profile%eval_at_rho(1.0_wp)
WRITE(Unit_stdOut,'(4X,A,3F12.4)')'  pressure     axis/middle/edge',vmec_pres_profile%eval_at_rho(0.0_wp),&
vmec_pres_profile%eval_at_rho(rho(nFluxVMEC/2)),vmec_pres_profile%eval_at_rho(1.0_wp)


WRITE(UNIT_stdOut,'(A)')'  ... INIT VMEC DONE.'
WRITE(UNIT_stdOut,fmt_sep)

END SUBROUTINE InitVMEC