Change psi=-psi_t to psi=sgn(psi_t)*psi_t (#126)

* Change psi=-psi_t to psi=sgn(psi_t)*psi_t

geo/stella_geometry.f90

@@ -27,6 +27,7 @@ module stella_geometry
    public :: zeta
    public :: zed_scalefac
    public :: dxdXcoord, dydalpha
+   public :: sign_torflux
    public :: aref, bref
    public :: twist_and_shift_geo_fac
    public :: q_as_x, get_x_to_rho, gfac
@@ -71,6 +72,7 @@ module stella_geometry
    real, dimension(:), allocatable :: alpha
    real, dimension(:, :), allocatable :: zeta
 
+   integer :: sign_torflux
    integer :: geo_option_switch
    integer, parameter :: geo_option_local = 1
    integer, parameter :: geo_option_inputprof = 2
@@ -120,7 +122,6 @@ subroutine init_geometry(nalpha, naky)
       real :: dpsidrho, dpsidrho_psi0
       real :: new_zeta_min
       integer :: iy, ia, iz
-      integer :: sign_torflux
       integer :: dxdXcoord_sign, dydalpha_sign
       real :: field_period_ratio, bmag_z0
       real, dimension(:, :), allocatable :: grad_alpha_grad_alpha
@@ -171,6 +172,7 @@ subroutine init_geometry(nalpha, naky)
             drhodpsi_psi0 = 1./dpsidrho_psi0
             ! dxdXcoord = a*Bref*dx/dpsi = sign(dx/dpsi) * a*q/r
             dxdXcoord_sign = 1
+            sign_torflux = -1
             if (q_as_x) then
                dxdXcoord = dxdXcoord_sign * dpsidrho
             else
@@ -224,6 +226,7 @@ subroutine init_geometry(nalpha, naky)
             drhodpsi_psi0 = 1./dpsidrho_psi0
             ! dxdXcoord = a*Bref*dx/dpsi = sign(dx/dpsi) * a*q/r
             dxdXcoord_sign = 1
+            sign_torflux = -1
             if (q_as_x) then
                dxdXcoord = dxdXcoord_sign / drhodpsi_psi0
             else
@@ -278,6 +281,7 @@ subroutine init_geometry(nalpha, naky)
             ! dxdXcoord = a*Bref*dx/dpsi = sign(dx/dpsi) * a*q/r
             dxdXcoord_sign = 1
             dxdXcoord = dxdXcoord_sign * geo_surf%qinp / geo_surf%rhoc
+            sign_torflux = -1
             ! dydalpha = (dy/dalpha) / a = sign(dydalpha) * (dpsi/dr) / (a*Bref)
             dydalpha_sign = 1
             dydalpha = dydalpha_sign * dpsidrho
@@ -348,20 +352,21 @@ subroutine init_geometry(nalpha, naky)
                ! Restart the variable twist_and_shift_geo_fac_full
                twist_and_shift_geo_fac_full = 0
             end if
-            !> Bref = 2*abs(psi_tor_LCFS)/a^2
-            !> a*Bref*dx/dpsi_tor = sign(psi_tor)/rhotor
-            !> psi = -psi_tor
-            !> dxdXcoord = a*Bref*dx/dpsi = -a*Bref*dx/dpsi_tor = -sign(psi_tor)/rhotor
-            dxdXcoord_sign = -1
-            dxdXcoord = dxdXcoord_sign * sign_torflux / geo_surf%rhotor
-            !> dydalpha = (dy/dalpha) / a = sign(dydalpha) * rhotor
-            dydalpha_sign = 1
-            dydalpha = dydalpha_sign * geo_surf%rhotor
-            !> if using vmec, rho = sqrt(psitor/psitor_lcfs)
-            !> psiN = -psitor/(aref**2*Bref)
-            !> so drho/dpsiN = -drho/d(rho**2) * (aref**2*Bref/psitor_lcfs) = -1.0/rho
-            drhodpsi = dxdXcoord_sign * sign_torflux / geo_surf%rhotor
+
+            !> Define <dxdXcoord> = (ρref/a)(dx̃/dψ̃) and use dx/dψ = 1/(a*ρ0*Bref)
+            !> dxdXcoord = (ρref/a)(dx̃/dψ̃) = (ρref/a)(d(x/ρref)/d(ψ/(a^2 Bref)) = a Bref (dx/dψ) = 1/ρ0
+            dxdXcoord = 1 / geo_surf%rhotor
+
+            !> Define <dydalpha> = (ρref/a)(dỹ/dα̃) and use dy/dα = a*ρ0
+            !> dydalpha = (ρref/a)(dỹ/dα̃) = (ρref/a)(d(y/ρref)/dα) = (1/a) (dy/dα) = ρ0
+            dydalpha = geo_surf%rhotor
+
+            !> Define <drhodpsi> = dρ0/dψ̃ and use ρ0 = sqrt(psi_t/psi_{t,LCFS}) and ψ = sgn(psi_t)*psi_t
+            !> Use dρ0/dpsi_t = d(sqrt(psi_t/psi_{t,LCFS}))/dpsi_t = sgn(psi_t)/(a^2 ρ0 Bref) with Bref = 2 |psi_{t,LCFS}| / a^2
+            !> drhodpsi = dρ0/dψ̃ = dρ0/d(ψ/(a^2 Bref)) = a^2 Bref sgn(psi_t) dρ0/dpsi_t = 1/ρ0
+            drhodpsi = 1 / geo_surf%rhotor
             drhodpsi_psi0 = drhodpsi
+
             bmag_psi0 = bmag
             grad_x_grad_y_end = grad_alpha_grad_psi(1, nzgrid) * (aref * aref * bref)
             select case (boundary_option_switch)
@@ -393,9 +398,12 @@ subroutine init_geometry(nalpha, naky)
             !> gds2 = |grad y|^2 = |grad alpha|^2 * (dy/dalpha)^2
             !> note that rhotor = sqrt(psi/psi_LCFS)
             gds2 = grad_alpha_grad_alpha * dydalpha**2
-            !> gds21 = shat * grad x . grad y = shat * dx/dpsi_t * dy/dalpha * grad alpha . grad psi_t
-            !> NB: psi = -psi_t and so dx/dpsi = = dx/dpsi_t, which is why there is a minus sign here
-            gds21 = -grad_alpha_grad_psi * geo_surf%shat * dxdXcoord * dydalpha
+
+            !> Define <gds21> = hat{s} ∇x . ∇y
+            !> Use (dx/dψ)*(dy/dα) = 1/(a ρ0 Bref) * (a ρ0) = 1/Bref
+            !> Use ∇x . ∇y = (dx/dψ)(dy/dα) ∇ψ . ∇α = (1/Bref) ∇ψ . ∇α = <grad_alpha_grad_psi>
+            gds21 = grad_alpha_grad_psi * geo_surf%shat
+
             !> gds22 = shat^2 * |grad x|^2 = shat^2 * |grad psi_t|^2 * (dx/dpsi_t)^2
             gds22 = (geo_surf%shat * grad_x)**2
             !gds22 = geo_surf%shat**2 * grad_psi_grad_psi * dxdXcoord**2
@@ -426,11 +434,11 @@ subroutine init_geometry(nalpha, naky)
          ! exb_nonlin_fac is equivalent to kxfac/2 in gs2
          if (q_as_x) then
             dqdrho = geo_surf%shat * geo_surf%qinp / geo_surf%rhoc
-            exb_nonlin_fac = 0.5 * dxdXcoord * dydalpha * drhodpsi * dqdrho
+            exb_nonlin_fac = -0.5 * sign_torflux * dxdXcoord * dydalpha * drhodpsi * dqdrho
             !the following will get multiplied by exb_nonlin_fac in advance_exb_nonlinearity
             exb_nonlin_fac_p = geo_surf%d2qdr2 / dqdrho - geo_surf%d2psidr2 * drhodpsi
          else
-            exb_nonlin_fac = 0.5 * dxdXcoord * dydalpha
+            exb_nonlin_fac = -0.5 * sign_torflux * dxdXcoord * dydalpha
             exb_nonlin_fac_p = 0.0
          end if
       end if
@@ -458,7 +466,7 @@ subroutine init_geometry(nalpha, naky)
       !> jacob is the Jacobian from Cartesian coordinates to (y,x,z) coordinates
       !> is ((grad y x grad x) . grad z)^(-1) = Lref*(dalpha/dy)*(dpsi/dx)/(Lref*Bref)*(B/Bref . grad z)^(-1)
       !> Lref*(dalpha/dy) = 1/dydalpha; (dpsi/dx)/(Lref*Bref) = 1 / dxdXcoord ; (B/Bref . grad z) = gradpar*bmag
-      jacob = 1.0 / (dydalpha * dxdXcoord * b_dot_grad_z * bmag)
+      jacob = -sign_torflux / (dydalpha * dxdXcoord * b_dot_grad_z * bmag)
       !    jacob = 1.0/(dydalpha*dxdXcoord*spread(gradpar,1,nalpha)*bmag)
 
       ! this is dl/B
@@ -846,6 +854,7 @@ subroutine broadcast_arrays
       call broadcast(zeta)
       call broadcast(dxdXcoord)
       call broadcast(dydalpha)
+      call broadcast(sign_torflux)
       call broadcast(twist_and_shift_geo_fac)
       call broadcast(grad_x_grad_y_end)
       call broadcast(vmec_chosen)

geo/vmec_interface/vmec_to_stella_geometry_interface.f90

@@ -1329,9 +1329,12 @@ subroutine vmec_to_stella_geometry_interface(nalpha, alpha0, nzgrid, &
       ! Compute the Cartesian components of other quantities we need:
       !*********************************************************************
 
-      grad_psi_X = grad_s_X * edge_toroidal_flux_over_2pi
-      grad_psi_Y = grad_s_Y * edge_toroidal_flux_over_2pi
-      grad_psi_Z = grad_s_Z * edge_toroidal_flux_over_2pi
+      ! In VMEC <s> = psi/psi_LCFS and in stella <psi> = sgn(psi_t)*psi_t
+      ! nabla psi = (dpsi/ds) * nabla s = d(psi_LCFS*s)/ds * nabla s
+      !           = psi_LCFS * nabla s = sgn(psi_t) * psi_{t,LCFS} * nabla s
+      grad_psi_X = grad_s_X * sign_toroidal_flux * edge_toroidal_flux_over_2pi
+      grad_psi_Y = grad_s_Y * sign_toroidal_flux * edge_toroidal_flux_over_2pi
+      grad_psi_Z = grad_s_Z * sign_toroidal_flux * edge_toroidal_flux_over_2pi
 
       ! Form grad alpha = grad (theta_vmec + Lambda - iota * zeta)
       do izeta = -nzgrid, nzgrid
@@ -1516,10 +1519,10 @@ subroutine vmec_to_stella_geometry_interface(nalpha, alpha0, nzgrid, &
 !    gds22 = (grad_psi_X * grad_psi_X + grad_psi_Y * grad_psi_Y + grad_psi_Z * grad_psi_Z) &
 !         * shat * shat / (L_reference * L_reference * B_reference * B_reference * normalized_toroidal_flux_used)
 
-      ! this is (grad zeta . grad x_stella) / bmag^2 = (grad zeta . grad psitor) * dx/dpsitor / bmag^2
-      ! = (grad zeta . grad psitor) * sign_torflux/rhotor/Lref/Bref/ bmag^2
-      gradzeta_gradx = sign_toroidal_flux &
-                       * (grad_zeta_X * grad_psi_X + grad_zeta_Y * grad_psi_Y + grad_zeta_Z * grad_psi_Z) &
+      ! Define <gradzeta_gradx> = (grad zeta . grad x) / (B/Bref)^2
+      ! gradzeta_gradx = (grad zeta . grad x) / bmag^2 = (grad zeta . grad psi) * dx/dpsi / bmag^2
+      !                = (grad zeta . grad psi) / (rhotor*Lref*Bref) / bmag^2
+      gradzeta_gradx = (grad_zeta_X * grad_psi_X + grad_zeta_Y * grad_psi_Y + grad_zeta_Z * grad_psi_Z) &
                        / (L_reference * B_reference * sqrt(normalized_toroidal_flux_used) * bmag**2)
 
       ! this is (grad zeta . grad y_stella) / bmag^2 = (grad zeta . grad alpha) * dy/dalpha / bmag^2
@@ -1537,37 +1540,51 @@ subroutine vmec_to_stella_geometry_interface(nalpha, alpha0, nzgrid, &
                          + grad_theta_pest_Z * grad_alpha_Z) &
                         * L_reference * sqrt(normalized_toroidal_flux_used) / bmag**2
 
-      ! this is psitor/|psitor|*((grad y_stella . grad zeta)*(grad x_stella . grad y_stella)
-      ! - (grad x_stella . grad zeta)*|grad y_stella|^2) / (B/Bref)^2
-      gds23 = sign_toroidal_flux * (gradzeta_grady * sign_toroidal_flux * grad_alpha_grad_psi &
-                                    - gradzeta_gradx * grad_alpha_grad_alpha * normalized_toroidal_flux_used)
-!    gds23 = sign_toroidal_flux*(gradzeta_grady * gds21/shat - gradzeta_gradx * gds2)
-
-      ! this is psitor/|psitor| * ((grad y_stella . grad zeta) * |grad x_stella|^2
-      ! - (grad x_stella . grad zeta)*(grad x_stella . grad y_stella)) / (B/Bref)^2
-      gds24 = (gradzeta_grady * grad_psi_grad_psi / normalized_toroidal_flux_used &
-               - gradzeta_gradx * sign_toroidal_flux * grad_alpha_grad_psi) * sign_toroidal_flux
-!    gds24 = (gradzeta_grady * gds22/shat**2 - gradzeta_gradx * gds21/shat) * sign_toroidal_flux
-
-      ! this is ((grad y_stella . grad theta_pest)*(grad x_stella . grad y_stella)
-      ! - (grad x_stella . grad theta_pest)*|grad y_stella|^2) / (B/Bref)^2
-      gds25 = gradtheta_grady * sign_toroidal_flux * grad_alpha_grad_psi &
-              - gradtheta_gradx * grad_alpha_grad_alpha * normalized_toroidal_flux_used
-!    gds25 = gradtheta_grady * gds21/shat - gradtheta_gradx * gds2
-
-      ! this is ((grad y_stella . grad theta_pest) * |grad x_stella|^2
-      ! - (grad x_stella . grad theta_pest)*(grad x_stella . grad y_stella)) / 2 / (B/Bref)^2
-      gds26 = 0.5 * (gradtheta_grady * grad_psi_grad_psi / normalized_toroidal_flux_used &
-                     - gradtheta_gradx * sign_toroidal_flux * grad_alpha_grad_psi)
-!    gds26 = 0.5*(gradtheta_grady * gds22/shat**2 - gradtheta_gradx * gds21/shat)
+      ! Define <gds23> = sgn(psi_t)/(B/Bref)^2 * [(∇y . ∇ζ) * (∇x . ∇y) - (∇x . ∇ζ) * |∇y|^2]
+      ! Use ∇x . ∇y = (dy/dα)(dx/dψ) ∇α . ∇ψ = (1/Bref) ∇α . ∇ψ = <grad_alpha_grad_psi>
+      ! Use |∇y|^2 = (dy/dα)^2 |∇α|^2 = a^2 ρ^2 |∇α|^2 = a^2 |∇α|^2 * ρ^2 = <grad_alpha_grad_alpha> * <rhotor>^2
+      ! Use (∇y . ∇ζ)/(B/Bref)^2 = <gradzeta_grady> and (∇x . ∇ζ)/(B/Bref)^2 = <gradzeta_gradx>
+      ! Use ρ = sqrt(psi_t/psi_{t,LCFS}) thus ρ^2 = s = psi_t/psi_{t,LCFS} = <normalized_toroidal_flux_used>
+      ! <gds23> = sgn(psi_t)*(gradzeta_grady*grad_alpha_grad_psi - gradzeta_gradx*grad_alpha_grad_alpha*rhotor^2)
+      gds23 = sign_toroidal_flux * (gradzeta_grady * grad_alpha_grad_psi - gradzeta_gradx * grad_alpha_grad_alpha * normalized_toroidal_flux_used)
+
+      ! Define <gds24> = sgn(psi_t)/(B/Bref)^2 * [(∇y . ∇ζ) * |∇x|^2 - (∇x . ∇ζ) * (∇x . ∇y)]
+      ! Use ∇x . ∇y = (dy/dα)(dx/dψ) ∇α . ∇ψ = (1/Bref) ∇α . ∇ψ = <grad_alpha_grad_psi>
+      ! Use |∇x|^2 = (dx/dψ)^2 |∇ψ|^2 = 1/(a^2 ρ^2 Bref^2) |∇ψ|^2 = 1/(a^2 Bref^2) |∇ψ|^2 * 1/ρ^2  = <grad_psi_grad_psi> / <rhotor>^2
+      ! Use (∇y . ∇ζ)/(B/Bref)^2 = <gradzeta_grady> and (∇x . ∇ζ)/(B/Bref)^2 = <gradzeta_gradx>
+      ! Use ρ = sqrt(psi_t/psi_{t,LCFS}) thus ρ^2 = s = psi_t/psi_{t,LCFS} = <normalized_toroidal_flux_used>
+      ! <gds24> = sgn(psi_t)*(gradzeta_grady*grad_psi_grad_psi/rhotor^2 - gradzeta_gradx*grad_alpha_grad_psi)
+      gds24 = sign_toroidal_flux * (gradzeta_grady * grad_psi_grad_psi / normalized_toroidal_flux_used - gradzeta_gradx * grad_alpha_grad_psi)
+
+      ! Define <gds25> = sgn(psi_t)/(B/Bref)^2 * [(∇y . ∇θp) * |∇x|^2 - (∇x . ∇θp) * (∇x . ∇y)]
+      ! Use ∇x . ∇y = (dy/dα)(dx/dψ) ∇α . ∇ψ = (1/Bref) ∇α . ∇ψ = <grad_alpha_grad_psi>
+      ! Use |∇y|^2 = (dy/dα)^2 |∇α|^2 = a^2 ρ^2 |∇α|^2 = a^2 |∇α|^2 * ρ^2 = <grad_alpha_grad_alpha> * <rhotor>^2
+      ! Use (∇y . ∇θp)/(B/Bref)^2 = <gradtheta_grady> and (∇x . ∇θp)/(B/Bref)^2 = <gradtheta_gradx>
+      ! Use ρ = sqrt(psi_t/psi_{t,LCFS}) thus ρ^2 = s = psi_t/psi_{t,LCFS} = <normalized_toroidal_flux_used>
+      ! <gds25> = sgn(psi_t)*(gradtheta_grady*grad_alpha_grad_psi - gradtheta_gradx*grad_alpha_grad_alpha*rhotor^2)
+      ! Note that <gds25> was wrong before 17/09/2023 (branch upgrade_clockwise_vmecs) for CCW equilibria
+      gds25 = gradtheta_grady * grad_alpha_grad_psi - gradtheta_gradx * grad_alpha_grad_alpha * normalized_toroidal_flux_used
+
+      ! Define <gds26> = 1/2 * sgn(psi_t)/(B/Bref)^2 * [(∇y . ∇θp) * |∇x|^2 - (∇x . ∇θp) * (∇x . ∇y)]
+      ! Use ∇x . ∇y = (dy/dα)(dx/dψ) ∇α . ∇ψ = (1/Bref) ∇α . ∇ψ = <grad_alpha_grad_psi>
+      ! Use |∇x|^2 = (dx/dψ)^2 |∇ψ|^2 = 1/(a^2 ρ^2 Bref^2) |∇ψ|^2 = 1/(a^2 Bref^2) |∇ψ|^2 * 1/ρ^2  = <grad_psi_grad_psi> / <rhotor>^2
+      ! Use (∇y . ∇θp)/(B/Bref)^2 = <gradtheta_grady> and (∇x . ∇θp)/(B/Bref)^2 = <gradtheta_gradx>
+      ! Use ρ = sqrt(psi_t/psi_{t,LCFS}) thus ρ^2 = s = psi_t/psi_{t,LCFS} = <normalized_toroidal_flux_used>
+      ! <gds26> = 0.5*sgn(psi_t)*(gradtheta_grady*grad_psi_grad_psi/rhotor^2 - gradtheta_gradx*grad_alpha_grad_psi)
+      ! Note that <gds26> was wrong before 17/09/2023 (branch upgrade_clockwise_vmecs) for CCW equilibria
+      gds26 = 0.5 * (gradtheta_grady * grad_psi_grad_psi / normalized_toroidal_flux_used - gradtheta_gradx * grad_alpha_grad_psi)
 
       gbdrift_alpha = 2 * B_reference * L_reference * L_reference * B_cross_grad_B_dot_grad_alpha &
                       / (B * B * B)
 !    gbdrift = 2 * B_reference * L_reference * L_reference * sqrt_s * B_cross_grad_B_dot_grad_alpha &
 
-      gbdrift0_psi = -edge_toroidal_flux_over_2pi &
-                     * (B_sub_theta_vmec * d_B_d_zeta - B_sub_zeta * d_B_d_theta_vmec) / sqrt_g &
-                     * 2 * shat / (B * B * B)
+      ! Define <gbdrift0_psi> = hat{s} * (2/B^3) * B x ∇B · ∇ψ
+      ! Use hat{s} * (2/B^3) = 2 * shat / (B * B * B)
+      ! Use B x ∇B · ∇ψ = psi_LCFS / sqrt(g) * ( ∂B/∂ζ * Btheta - ∂B/∂θ * Bzeta )
+      ! Use psi_LCFS = sgn(psi_t) * psi_{t,LCFS} = <sign_toroidal_flux> * <edge_toroidal_flux_over_2pi>
+      gbdrift0_psi = sign_toroidal_flux * edge_toroidal_flux_over_2pi * (B_sub_theta_vmec * d_B_d_zeta &
+                                                                         - B_sub_zeta * d_B_d_theta_vmec) / sqrt_g * 2 * shat / (B * B * B)
+
 !    gbdrift0 = (B_sub_theta_vmec * d_B_d_zeta - B_sub_zeta * d_B_d_theta_vmec) / sqrt_g * edge_toroidal_flux_over_2pi &
 !    gbdrift0 = abs(edge_toroidal_flux_over_2pi) &
 !         * (B_sub_theta_vmec * d_B_d_zeta - B_sub_zeta * d_B_d_theta_vmec) / sqrt_g &

init_g.f90

@@ -391,6 +391,7 @@ subroutine ginit_noise
       use mp, only: scope, crossdomprocs, subprocs
       use file_utils, only: runtype_option_switch, runtype_multibox
       use physics_flags, only: nonlinear
+      use stella_geometry, only: sign_torflux
       use ran
 
       implicit none
@@ -432,8 +433,10 @@ subroutine ginit_noise
             do iky = 1, naky
                do it = 1, ntubes
                   do iz = -nzgrid, nzgrid
+
+                     ! For the same rng-seed, the <-sign_torflux> will make the time trace of CCW and CW more similar
                      a = ranf() - 0.5
-                     b = ranf() - 0.5
+                     b = -sign_torflux * (ranf() - 0.5)
                      ! do not populate high k modes with large amplitudes
                      if ((ikx > 1 .or. iky > 1) .and. (kperp2(iky, ikx, ia, iz) >= kmin)) then
                         !the following as an extra factor of kmin to offset the Gamma-1 in quasineutrality