fixaray.f90

Go to the documentation of this file.

 
 
SUBROUTINE fixaray
 
  USE vmec_main, p5 => cp5
  USE vmec_params, ONLY: jmin2, mscale, nscale, mnyq, nnyq, signgs
 
  use dbgout
 
  IMPLICIT NONE
 
  REAL(rprec), PARAMETER :: two=2.0_dp
 
  REAL(rprec), PARAMETER :: pexp=4.0_dp ! for <M> spectral width screen diagnostic
 
  INTEGER :: i, m, j, n, mn, mn1, nmin0, istat1, istat2
  INTEGER :: mnyq0, nnyq0
  REAL(rprec):: argi, arg, argj, dnorm, dnorm3
  logical :: dbg_fixaray
  real(rprec), allocatable :: arg_mu(:,:), arg_nv(:,:)
 
 ! COMPUTE TRIGONOMETRIC FUNCTION ARRAYS
 ! NOTE: ARRAYS ALLOCATED HERE ARE GLOBAL AND ARE DEALLOCATED IN FILEOUT
 ! NOTE: NEED 2 X NYQUIST FOR FAST HESSIAN CALCULATIONS
  mnyq0  = ntheta1/2 ! maximum mode numbers supported by grid
  nnyq0  = nzeta/2
 
  ! make sure that mnyq, nnyq are at least twice mpol-1, ntor
  ! or large enough to fully represent the information held in realspace (mnyq0, nnyq0)
  mnyq = max(0, 2*mnyq0, 2*mpol1)
  nnyq = max(0, 2*nnyq0, 2*ntor)
 
  mnmax_nyq = nnyq/2 + 1 + mnyq/2 * (2 * nnyq/2 + 1)
 
  ALLOCATE(cosmu(ntheta3,0:mnyq),  sinmu(ntheta3,0:mnyq),           &
           cosmum(ntheta3,0:mnyq), sinmum(ntheta3,0:mnyq),          &
           cosmui(ntheta3,0:mnyq), cosmumi(ntheta3,0:mnyq),         &
           cosmui3(ntheta3,0:mnyq),cosmumi3(ntheta3,0:mnyq),        &
           sinmui(ntheta3,0:mnyq), sinmumi(ntheta3,0:mnyq),         &
           cosnv(nzeta,0:nnyq),    sinnv(nzeta,0:nnyq),             &
           cosnvn(nzeta,0:nnyq),   sinnvn(nzeta,0:nnyq),            &
           cos01(nznt),            sin01(nznt),            stat=istat1 )
  ALLOCATE(xm(mnmax), xn(mnmax), ixm(mnsize), jmin3(0:mnsize-1),    &
           xm_nyq(mnmax_nyq), xn_nyq(mnmax_nyq),                    &
           mscale(0:mnyq), nscale(0:nnyq), stat=istat2)
  allocate(arg_mu(ntheta3,0:mnyq), arg_nv(nzeta,0:nnyq))
 
  IF (istat1.ne.0) stop 'allocation error in fixaray: istat1'
  IF (istat2.ne.0) stop 'allocation error in fixaray: istat2'
 
  ! Normalization factor for forward Fourier transforms (e.g. tomnsp*)
  ! In any case (symmetric or asymmetric), the Fourier integrals in tomnsp* and alias
  ! (the only place where dnorm is used through cosmui, sinmui, etc.)
  ! are ever only taken over [0, pi] corresponding to 1, ..., ntheta2.
  dnorm = one/(nzeta*(ntheta2-1.0_dp))
 
  ! Normalization factor for surface integrals/averages (wint is based on cosmui3).
  ! For the asymmetric case, the norm in wint thus has to be 1/(nzeta*ntheta3).
  IF (lasym) then
    dnorm3 = one/(nzeta*ntheta1)
  else
    dnorm3 = one/(nzeta*(ntheta2-1.0_dp))
  end if
 
  ! (from vmec_params)
  ! array for norming theta-trig functions (internal use only)
  ! so that the discrete SUM[cos(mu)*cos(m'u)] = .5 delta(m,m')
  ! and analogously for zeta/v
  mscale(0) = 1.0_dp
  nscale(0) = 1.0_dp
  mscale(1:mnyq) = mscale(0)/osqrt2 ! == sqrt(2)
  nscale(1:nnyq) = nscale(0)/osqrt2 ! == sqrt(2)
 
  r0scale = mscale(0)*nscale(0) ! == 1
 
  ! GENERALLY, ONLY NEED THIS FROM 1, ntheta2 EXCEPT IN GETBRHO ROUTINE
  DO i = 1, ntheta3
     argi = twopi*(i-1)/ntheta1
     DO m = 0, mnyq
        arg = argi*m
 
        arg_mu(i,m) = arg
 
        cosmu(i,m) = cos(arg)*mscale(m)
        sinmu(i,m) = sin(arg)*mscale(m)
 
        cosmui(i,m) = dnorm*cosmu(i,m)
        sinmui(i,m) = dnorm*sinmu(i,m)
        IF (i.EQ.1 .OR. i.EQ.ntheta2) then
           ! Trapezoidal integration requires a factor of 1/2 for the first and the last point.
           ! Note that this is done also in the case of an asymmetric run!
           ! There, cosmui is only used up to ntheta2 anyway...
           cosmui(i,m)=cosmui(i,m)/2.0_dp
        end if
 
        ! Use this if integration over FULL 1,ntheta3 interval
        cosmui3(i,m) = dnorm3*cosmu(i,m)
        IF (.not.lasym .and. (i.eq.1 .or. i.eq.ntheta2)) then
           ! Note that cosmui3 is only ever used to construct the wint array.
           ! where only the m=0 component of cosmui3 enters.
           cosmui3(i,m) = cosmui3(i,m)/2.0_dp
        end if
 
        cosmum(i,m)   = cosmu(i,m)*(m)
        sinmum(i,m)   =-sinmu(i,m)*(m)
        cosmumi(i,m)  = cosmui(i,m)*(m)
        sinmumi(i,m)  =-sinmui(i,m)*(m)
 
        cosmumi3(i,m) = cosmui3(i,m)*m ! not used... --> remove?
     END DO
  END DO
 
  DO j = 1, nzeta
     argj = twopi*(j-1)/nzeta
     DO n = 0, nnyq
        arg = argj*(n)
 
        arg_nv(j,n) = arg
 
        cosnv(j,n) = cos(arg)*nscale(n)
        sinnv(j,n) = sin(arg)*nscale(n)
        cosnvn(j,n) =  cosnv(j,n)*(n*nfp)
        sinnvn(j,n) = -sinnv(j,n)*(n*nfp)
     END DO
  END DO
 
  ! R,Z,L / s**(m/2) ARE LINEAR NEAR ORIGIN
  mn = 0
  mn1 = 0
  DO m = 0, mpol1
     xmpq(m,1) = m*(m - 1.0_dp)   ! used for spectral constraint force --> m^2-m
 
     ! xmpq(m,2:3) are ONLY used for screen diagnostic <M> !!!
     xmpq(m,2) = m**pexp     ! m^p     with p = pexp = 4
     xmpq(m,3) = m**(pexp+1.0_dp) ! m^(p+q) with q = 1
 
     ! compute ixm == _i_nteger version of xm
     DO n = 0, ntor
        jmin3(mn) = jmin2(m) ! vmec_params: starting js(m) values for which R,Z are evolved
        mn = mn + 1
        ixm(mn) = m
     END DO
 
     ! use this loop also to compute xm, xn array contents
     nmin0 = -ntor
     IF (m .eq. 0) nmin0 = 0
     DO n = nmin0, ntor
        mn1 = mn1 + 1
        xm(mn1) = m
        xn(mn1) = n*nfp
     END DO
  END DO
 
  IF (mn1 .ne. mnmax) stop 'mn1 != mnmax'
 
  dbg_fixaray = open_dbg_context("fixaray")
  if (dbg_fixaray) then
 
    call add_real_2d("arg_mu",    ntheta3, mnyq+1, arg_mu(1:ntheta3, 0:mnyq))
    call add_real_2d("arg_nv",    nzeta,   nnyq+1, arg_nv(1:nzeta,   0:nnyq))
 
    call add_real_2d("cosmu",    ntheta2, mnyq+1, cosmu(1:ntheta2,:))
    call add_real_2d("sinmu",    ntheta2, mnyq+1, sinmu(1:ntheta2,:))
    call add_real_2d("cosmum",   ntheta2, mnyq+1, cosmum(1:ntheta2,:))
    call add_real_2d("sinmum",   ntheta2, mnyq+1, sinmum(1:ntheta2,:))
    call add_real_2d("cosmui",   ntheta2, mnyq+1, cosmui(1:ntheta2,:))
    call add_real_2d("sinmui",   ntheta2, mnyq+1, sinmui(1:ntheta2,:))
    call add_real_2d("cosmui3",  ntheta2, mnyq+1, cosmui3(1:ntheta2,:))
    call add_real_2d("cosmumi",  ntheta2, mnyq+1, cosmumi(1:ntheta2,:))
    call add_real_2d("sinmumi",  ntheta2, mnyq+1, sinmumi(1:ntheta2,:))
    call add_real_2d("cosmumi3", ntheta2, mnyq+1, cosmumi3(1:ntheta2,:))
 
    call add_real_2d("cosnv",  nzeta, nnyq+1, cosnv)
    call add_real_2d("sinnv",  nzeta, nnyq+1, sinnv)
    call add_real_2d("cosnvn", nzeta, nnyq+1, cosnvn)
    call add_real_2d("sinnvn", nzeta, nnyq+1, sinnvn)
 
    call add_real_1d("mscale", mnyq+1, mscale)
    call add_real_1d("nscale", nnyq+1, nscale)
 
    ! TODO: add r0scale output (currently in residue...)
 
  end if
 
  ! COMPUTE NYQUIST-SIZED ARRAYS FOR OUTPUT.
  ! RESTORE m,n Nyquist TO 1 X ... (USED IN WROUT, JXBFORCE)
  !  mnyq = mnyq0;  nnyq = nnyq0
  mnyq = mnyq/2
  nnyq = nnyq/2
 
  mn1 = 0
  DO m = 0, mnyq
     nmin0 = -nnyq
     IF (m .eq. 0) nmin0 = 0
     DO n = nmin0, nnyq
        mn1 = mn1 + 1
        xm_nyq(mn1) = m
        xn_nyq(mn1) = n*nfp
     END DO
  END DO
 
  IF (mn1 .ne. mnmax_nyq) stop 'mn1 != mnmax_nyq'
 
  ! cos01 and sin01 are used as trigmult for precondn in bcovar
  ! cos01:  m*cos(m u - n v) for m=1, n=0 --> d(sin(mu))/du
  ! sin01: -m*sin(m u - n v) for m=1, n=0 --> d(cos(mu))/du
  ! evaluated on (ntheta3 x nzeta) grid
  mn = 0 ! mn is re-used as linear grid index here
  m = 1
  DO i = 1, ntheta3
     argi = twopi*(i - 1)/ntheta1
     DO j = 1, nzeta
        mn = mn + 1
        cos01(mn) = m*cos(m*argi)*mscale(m)
        sin01(mn) =-m*sin(m*argi)*mscale(m)
     END DO
  END DO
 
  ! _fac_tor for _con_straint
  faccon(0) = zero
  faccon(1:mpol1-1) = -0.25_dp*signgs/xmpq(2:mpol1,1)**2.0_dp
  faccon(mpol1) = zero
 
  if (dbg_fixaray) then
 
    call add_int("ntheta3", ntheta3)
    call add_int("mnyq", mnyq)
    call add_int("nzeta", nzeta)
    call add_int("nnyq", nnyq)
    call add_int("nznt", nznt)
    call add_int("mnmax", mnmax)
    call add_int("mnsize", mnsize)
    call add_int("mnmax_nyq", mnmax_nyq)
 
    call add_real_1d("cos01", nznt, cos01)
    call add_real_1d("sin01", nznt, sin01)
 
    call add_real_1d("xm", mnmax, xm)
    call add_real_1d("xn", mnmax, xn)
 
    call add_real_1d("xm_nyq", mnmax_nyq, xm_nyq)
    call add_real_1d("xn_nyq", mnmax_nyq, xn_nyq)
 
    call add_int_1d("ixm", mnsize, ixm)
 
    call close_dbg_out()
  end if
 
  if (open_dbg_context("spectral_constraint", num_eqsolve_retries)) then
 
    ! xmpq is allocated statically, so need size here explicitly!
    call add_real_2d("xmpq", mpol, 3, xmpq(0:mpol1,1:3))
 
    call add_real_1d("faccon", mpol, faccon)
 
    call close_dbg_out()
  end if
 
END SUBROUTINE fixaray