From b44e6c08d590abfd8ec60d12677c6e2cfd381caf Mon Sep 17 00:00:00 2001 From: Michael Barnes Date: Fri, 5 Jul 2024 15:47:55 +0100 Subject: [PATCH] changed from gds2, gds21 and gds22 to grad_y_dot_grad_y, grad_x_dot_grad_y and grad_x_dot_grad_x. also, implemented a simple z-pinch equilibrium geometry option. --- CMakeLists.txt | 1 + dist_fn.f90 | 18 +-- geo/millerlocal.f90 | 64 ++++---- geo/stella_geometry.f90 | 142 +++++++++++------- .../compare_GIST_to_direct_vmec_gs2_interface | 40 ++--- ..._vmec_to_stella_geometry_interface_results | 36 ++--- ...test_vmec_to_stella_geometry_interface.f90 | 28 ++-- geo/zpinch.f90 | 82 ++++++++++ post_processing/stella_data.py | 6 +- stella_diagnostics.f90 | 56 +++++-- stella_io.fpp | 8 +- 11 files changed, 312 insertions(+), 169 deletions(-) create mode 100644 geo/zpinch.f90 diff --git a/CMakeLists.txt b/CMakeLists.txt index 47f505246b..cb847dc9d7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -75,6 +75,7 @@ set(STELLA_SOURCES_f90 geo/vmec_geo.f90 geo/vmec_interface/vmec_to_stella_geometry_interface.f90 geo/vmec_interface/fzero.f90 + geo/zpinch.f90 ) # Sources that _do_ need preprocessing diff --git a/dist_fn.f90 b/dist_fn.f90 index 20deadd683..e5c8f10613 100644 --- a/dist_fn.f90 +++ b/dist_fn.f90 @@ -127,10 +127,10 @@ end subroutine init_dist_fn subroutine init_kperp2 use dist_fn_arrays, only: kperp2 - use stella_geometry, only: gds2, gds21, gds22 + use stella_geometry, only: grad_y_dot_grad_y, grad_x_dot_grad_y, grad_x_dot_grad_x use stella_geometry, only: geo_surf, q_as_x use zgrid, only: nzgrid - use kt_grids, only: naky, nakx, theta0 + use kt_grids, only: naky, nakx use kt_grids, only: akx, aky use kt_grids, only: zonal_mode use kt_grids, only: nalpha @@ -146,19 +146,15 @@ subroutine init_kperp2 allocate (kperp2(naky, nakx, nalpha, -nzgrid:nzgrid)) do iky = 1, naky - if (zonal_mode(iky)) then + if (zonal_mode(iky) .and. q_as_x) then do ikx = 1, nakx - if (q_as_x) then - kperp2(iky, ikx, :, :) = akx(ikx) * akx(ikx) * gds22 - else - kperp2(iky, ikx, :, :) = akx(ikx) * akx(ikx) * gds22 / (geo_surf%shat**2) - end if + kperp2(iky, ikx, :, :) = akx(ikx) * akx(ikx) * grad_x_dot_grad_x * geo_surf%shat**2 end do else do ikx = 1, nakx - kperp2(iky, ikx, :, :) = aky(iky) * aky(iky) & - * (gds2 + 2.0 * theta0(iky, ikx) * gds21 & - + theta0(iky, ikx) * theta0(iky, ikx) * gds22) + kperp2(iky, ikx, :, :) = aky(iky) * aky(iky) * grad_y_dot_grad_y & + + 2.0 * aky(iky) * akx(ikx) * grad_x_dot_grad_y & + + akx(ikx) * akx(ikx) * grad_x_dot_grad_x end do end if end do diff --git a/geo/millerlocal.f90 b/geo/millerlocal.f90 index bea47de47f..498b628261 100644 --- a/geo/millerlocal.f90 +++ b/geo/millerlocal.f90 @@ -32,7 +32,7 @@ module millerlocal real :: bi, dqdr, d2Idr2 real, dimension(:), allocatable :: grho, bmag, grho_psi0, bmag_psi0, gradpar real, dimension(:), allocatable :: gradpararc, arc - real, dimension(:), allocatable :: gds2, gds21, gds22 + real, dimension(:), allocatable :: grad_y_dot_grad_y, grad_x_dot_grad_y, grad_x_dot_grad_x real, dimension(:), allocatable :: gds23, gds24 real, dimension(:), allocatable :: gbdrift0, gbdrift real, dimension(:), allocatable :: cvdrift0, cvdrift @@ -353,7 +353,7 @@ end subroutine communicate_parameters_multibox subroutine get_local_geo(nzed, nzgrid, zed_in, zed_equal_arc, & dpsidrho_out, dpsidrho_psi0_out, dIdrho_out, grho_out, & bmag_out, bmag_psi0_out, & - gds2_out, gds21_out, gds22_out, gds23_out, gds24_out, gradpar_out, & + grad_y_dot_grad_y_out, grad_x_dot_grad_y_out, grad_x_dot_grad_x_out, gds23_out, gds24_out, gradpar_out, & gbdrift0_out, gbdrift_out, cvdrift0_out, cvdrift_out, & dBdrho_out, d2Bdrdth_out, dgradpardrho_out, & btor_out, rmajor_out, & @@ -374,7 +374,7 @@ subroutine get_local_geo(nzed, nzgrid, zed_in, zed_equal_arc, & real, intent(out) :: dpsidrho_out, dpsidrho_psi0_out, dIdrho_out real, dimension(-nzgrid:), intent(out) :: grho_out, & bmag_out, bmag_psi0_out, & - gds2_out, gds21_out, gds22_out, gds23_out, gds24_out, & + grad_y_dot_grad_y_out, grad_x_dot_grad_y_out, grad_x_dot_grad_x_out, gds23_out, gds24_out, & gradpar_out, gbdrift0_out, & gbdrift_out, cvdrift0_out, cvdrift_out, & dBdrho_out, d2Bdrdth_out, dgradpardrho_out, & @@ -555,7 +555,7 @@ subroutine get_local_geo(nzed, nzgrid, zed_in, zed_equal_arc, & ! get |grad theta|^2, grad r . grad theta, grad alpha . grad theta, etc. call get_graddotgrad(dpsidrho, grho) - call get_gds(gds2, gds21, gds22, gds23, gds24) + call get_gds(grad_y_dot_grad_y, grad_x_dot_grad_y, grad_x_dot_grad_x, gds23, gds24) ! this is (grad alpha x B) . grad theta cross = dpsidrho * (gradrho_gradalph * gradalph_gradthet - gradalph2 * gradrho_gradthet) @@ -635,11 +635,11 @@ subroutine get_local_geo(nzed, nzgrid, zed_in, zed_equal_arc, & call geo_spline(theta, grho_psi0, zed_arc, grho_out) !grho is used to normalize fluxes call geo_spline(theta, bmag, zed_arc, bmag_out) call geo_spline(theta, bmag_psi0, zed_arc, bmag_psi0_out) - call geo_spline(theta, gds2, zed_arc, gds2_out) - call geo_spline(theta, gds21, zed_arc, gds21_out) - call geo_spline(theta, gds22, zed_arc, gds22_out) - call geo_spline(theta, gds21, zed_arc, gds23_out) - call geo_spline(theta, gds21, zed_arc, gds24_out) + call geo_spline(theta, grad_y_dot_grad_y, zed_arc, grad_y_dot_grad_y_out) + call geo_spline(theta, grad_x_dot_grad_y, zed_arc, grad_x_dot_grad_y_out) + call geo_spline(theta, grad_x_dot_grad_x, zed_arc, grad_x_dot_grad_x_out) + call geo_spline(theta, gds23, zed_arc, gds23_out) + call geo_spline(theta, gds24, zed_arc, gds24_out) call geo_spline(theta, gradpararc, zed_arc, gradpar_out) call geo_spline(theta, gbdrift, zed_arc, gbdrift_out) call geo_spline(theta, gbdrift0, zed_arc, gbdrift0_out) @@ -663,11 +663,11 @@ subroutine get_local_geo(nzed, nzgrid, zed_in, zed_equal_arc, & call geo_spline(theta, grho_psi0, zed_in, grho_out) !grho is used to normalize fluxes call geo_spline(theta, bmag, zed_in, bmag_out) call geo_spline(theta, bmag_psi0, zed_in, bmag_psi0_out) - call geo_spline(theta, gds2, zed_in, gds2_out) - call geo_spline(theta, gds21, zed_in, gds21_out) - call geo_spline(theta, gds22, zed_in, gds22_out) - call geo_spline(theta, gds21, zed_in, gds23_out) - call geo_spline(theta, gds21, zed_in, gds24_out) + call geo_spline(theta, grad_y_dot_grad_y, zed_in, grad_y_dot_grad_y_out) + call geo_spline(theta, grad_x_dot_grad_y, zed_in, grad_x_dot_grad_y_out) + call geo_spline(theta, grad_x_dot_grad_x, zed_in, grad_x_dot_grad_x_out) + call geo_spline(theta, gds23, zed_in, gds23_out) + call geo_spline(theta, gds24, zed_in, gds24_out) call geo_spline(theta, gradpar, zed_in, gradpar_out) call geo_spline(theta, gbdrift, zed_in, gbdrift_out) call geo_spline(theta, gbdrift0, zed_in, gbdrift0_out) @@ -722,8 +722,8 @@ subroutine get_local_geo(nzed, nzgrid, zed_in, zed_equal_arc, & '31.grgalph', '32.dgagr', '33.galph2', '34.dga2', '35.cross', & '36.dcrossdr', '37.gbdrift0', '38.dgbdrift0', '39.cvdrift0', '40.dcvdrift0', & '41.gbdrift', '42.dgbdrift', '43.cvdrift', '44.dcvdrift', '45.drzdth', & - '46.gradpar', '47.dgpardr', '48.gradparB', '49.dgparBdr', '50.gds2', & - '51.dgds2dr', '52.gds21', '53.dgds21dr', '54.gds22', '55.dgds22dr', & + '46.gradpar', '47.dgpardr', '48.gradparB', '49.dgparBdr', '50.grad_y_dot_grad_y', & + '51.dgds2dr', '52.grad_x_dot_grad_y', '53.dgds21dr', '54.grad_x_dot_grad_x', '55.dgds22dr', & '56.gds23', '57.gds24', '58.Zr' do i = -nz, nz @@ -736,8 +736,8 @@ subroutine get_local_geo(nzed, nzgrid, zed_in, zed_equal_arc, & gradrho_gradalph(i), dgagr(i), gradalph2(i), dga2(i), cross(i), & dcrossdr(i), gbdrift0(i), dgbdrift0drho(i), cvdrift0(i), dcvdrift0drho(i), & gbdrift(i), dgbdriftdrho(i), cvdrift(i), dcvdriftdrho(i), drzdth(i), & - gradpar(i), dgradpardrho(i), gradparB(i), dgradparBdrho(i), gds2(i), & - dgds2dr(i), gds21(i), dgds21dr(i), gds22(i), dgds22dr(i), gds23(i), gds24(i), & + gradpar(i), dgradpardrho(i), gradparB(i), dgradparBdrho(i), grad_y_dot_grad_y(i), & + dgds2dr(i), grad_x_dot_grad_y(i), dgds21dr(i), grad_x_dot_grad_x(i), dgds22dr(i), gds23(i), gds24(i), & Zr(2, i) end do close (1001) @@ -754,7 +754,7 @@ subroutine allocate_arrays(nr, nz) ! periodic quantities can be computed on 2*pi grid and replicated allocate (grho(-nz:nz), bmag(-nz:nz), gradpar(-nz:nz)) - allocate (gds2(-nz:nz), gds21(-nz:nz), gds22(-nz:nz), gds23(-nz:nz), gds24(-nz:nz)) + allocate (grad_y_dot_grad_y(-nz:nz), grad_x_dot_grad_y(-nz:nz), grad_x_dot_grad_x(-nz:nz), gds23(-nz:nz), gds24(-nz:nz)) allocate (gbdrift0(-nz:nz), gbdrift(-nz:nz)) allocate (cvdrift0(-nz:nz), cvdrift(-nz:nz)) allocate (Rr(nr, -nz:nz), Zr(nr, -nz:nz)) @@ -791,9 +791,9 @@ subroutine deallocate_arrays deallocate (bmag) deallocate (gradpar) - deallocate (gds2) - deallocate (gds21) - deallocate (gds22) + deallocate (grad_y_dot_grad_y) + deallocate (grad_x_dot_grad_y) + deallocate (grad_x_dot_grad_x) deallocate (gds23) deallocate (gds24) deallocate (gbdrift0) @@ -1029,25 +1029,29 @@ subroutine get_graddotgrad(dpsidrho, grho) end subroutine get_graddotgrad - subroutine get_gds(gds2, gds21, gds22, gds23, gds24) + subroutine get_gds(grad_y_dot_grad_y, grad_x_dot_grad_y, grad_x_dot_grad_x, gds23, gds24) implicit none - real, dimension(-nz:), intent(out) :: gds2, gds21, gds22, gds23, gds24 + real, dimension(-nz:), intent(out) :: grad_y_dot_grad_y, grad_x_dot_grad_y, grad_x_dot_grad_x, gds23, gds24 + ! NB: dx/dpsi = q/(r*Br), dy/dalpha = (dpsi/dr) / Br and psi_norm = psi/(Br*a^2) + ! |grad alpha|^2 * (dpsiN/drho)^2 (dpsiN/drho factor accounts for ky normalization) - gds2 = gradalph2 * dpsidrho_psi0**2 - ! (grad q . grad alpha) * (dpsiN/drho)^2 - gds21 = gradrho_gradalph * dqdr * dpsidrho_psi0**2 - ! |grad q|^2 * (dpsiN/drho)^2 - gds22 = (grho * dpsidrho_psi0 * dqdr)**2 + grad_y_dot_grad_y = gradalph2 * dpsidrho_psi0**2 + ! grad x . grad y = (grad_norm rho . grad_norm alpha) * (dx_norm/drho) * (dy_norm/dalpha), + ! with dx_norm/drho = (q/r*Bref) * (dpsi/drho) / a = (q/rho) * (dpsi_norm/drho) + ! and dy_norm/dalpha = dpsi_norm/drho + grad_x_dot_grad_y = gradrho_gradalph * dpsidrho_psi0**2 * qinp / rhoc + ! |grad x|^2 = |grad_norm rho * dx_norm/drho|^2 = (grho * q/(r*Br*a) * dpsi/drho)^2 + grad_x_dot_grad_x = (grho * qinp / rhoc * dpsidrho_psi0)**2 ! (grad rho . grad theta * |grad alpha|^2 - grad alpha . grad theta * grad rho . grad alpha) * (dpsiN/drho)^2 / B^2 gds23 = (gradrho_gradthet * gradalph2 - gradalph_gradthet * gradrho_gradalph) * (dpsidrho_psi0 / bmag)**2 ! (grad rho . grad theta * grad rho . grad alpha - grad alpha . grad theta * |grad rho|^2) * (dpsiN/drho)^2 / B^2 * q/rho gds24 = (gradrho_gradthet * gradrho_gradalph - gradalph_gradthet * grho**2) & * (dpsidrho_psi0 / bmag)**2 * (local%qinp_psi0 / local%rhoc_psi0) - ! note that kperp2 = (n0/a)^2*(drho/dpsiN)^2*(gds2 + 2*theta0*gds21 + theta0^2*gds22) + ! note that kperp2 = (n0/a)^2*(drho/dpsiN)^2*(grad_y_dot_grad_y + 2*theta0*gds21 + theta0^2*gds22) ! theta0 = kx/(ky*shat) end subroutine get_gds diff --git a/geo/stella_geometry.f90 b/geo/stella_geometry.f90 index 37d80bc4cb..a176cec226 100644 --- a/geo/stella_geometry.f90 +++ b/geo/stella_geometry.f90 @@ -13,7 +13,7 @@ module stella_geometry public :: gbdrift, gbdrift0 public :: dcvdriftdrho, dcvdrift0drho public :: dgbdriftdrho, dgbdrift0drho - public :: gds2, gds21, gds22, gds23, gds24, gds25, gds26 + public :: grad_y_dot_grad_y, grad_x_dot_grad_y, grad_x_dot_grad_x, gds23, gds24, gds25, gds26 public :: dgds2dr, dgds21dr, dgds22dr public :: exb_nonlin_fac, exb_nonlin_fac_p public :: jacob, djacdrho @@ -59,7 +59,7 @@ module stella_geometry real, dimension(:, :), allocatable :: gbdrift, gbdrift0 real, dimension(:, :), allocatable :: dcvdriftdrho, dcvdrift0drho real, dimension(:, :), allocatable :: dgbdriftdrho, dgbdrift0drho - real, dimension(:, :), allocatable :: gds2, gds21, gds22, gds23, gds24, gds25, gds26 + real, dimension(:, :), allocatable :: grad_y_dot_grad_y, grad_x_dot_grad_y, grad_x_dot_grad_x, gds23, gds24, gds25, gds26 real, dimension(:, :), allocatable :: dgds2dr, dgds21dr real, dimension(:, :), allocatable :: dgds22dr real, dimension(:, :), allocatable :: theta_vmec @@ -78,10 +78,11 @@ module stella_geometry integer, parameter :: geo_option_inputprof = 2 integer, parameter :: geo_option_vmec = 3 integer, parameter :: geo_option_multibox = 4 + integer, parameter :: geo_option_zpinch = 5 logical :: overwrite_geometry logical :: overwrite_bmag, overwrite_gradpar - logical :: overwrite_gds2, overwrite_gds21, overwrite_gds22 + logical :: overwrite_grad_y_dot_grad_y, overwrite_grad_x_dot_grad_y, overwrite_grad_x_dot_grad_x logical :: overwrite_gds23, overwrite_gds24 logical :: overwrite_gbdrift, overwrite_cvdrift, overwrite_gbdrift0 logical :: q_as_x @@ -99,6 +100,7 @@ subroutine init_geometry(nalpha, naky) use mp, only: proc0 use millerlocal, only: read_local_parameters, get_local_geo use millerlocal, only: communicate_parameters_multibox + use zpinch, only: get_zpinch_geometry_coefficients use vmec_geo, only: read_vmec_parameters, get_vmec_geo use vmec_to_stella_geometry_interface_mod, only: desired_zmin use inputprofiles_interface, only: read_inputprof_geo @@ -155,7 +157,7 @@ subroutine init_geometry(nalpha, naky) call get_local_geo(nzed, nzgrid, zed, zed_equal_arc, & dpsidrho, dpsidrho_psi0, dIdrho, grho(1, :), & bmag(1, :), bmag_psi0(1, :), & - gds2(1, :), gds21(1, :), gds22(1, :), & + grad_y_dot_grad_y(1, :), grad_x_dot_grad_y(1, :), grad_x_dot_grad_x(1, :), & gds23(1, :), gds24(1, :), gradpar, & gbdrift0(1, :), gbdrift(1, :), cvdrift0(1, :), cvdrift(1, :), & dBdrho, d2Bdrdth, dgradpardrho, btor, rmajor, & @@ -195,6 +197,37 @@ subroutine init_geometry(nalpha, naky) ! aref and bref should not be needed, so set to 1 aref = 1.0; bref = 1.0 zeta(1, :) = zed * geo_surf%qinp + case (geo_option_zpinch) + call allocate_arrays(nalpha, nzgrid) + call get_zpinch_geometry_coefficients(nzgrid, bmag(1, :), gradpar, grho(1, :), & + grad_y_dot_grad_y(1, :), grad_x_dot_grad_y(1, :), grad_x_dot_grad_x(1, :), & + gbdrift0(1, :), gbdrift(1, :), cvdrift0(1, :), cvdrift(1, :), btor, rmajor) + + sign_torflux = -1 + ! effectively choose psi = x * B * a_ref = x * B * r0 + dpsidrho = 1.0 ; dpsidrho_psi0 = 1.0 + bmag_psi0 = bmag + !> b_dot_grad_z is the alpha-dependent b . grad z, + !> and gradpar is the constant-in-alpha part of it. + !> for a z-pinch, b_dot_grad_z is independent of alpha. + b_dot_grad_z(1, :) = gradpar + ! note that psi here is meaningless + drhodpsi = 1./dpsidrho + drhodpsi_psi0 = 1./dpsidrho_psi0 + ! dxdXcoord = a*Bref*dx/dpsi = sign(dx/dpsi) * a*q/r + dxdXcoord_sign = 1 + dxdXcoord = 1.0 + ! dydalpha = (dy/dalpha) / a = sign(dydalpha) * (dpsi/dr) / (a*Bref) + dydalpha_sign = 1 + dydalpha = dydalpha_sign * dpsidrho + grad_x = sqrt(grad_x_dot_grad_x) + ! there is no magnetic shear in the z-pinch and thus no need for twist-and-shift + twist_and_shift_geo_fac = 1.0 + ! aref and bref should not be needed, so set to 1 + aref = 1.0; bref = 1.0 + ! zeta should not be needed + zeta(1, :) = 0.0 + case (geo_option_multibox) ! read in Miller local parameters call read_local_parameters(nzed, nzgrid, geo_surf) @@ -209,7 +242,7 @@ subroutine init_geometry(nalpha, naky) call get_local_geo(nzed, nzgrid, zed, zed_equal_arc, & dpsidrho, dpsidrho_psi0, dIdrho, grho(1, :), & bmag(1, :), bmag_psi0(1, :), & - gds2(1, :), gds21(1, :), gds22(1, :), & + grad_y_dot_grad_y(1, :), grad_x_dot_grad_y(1, :), grad_x_dot_grad_x(1, :), & gds23(1, :), gds24(1, :), gradpar, & gbdrift0(1, :), gbdrift(1, :), cvdrift0(1, :), cvdrift(1, :), & dBdrho, d2Bdrdth, dgradpardrho, btor, rmajor, & @@ -263,7 +296,7 @@ subroutine init_geometry(nalpha, naky) call get_local_geo(nzed, nzgrid, zed, zed_equal_arc, & dpsidrho, dpsidrho_psi0, dIdrho, grho(1, :), & bmag(1, :), bmag_psi0(1, :), & - gds2(1, :), gds21(1, :), gds22(1, :), & + grad_y_dot_grad_y(1, :), grad_x_dot_grad_y(1, :), grad_x_dot_grad_x(1, :), & gds23(1, :), gds24(1, :), gradpar, & gbdrift0(1, :), gbdrift(1, :), cvdrift0(1, :), cvdrift(1, :), & dBdrho, d2Bdrdth, dgradpardrho, btor, rmajor, & @@ -394,18 +427,17 @@ subroutine init_geometry(nalpha, naky) !> grad_x = | grad x | grad_x = sqrt(abs(grad_psi_grad_psi * dxdXcoord**2)) - !> gds2 = |grad y|^2 = |grad alpha|^2 * (dy/dalpha)^2 + !> grad_y_dot_grad_y = |grad y|^2 = |grad alpha|^2 * (dy/dalpha)^2 !> note that rhotor = sqrt(psi/psi_LCFS) - gds2 = grad_alpha_grad_alpha * dydalpha**2 + grad_y_dot_grad_y = grad_alpha_grad_alpha * dydalpha**2 - !> Define = hat{s} ∇x . ∇y + !> Define grad_x_dot_grad_y = ∇x . ∇y !> Use (dx/dψ)*(dy/dα) = 1/(a ρ0 Bref) * (a ρ0) = 1/Bref !> Use ∇x . ∇y = (dx/dψ)(dy/dα) ∇ψ . ∇α = (1/Bref) ∇ψ . ∇α = - gds21 = grad_alpha_grad_psi * geo_surf%shat + grad_x_dot_grad_y = grad_alpha_grad_psi - !> 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 + !> grad_x_dot_grad_x = |grad x|^2 = |grad psi_t|^2 * (dx/dpsi_t)^2 + grad_x_dot_grad_x = grad_x**2 !> gbdrift_alpha and cvdrift_alpha contain !> the grad-B and curvature drifts projected onto @@ -588,7 +620,7 @@ subroutine overwrite_selected_geometric_coefficients(nalpha) integer :: ia, iz real :: bmag_file, gradpar_file - real :: gds2_file, gds21_file, gds22_file, gds23_file, gds24_file + real :: grad_y_dot_grad_y_file, grad_x_dot_grad_y_file, grad_x_dot_grad_x_file, gds23_file, gds24_file real :: gbdrift_file, cvdrift_file, gbdrift0_file call get_unused_unit(geofile_unit) @@ -598,19 +630,19 @@ subroutine overwrite_selected_geometric_coefficients(nalpha) read (geofile_unit, fmt=*) dum_char read (geofile_unit, fmt=*) dum_char - ! overwrite bmag, gradpar, gds2, gds21, gds22, gds23, gds24, gbdrift, cvdrift, gbdrift0, and cvdrift0 + ! overwrite bmag, gradpar, grad_y_dot_grad_y, grad_x_dot_grad_y, grad_x_dot_grad_x, gds23, gds24, gbdrift, cvdrift, gbdrift0, and cvdrift0 ! with values from file do ia = 1, nalpha do iz = -nzgrid, nzgrid read (geofile_unit, fmt='(13e12.4)') dum_real, dum_real, dum_real, bmag_file, gradpar_file, & - gds2_file, gds21_file, gds22_file, gds23_file, & + grad_y_dot_grad_y_file, grad_x_dot_grad_y_file, grad_x_dot_grad_x_file, gds23_file, & gds24_file, gbdrift_file, cvdrift_file, gbdrift0_file if (overwrite_bmag) bmag(ia, iz) = bmag_file if (overwrite_gradpar) gradpar(iz) = gradpar_file if (overwrite_gradpar) b_dot_grad_z(1, iz) = gradpar_file ! assuming we are only reading in for a single alpha. Usually, gradpar is the average of all b_dot_grad_z values. - if (overwrite_gds2) gds2(ia, iz) = gds2_file - if (overwrite_gds21) gds21(ia, iz) = gds21_file - if (overwrite_gds22) gds22(ia, iz) = gds22_file + if (overwrite_grad_y_dot_grad_y) grad_y_dot_grad_y(ia, iz) = grad_y_dot_grad_y_file + if (overwrite_grad_x_dot_grad_y) grad_x_dot_grad_y(ia, iz) = grad_x_dot_grad_y_file + if (overwrite_grad_x_dot_grad_x) grad_x_dot_grad_x(ia, iz) = grad_x_dot_grad_x_file if (overwrite_gds23) gds23(ia, iz) = gds23_file if (overwrite_gds24) gds24(ia, iz) = gds24_file if (overwrite_gbdrift) gbdrift(ia, iz) = gbdrift_file @@ -638,9 +670,9 @@ subroutine set_ffs_geo_coefs_constant(nalpha) call set_coef_constant(grho, nalpha) call set_coef_constant(bmag, nalpha) call set_coef_constant(bmag_psi0, nalpha) - call set_coef_constant(gds2, nalpha) - call set_coef_constant(gds21, nalpha) - call set_coef_constant(gds22, nalpha) + call set_coef_constant(grad_y_dot_grad_y, nalpha) + call set_coef_constant(grad_x_dot_grad_y, nalpha) + call set_coef_constant(grad_x_dot_grad_x, nalpha) call set_coef_constant(gds23, nalpha) call set_coef_constant(gds24, nalpha) call set_coef_constant(gds25, nalpha) @@ -677,28 +709,28 @@ subroutine allocate_arrays(nalpha, nzgrid) if (.not. allocated(bmag)) allocate (bmag(nalpha, -nzgrid:nzgrid)) if (.not. allocated(bmag_psi0)) allocate (bmag_psi0(nalpha, -nzgrid:nzgrid)) - if (.not. allocated(gds2)) allocate (gds2(nalpha, -nzgrid:nzgrid)) - if (.not. allocated(gds21)) allocate (gds21(nalpha, -nzgrid:nzgrid)) - if (.not. allocated(gds22)) allocate (gds22(nalpha, -nzgrid:nzgrid)) - if (.not. allocated(gds23)) allocate (gds23(nalpha, -nzgrid:nzgrid)) - if (.not. allocated(gds24)) allocate (gds24(nalpha, -nzgrid:nzgrid)) - if (.not. allocated(gds25)) allocate (gds25(nalpha, -nzgrid:nzgrid)) - if (.not. allocated(gds26)) allocate (gds26(nalpha, -nzgrid:nzgrid)) - if (.not. allocated(dgds2dr)) allocate (dgds2dr(nalpha, -nzgrid:nzgrid)) - if (.not. allocated(dgds21dr)) allocate (dgds21dr(nalpha, -nzgrid:nzgrid)) - if (.not. allocated(dgds22dr)) allocate (dgds22dr(nalpha, -nzgrid:nzgrid)) + if (.not. allocated(grad_y_dot_grad_y)) allocate (grad_y_dot_grad_y(nalpha, -nzgrid:nzgrid)) + if (.not. allocated(grad_x_dot_grad_y)) allocate (grad_x_dot_grad_y(nalpha, -nzgrid:nzgrid)) + if (.not. allocated(grad_x_dot_grad_x)) allocate (grad_x_dot_grad_x(nalpha, -nzgrid:nzgrid)) + if (.not. allocated(gds23)) allocate (gds23(nalpha, -nzgrid:nzgrid)) ; gds23 = 0.0 + if (.not. allocated(gds24)) allocate (gds24(nalpha, -nzgrid:nzgrid)) ; gds24 = 0.0 + if (.not. allocated(gds25)) allocate (gds25(nalpha, -nzgrid:nzgrid)) ; gds25 = 0.0 + if (.not. allocated(gds26)) allocate (gds26(nalpha, -nzgrid:nzgrid)) ; gds26 = 0.0 + if (.not. allocated(dgds2dr)) allocate (dgds2dr(nalpha, -nzgrid:nzgrid)) ; dgds2dr = 0.0 + if (.not. allocated(dgds21dr)) allocate (dgds21dr(nalpha, -nzgrid:nzgrid)) ; dgds21dr = 0.0 + if (.not. allocated(dgds22dr)) allocate (dgds22dr(nalpha, -nzgrid:nzgrid)) ; dgds22dr = 0.0 if (.not. allocated(gbdrift)) allocate (gbdrift(nalpha, -nzgrid:nzgrid)) if (.not. allocated(gbdrift0)) allocate (gbdrift0(nalpha, -nzgrid:nzgrid)) if (.not. allocated(cvdrift)) allocate (cvdrift(nalpha, -nzgrid:nzgrid)) if (.not. allocated(cvdrift0)) allocate (cvdrift0(nalpha, -nzgrid:nzgrid)) - if (.not. allocated(dgbdriftdrho)) allocate (dgbdriftdrho(nalpha, -nzgrid:nzgrid)) - if (.not. allocated(dcvdriftdrho)) allocate (dcvdriftdrho(nalpha, -nzgrid:nzgrid)) - if (.not. allocated(dgbdrift0drho)) allocate (dgbdrift0drho(nalpha, -nzgrid:nzgrid)) - if (.not. allocated(dcvdrift0drho)) allocate (dcvdrift0drho(nalpha, -nzgrid:nzgrid)) + if (.not. allocated(dgbdriftdrho)) allocate (dgbdriftdrho(nalpha, -nzgrid:nzgrid)) ; dgbdriftdrho = 0.0 + if (.not. allocated(dcvdriftdrho)) allocate (dcvdriftdrho(nalpha, -nzgrid:nzgrid)) ; dcvdriftdrho = 0.0 + if (.not. allocated(dgbdrift0drho)) allocate (dgbdrift0drho(nalpha, -nzgrid:nzgrid)) ; dgbdrift0drho = 0.0 + if (.not. allocated(dcvdrift0drho)) allocate (dcvdrift0drho(nalpha, -nzgrid:nzgrid)) ; dcvdrift0drho = 0.0 if (.not. allocated(dbdzed)) allocate (dbdzed(nalpha, -nzgrid:nzgrid)) if (.not. allocated(theta_vmec)) allocate (theta_vmec(nalpha, -nzgrid:nzgrid)) if (.not. allocated(jacob)) allocate (jacob(nalpha, -nzgrid:nzgrid)) - if (.not. allocated(djacdrho)) allocate (djacdrho(nalpha, -nzgrid:nzgrid)) + if (.not. allocated(djacdrho)) allocate (djacdrho(nalpha, -nzgrid:nzgrid)) ; djacdrho = 0.0 if (.not. allocated(grho)) allocate (grho(nalpha, -nzgrid:nzgrid)) if (.not. allocated(grad_x)) allocate (grad_x(nalpha, -nzgrid:nzgrid)) if (.not. allocated(dl_over_b)) allocate (dl_over_b(nalpha, -nzgrid:nzgrid)) @@ -709,9 +741,9 @@ subroutine allocate_arrays(nalpha, nzgrid) if (.not. allocated(zed_eqarc)) allocate (zed_eqarc(-nzgrid:nzgrid)) if (.not. allocated(btor)) allocate (btor(-nzgrid:nzgrid)) if (.not. allocated(rmajor)) allocate (rmajor(-nzgrid:nzgrid)) - if (.not. allocated(dBdrho)) allocate (dBdrho(-nzgrid:nzgrid)) - if (.not. allocated(d2Bdrdth)) allocate (d2Bdrdth(-nzgrid:nzgrid)) - if (.not. allocated(dgradpardrho)) allocate (dgradpardrho(-nzgrid:nzgrid)) + if (.not. allocated(dBdrho)) allocate (dBdrho(-nzgrid:nzgrid)) ; dBdrho = 0.0 + if (.not. allocated(d2Bdrdth)) allocate (d2Bdrdth(-nzgrid:nzgrid)) ; d2Bdrdth = 0.0 + if (.not. allocated(dgradpardrho)) allocate (dgradpardrho(-nzgrid:nzgrid)) ; dgradpardrho = 0.0 if (.not. allocated(alpha)) allocate (alpha(nalpha)); alpha = 0. if (.not. allocated(zeta)) allocate (zeta(nalpha, -nzgrid:nzgrid)); zeta = 0. @@ -734,23 +766,24 @@ subroutine read_parameters logical :: exist character(20) :: geo_option - type(text_option), dimension(5), parameter :: geoopts = (/ & + type(text_option), dimension(6), parameter :: geoopts = (/ & text_option('default', geo_option_local), & text_option('miller', geo_option_local), & text_option('local', geo_option_local), & + text_option('zpinch', geo_option_zpinch), & text_option('input.profiles', geo_option_inputprof), & text_option('vmec', geo_option_vmec)/) namelist /geo_knobs/ geo_option, geo_file, overwrite_bmag, overwrite_gradpar, & - overwrite_gds2, overwrite_gds21, overwrite_gds22, overwrite_gds23, overwrite_gds24, & + overwrite_grad_y_dot_grad_y, overwrite_grad_x_dot_grad_y, overwrite_grad_x_dot_grad_x, overwrite_gds23, overwrite_gds24, & overwrite_gbdrift, overwrite_cvdrift, overwrite_gbdrift0, q_as_x, set_bmag_const geo_option = 'local' overwrite_bmag = .false. overwrite_gradpar = .false. - overwrite_gds2 = .false. - overwrite_gds21 = .false. - overwrite_gds22 = .false. + overwrite_grad_y_dot_grad_y = .false. + overwrite_grad_x_dot_grad_y = .false. + overwrite_grad_x_dot_grad_x = .false. overwrite_gds23 = .false. overwrite_gds24 = .false. overwrite_gbdrift = .false. @@ -773,7 +806,7 @@ subroutine read_parameters end if overwrite_geometry = overwrite_bmag .or. overwrite_gradpar & - .or. overwrite_gds2 .or. overwrite_gds21 .or. overwrite_gds22 & + .or. overwrite_grad_y_dot_grad_y .or. overwrite_grad_x_dot_grad_y .or. overwrite_grad_x_dot_grad_x & .or. overwrite_gds23 .or. overwrite_gds24 & .or. overwrite_cvdrift .or. overwrite_gbdrift .or. overwrite_gbdrift0 @@ -800,9 +833,9 @@ subroutine broadcast_arrays call broadcast(rmajor) call broadcast(gradpar) call broadcast(b_dot_grad_z) - call broadcast(gds2) - call broadcast(gds21) - call broadcast(gds22) + call broadcast(grad_y_dot_grad_y) + call broadcast(grad_x_dot_grad_y) + call broadcast(grad_x_dot_grad_x) call broadcast(gds23) call broadcast(gds24) call broadcast(gds25) @@ -1015,13 +1048,13 @@ subroutine write_geometric_coefficients(nalpha) dxdXcoord, dydalpha, exb_nonlin_fac, exb_nonlin_fac_p * exb_nonlin_fac write (geometry_unit, *) - write (geometry_unit, '(15a12)') '# alpha', 'zed', 'zeta', 'bmag', 'bdot_grad_z', 'gds2', 'gds21', 'gds22', & + write (geometry_unit, '(15a12)') '# alpha', 'zed', 'zeta', 'bmag', 'bdot_grad_z', 'grad_y_dot_grad_y', 'grad_x_dot_grad_y', 'grad_x_dot_grad_x', & 'gds23', 'gds24', 'gbdrift', 'cvdrift', 'gbdrift0', 'bmag_psi0', 'btor' do ia = 1, nalpha do iz = -nzgrid, nzgrid ! write (geometry_unit,'(15e12.4)') alpha(ia), zed(iz), zeta(ia,iz), bmag(ia,iz), gradpar(iz), & write (geometry_unit, '(15e12.4)') alpha(ia), zed(iz), zeta(ia, iz), bmag(ia, iz), b_dot_grad_z(ia, iz), & - gds2(ia, iz), gds21(ia, iz), gds22(ia, iz), gds23(ia, iz), & + grad_y_dot_grad_y(ia, iz), grad_x_dot_grad_y(ia, iz), grad_x_dot_grad_x(ia, iz), gds23(ia, iz), & gds24(ia, iz), gbdrift(ia, iz), cvdrift(ia, iz), gbdrift0(ia, iz), & bmag_psi0(ia, iz), btor(iz) end do @@ -1043,6 +1076,7 @@ subroutine finish_init_geometry select case (geo_option_switch) case (geo_option_local) call finish_local_geo + case (geo_option_zpinch) case (geo_option_multibox) call finish_local_geo case (geo_option_inputprof) @@ -1071,9 +1105,9 @@ subroutine finish_geometry if (allocated(b_dot_grad_z)) deallocate (b_dot_grad_z) if (allocated(dl_over_b)) deallocate (dl_over_b) if (allocated(d_dl_over_b_drho)) deallocate (d_dl_over_b_drho) - if (allocated(gds2)) deallocate (gds2) - if (allocated(gds21)) deallocate (gds21) - if (allocated(gds22)) deallocate (gds22) + if (allocated(grad_y_dot_grad_y)) deallocate (grad_y_dot_grad_y) + if (allocated(grad_x_dot_grad_y)) deallocate (grad_x_dot_grad_y) + if (allocated(grad_x_dot_grad_x)) deallocate (grad_x_dot_grad_x) if (allocated(gds23)) deallocate (gds23) if (allocated(gds24)) deallocate (gds24) if (allocated(gds25)) deallocate (gds25) diff --git a/geo/vmec_interface/compare_GIST_to_direct_vmec_gs2_interface b/geo/vmec_interface/compare_GIST_to_direct_vmec_gs2_interface index b67d1c77fa..6481d437f4 100755 --- a/geo/vmec_interface/compare_GIST_to_direct_vmec_gs2_interface +++ b/geo/vmec_interface/compare_GIST_to_direct_vmec_gs2_interface @@ -45,9 +45,9 @@ data = np.loadtxt(f) gist_theta = data[:,0] gist_B = data[:,1] gist_gradpar = data[:,2] -gist_gds2 = data[:,3] -gist_gds21 = data[:,4] -gist_gds22 = data[:,5] +gist_grad_y_dot_grad_y = data[:,3] +gist_grad_x_dot_grad_y = data[:,4] +gist_grad_x_dot_grad_x = data[:,5] gist_cvdrift = data[:,6] gist_cvdrift0 = data[:,7] gist_gbdrift = data[:,8] @@ -128,8 +128,8 @@ for j in range(next_index, next_index+nalpha): splitline = lines[j].split() data.append([float(item) for item in splitline]) -gds2 = np.array(data) -if verbose: print "gds2:",gds2 +grad_y_dot_grad_y = np.array(data) +if verbose: print "grad_y_dot_grad_y:",grad_y_dot_grad_y next_index += nalpha+1 data = [] @@ -137,8 +137,8 @@ for j in range(next_index, next_index+nalpha): splitline = lines[j].split() data.append([float(item) for item in splitline]) -gds21 = np.array(data) -if verbose: print "gds21:",gds21 +grad_x_dot_grad_y = np.array(data) +if verbose: print "grad_x_dot_grad_y:",grad_x_dot_grad_y next_index += nalpha+1 data = [] @@ -146,8 +146,8 @@ for j in range(next_index, next_index+nalpha): splitline = lines[j].split() data.append([float(item) for item in splitline]) -gds22 = np.array(data) -if verbose: print "gds22:",gds22 +grad_x_dot_grad_x = np.array(data) +if verbose: print "grad_x_dot_grad_x:",grad_x_dot_grad_x next_index += nalpha+1 data = [] @@ -208,8 +208,8 @@ linespec_gist='.-' linespec_new = '.-' markerSize = 3; -# Flipping the sign of (d x / d psi) causes gbdrift0=cvdrift0 and gds21 to flip sign. -# Furthermore, GIST gets the sign of gds21 backwards. +# Flipping the sign of (d x / d psi) causes gbdrift0=cvdrift0 and grad_x_dot_grad_y to flip sign. +# Furthermore, GIST gets the sign of grad_x_dot_grad_y backwards. plt.subplot(numRows,numCols,plotNum) plotNum += 1 @@ -228,24 +228,24 @@ plt.title('gradpar') plt.subplot(numRows,numCols,plotNum) plotNum += 1 -plt.plot(gist_zeta, gist_gds2, linespec_gist) -plt.plot(zeta, gds2[ialpha,:],linespec_new,label='New interface',ms=markerSize) +plt.plot(gist_zeta, gist_grad_y_dot_grad_y, linespec_gist) +plt.plot(zeta, grad_y_dot_grad_y[ialpha,:],linespec_new,label='New interface',ms=markerSize) plt.xlabel('zeta') -plt.title('gds2') +plt.title('grad_y_dot_grad_y') plt.subplot(numRows,numCols,plotNum) plotNum += 1 -plt.plot(gist_zeta, gist_gds21, linespec_gist) -plt.plot(zeta, gds21[ialpha,:],linespec_new,label='New interface',ms=markerSize) +plt.plot(gist_zeta, gist_grad_x_dot_grad_y, linespec_gist) +plt.plot(zeta, grad_x_dot_grad_y[ialpha,:],linespec_new,label='New interface',ms=markerSize) plt.xlabel('zeta') -plt.title('gds21') +plt.title('grad_x_dot_grad_y') plt.subplot(numRows,numCols,plotNum) plotNum += 1 -plt.plot(gist_zeta, gist_gds22, linespec_gist) -plt.plot(zeta, gds22[ialpha,:],linespec_new,label='New interface',ms=markerSize) +plt.plot(gist_zeta, gist_grad_x_dot_grad_x, linespec_gist) +plt.plot(zeta, grad_x_dot_grad_x[ialpha,:],linespec_new,label='New interface',ms=markerSize) plt.xlabel('zeta') -plt.title('gds22') +plt.title('grad_x_dot_grad_x') plt.subplot(numRows,numCols,plotNum) plotNum += 1 diff --git a/geo/vmec_interface/plot_vmec_to_stella_geometry_interface_results b/geo/vmec_interface/plot_vmec_to_stella_geometry_interface_results index 57c6c9a2d8..77cb3c96be 100755 --- a/geo/vmec_interface/plot_vmec_to_stella_geometry_interface_results +++ b/geo/vmec_interface/plot_vmec_to_stella_geometry_interface_results @@ -52,8 +52,8 @@ for j in range(next_index, next_index+nalpha): splitline = lines[j].split() data.append([float(item) for item in splitline]) -gds2 = np.array(data) -if verbose: print "gds2:",gds2 +grad_y_dot_grad_y = np.array(data) +if verbose: print "grad_y_dot_grad_y:",grad_y_dot_grad_y next_index += nalpha+1 data = [] @@ -61,8 +61,8 @@ for j in range(next_index, next_index+nalpha): splitline = lines[j].split() data.append([float(item) for item in splitline]) -gds21 = np.array(data) -if verbose: print "gds21:",gds21 +grad_x_dot_grad_y = np.array(data) +if verbose: print "grad_x_dot_grad_y:",grad_x_dot_grad_y next_index += nalpha+1 data = [] @@ -70,8 +70,8 @@ for j in range(next_index, next_index+nalpha): splitline = lines[j].split() data.append([float(item) for item in splitline]) -gds22 = np.array(data) -if verbose: print "gds22:",gds22 +grad_x_dot_grad_x = np.array(data) +if verbose: print "grad_x_dot_grad_x:",grad_x_dot_grad_x next_index += nalpha+1 data = [] @@ -130,18 +130,18 @@ plt.title('gradpar') plt.xlabel('zeta') plt.subplot(numRows,numCols,3) -plt.plot(zeta,gds2[ialpha,:],'.-') -plt.title('gds2') +plt.plot(zeta,grad_y_dot_grad_y[ialpha,:],'.-') +plt.title('grad_y_dot_grad_y') plt.xlabel('zeta') plt.subplot(numRows,numCols,4) -plt.plot(zeta,gds21[ialpha,:],'.-') -plt.title('gds21') +plt.plot(zeta,grad_x_dot_grad_y[ialpha,:],'.-') +plt.title('grad_x_dot_grad_y') plt.xlabel('zeta') plt.subplot(numRows,numCols,5) -plt.plot(zeta,gds22[ialpha,:],'.-') -plt.title('gds22') +plt.plot(zeta,grad_x_dot_grad_x[ialpha,:],'.-') +plt.title('grad_x_dot_grad_x') plt.xlabel('zeta') plt.subplot(numRows,numCols,6) @@ -190,22 +190,22 @@ plt.ylabel('alpha') plt.colorbar() plt.subplot(numRows,numCols,3) -plt.contourf(zeta,alpha,gds2,numContours) -plt.title('gds2') +plt.contourf(zeta,alpha,grad_y_dot_grad_y,numContours) +plt.title('grad_y_dot_grad_y') plt.xlabel('zeta') plt.ylabel('alpha') plt.colorbar() plt.subplot(numRows,numCols,4) -plt.contourf(zeta,alpha,gds21,numContours) -plt.title('gds21') +plt.contourf(zeta,alpha,grad_x_dot_grad_y,numContours) +plt.title('grad_x_dot_grad_y') plt.xlabel('zeta') plt.ylabel('alpha') plt.colorbar() plt.subplot(numRows,numCols,5) -plt.contourf(zeta,alpha,gds22,numContours) -plt.title('gds22') +plt.contourf(zeta,alpha,grad_x_dot_grad_x,numContours) +plt.title('grad_x_dot_grad_x') plt.xlabel('zeta') plt.ylabel('alpha') plt.colorbar() diff --git a/geo/vmec_interface/test_vmec_to_stella_geometry_interface.f90 b/geo/vmec_interface/test_vmec_to_stella_geometry_interface.f90 index 36778896d3..dc7b48f75b 100644 --- a/geo/vmec_interface/test_vmec_to_stella_geometry_interface.f90 +++ b/geo/vmec_interface/test_vmec_to_stella_geometry_interface.f90 @@ -28,7 +28,7 @@ program test_vmec_to_stella_geometry_interface integer :: sign_toroidal_flux real, dimension(nalpha) :: alpha real, dimension(-nzgrid:nzgrid) :: zeta - real, dimension(nalpha, -nzgrid:nzgrid) :: bmag, gradpar, gds2, gds21, gds22, gds23, gds24, gds25, gds26 + real, dimension(nalpha, -nzgrid:nzgrid) :: bmag, gradpar, grad_y_dot_grad_y, grad_x_dot_grad_y, grad_x_dot_grad_x, gds23, gds24, gds25, gds26 real, dimension(nalpha, -nzgrid:nzgrid) :: gbdrift, gbdrift0, cvdrift, cvdrift0 real, dimension(nalpha, -nzgrid:nzgrid) :: theta_vmec real, dimension(nalpha, -nzgrid:nzgrid) :: B_sub_zeta, B_sub_theta_vmec, displacement @@ -51,7 +51,7 @@ program test_vmec_to_stella_geometry_interface desired_normalized_toroidal_flux, vmec_surface_option, verbose, & normalized_toroidal_flux_used, safety_factor_q, shat, L_reference, B_reference, nfp, & sign_toroidal_flux, & - alpha, zeta, bmag, gradpar, gds2, gds21, gds22, gds23, gds24, gds25, gds26, & + alpha, zeta, bmag, gradpar, grad_y_dot_grad_y, grad_x_dot_grad_y, grad_x_dot_grad_x, gds23, gds24, gds25, gds26, & gbdrift, gbdrift0, cvdrift, cvdrift0, & theta_vmec, B_sub_zeta, B_sub_theta_vmec, displacement, gradpar_zeta_prefac, & ierr) @@ -89,19 +89,19 @@ program test_vmec_to_stella_geometry_interface print *, gradpar(j, :) end do - print *, "gds2:" + print *, "grad_y_dot_grad_y:" do j = 1, nalpha - print *, gds2(j, :) + print *, grad_y_dot_grad_y(j, :) end do - print *, "gds21:" + print *, "grad_x_dot_grad_y:" do j = 1, nalpha - print *, gds21(j, :) + print *, grad_x_dot_grad_y(j, :) end do - print *, "gds22:" + print *, "grad_x_dot_grad_x:" do j = 1, nalpha - print *, gds22(j, :) + print *, grad_x_dot_grad_x(j, :) end do print *, "gds23:" @@ -168,19 +168,19 @@ program test_vmec_to_stella_geometry_interface write (iunit, *) gradpar(j, :) end do - write (iunit, *) 'gds2' + write (iunit, *) 'grad_y_dot_grad_y' do j = 1, nalpha - write (iunit, *) gds2(j, :) + write (iunit, *) grad_y_dot_grad_y(j, :) end do - write (iunit, *) 'gds21' + write (iunit, *) 'grad_x_dot_grad_y' do j = 1, nalpha - write (iunit, *) gds21(j, :) + write (iunit, *) grad_x_dot_grad_y(j, :) end do - write (iunit, *) 'gds22' + write (iunit, *) 'grad_x_dot_grad_x' do j = 1, nalpha - write (iunit, *) gds22(j, :) + write (iunit, *) grad_x_dot_grad_x(j, :) end do write (iunit, *) 'gds23' diff --git a/geo/zpinch.f90 b/geo/zpinch.f90 new file mode 100644 index 0000000000..3ff95aee48 --- /dev/null +++ b/geo/zpinch.f90 @@ -0,0 +1,82 @@ +module zpinch + + implicit none + + public :: get_zpinch_geometry_coefficients + + private + +contains + + ! subroutine read_zpinch_parameters (radius) + + ! use file_utils, only: input_unit_exist + + ! implicit none + + ! real, intent (out) :: radius + + ! integer :: in_file + ! logical :: exist + + ! namelist /zpinch_parameters/ radius + + ! ! default value for radius + ! radius = 0.5 + + ! in_file = input_unit_exist("zpinch_parameters", exist) + ! if (exist) read (unit=in_file, nml=zpinch_parameters) + + ! end subroutine read_zpinch_parameters + + ! use Z-pinch equilibrium. + ! the parallel coordinate, z, is chosen to be the arc-length, + ! i.e., the physical poloidal angle times the radius of the chosen magnetic field, r0. + ! the radial coordinate, x, is directed away from the middle of the circle + ! and is normalised by the local radius r0. + ! the bi-normal coordinate, y, is chosen to form an orthogonal, right-handed coordinate + ! system with (y,x,z) and is also normalised by r0. + subroutine get_zpinch_geometry_coefficients (nzgrid, bmag, gradpar, grad_rho, & + grad_y_dot_grad_y, grad_x_dot_grad_y, grad_x_dot_grad_x, & + gbdrift0, gbdrift, cvdrift0, cvdrift, btor, rmajor) + + implicit none + + integer, intent(in) :: nzgrid + real, dimension (-nzgrid:), intent(out) :: bmag, gradpar, grad_rho, & + grad_y_dot_grad_y, grad_x_dot_grad_y, grad_x_dot_grad_x, & + gbdrift0, gbdrift, cvdrift0, cvdrift, btor, rmajor + + ! bmag = B(r0) / B_ref + ! as B is constant along radius r0, choose B_ref = B(r0), so bmag = 1 + bmag = 1.0 + ! gradpar = bhat . grad z = b . (r0*grad) theta = 1 + gradpar = 1.0 + ! grad_rho = | r0 * grad (r / r0) | = 1 + grad_rho = 1.0 + ! grad_y_dot_grad_y = 1 + grad_y_dot_grad_y = 1.0 + ! grad_x_dot_grad_y = 0 because (x,y) are orthogonal coordinates + grad_x_dot_grad_y = 0.0 + ! grad_x_dot_grad_x = 1 + grad_x_dot_grad_x = 1.0 + ! the x-component of the grad-B drift is proportional to gbdrift0; zero in a z-pinch + gbdrift0 = 0.0 + ! the x-component of the curvature drift is proportional to cvdrift0; zero in a z-pinch + cvdrift0 = 0.0 + ! gbdrift = 2 * bhat / B_norm x (grad_norm B_norm / B_norm) . grad y + ! = 2 * bhat / B_norm x (grad_norm x * d ln B_norm / dx) . grad y + ! = -2 * d(ln B_norm) / dx_norm = 2 * r0 / L_B = 1 (from MHD equilibrium with beta=0) + gbdrift = 2.0 + ! cvdrift = 2 * bhat / B_norm x (bhat . grad_norm bhat) . grad y = 2 + cvdrift = gbdrift + + ! btor and rmajor used for flow shear calculations; + ! choosing btor = 0 means all flow shear is perpendicular rather than parallel + btor = 0.0 + ! rmajor set to one to avoid divide by zero and otherwise is unused + rmajor = 1.0 + + end subroutine get_zpinch_geometry_coefficients + +end module zpinch diff --git a/post_processing/stella_data.py b/post_processing/stella_data.py index 3e4a15126a..3da8095397 100644 --- a/post_processing/stella_data.py +++ b/post_processing/stella_data.py @@ -51,9 +51,9 @@ gbdrift0 = np.copy(ncfile.variables['gbdrift0'][:]) cvdrift = np.copy(ncfile.variables['cvdrift'][:]) cvdrift0 = np.copy(ncfile.variables['cvdrift0'][:]) -gds2 = np.copy(ncfile.variables['gds2'][:]) -gds21 = np.copy(ncfile.variables['gds21'][:]) -gds22 = np.copy(ncfile.variables['gds22'][:]) +grad_y_dot_grad_y = np.copy(ncfile.variables['grad_y_dot_grad_y'][:]) +grad_x_dot_grad_y = np.copy(ncfile.variables['grad_x_dot_grad_y'][:]) +grad_x_dot_grad_x = np.copy(ncfile.variables['grad_x_dot_grad_x'][:]) def read_stella_float(var): diff --git a/stella_diagnostics.f90 b/stella_diagnostics.f90 index c0e9570cd4..483b70f202 100644 --- a/stella_diagnostics.f90 +++ b/stella_diagnostics.f90 @@ -513,7 +513,7 @@ subroutine get_fluxes(g, pflx, vflx, qflx, & use stella_layouts, only: iky_idx, ikx_idx, iz_idx, it_idx, is_idx use species, only: spec, nspec use stella_geometry, only: jacob, grho, bmag, btor - use stella_geometry, only: gds21, gds22 + use stella_geometry, only: grad_x_dot_grad_y, grad_x_dot_grad_x use stella_geometry, only: geo_surf use zgrid, only: delzed, nzgrid, ntubes use vpamu_grids, only: nvpa, nmu @@ -578,9 +578,12 @@ subroutine get_fluxes(g, pflx, vflx, qflx, & ! parallel component gtmp1 = g(:, :, ikxkyz) * spread(vpa, 2, nmu) * geo_surf%rmaj * btor(iz) / bmag(ia, iz) call gyro_average(gtmp1, ikxkyz, gtmp2) +! gtmp1 = -g(:, :, ikxkyz) * zi * aky(iky) * spread(vperp2(ia, iz, :), 1, nvpa) * geo_surf%rhoc & +! * (gds21(ia, iz) + theta0(iky, ikx) * gds22(ia, iz)) * spec(is)%smz & +! / (geo_surf%qinp * geo_surf%shat * bmag(ia, iz)**2) gtmp1 = -g(:, :, ikxkyz) * zi * aky(iky) * spread(vperp2(ia, iz, :), 1, nvpa) * geo_surf%rhoc & - * (gds21(ia, iz) + theta0(iky, ikx) * gds22(ia, iz)) * spec(is)%smz & - / (geo_surf%qinp * geo_surf%shat * bmag(ia, iz)**2) + * (grad_x_dot_grad_y(ia, iz) + theta0(iky, ikx) * geo_surf%shat * grad_x_dot_grad_x(ia, iz)) * spec(is)%smz & + / (geo_surf%qinp * bmag(ia, iz)**2) call gyro_average_j1(gtmp1, ikxkyz, gtmp3) gtmp1 = gtmp2 + gtmp3 @@ -613,10 +616,14 @@ subroutine get_fluxes(g, pflx, vflx, qflx, & * geo_surf%rmaj * btor(iz) / bmag(1, iz) call gyro_average(gtmp1, ikxkyz, gtmp2) ! perp component +! gtmp1 = spread(vpa, 2, nmu) * spec(is)%stm * g(:, :, ikxkyz) & +! * zi * aky(iky) * spread(vperp2(ia, iz, :), 1, nvpa) * geo_surf%rhoc & +! * (gds21(ia, iz) + theta0(iky, ikx) * gds22(ia, iz)) * spec(is)%smz & +! / (geo_surf%qinp * geo_surf%shat * bmag(ia, iz)**2) gtmp1 = spread(vpa, 2, nmu) * spec(is)%stm * g(:, :, ikxkyz) & - * zi * aky(iky) * spread(vperp2(ia, iz, :), 1, nvpa) * geo_surf%rhoc & - * (gds21(ia, iz) + theta0(iky, ikx) * gds22(ia, iz)) * spec(is)%smz & - / (geo_surf%qinp * geo_surf%shat * bmag(ia, iz)**2) + * zi * aky(iky) * spread(vperp2(ia, iz, :), 1, nvpa) * geo_surf%rhoc & + * (grad_x_dot_grad_y(ia, iz) + theta0(iky, ikx) * geo_surf%shat * grad_x_dot_grad_x(ia, iz)) * spec(is)%smz & + / (geo_surf%qinp * bmag(ia, iz)**2) call gyro_average_j1(gtmp1, ikxkyz, gtmp3) ! FLAG -- NEED TO ADD IN CONTRIBUTION FROM BOLTZMANN PIECE !! @@ -666,7 +673,7 @@ subroutine get_fluxes_vmulo(g, phi, pflx, vflx, qflx, & use species, only: spec use stella_geometry, only: grho_norm, bmag, btor use stella_geometry, only: drhodpsi - use stella_geometry, only: gds21, gds22 + use stella_geometry, only: grad_x_dot_grad_y, grad_x_dot_grad_x use stella_geometry, only: dgds21dr, dgds22dr use stella_geometry, only: geo_surf use stella_geometry, only: dBdrho, dIdrho @@ -861,9 +868,12 @@ subroutine get_fluxes_vmulo(g, phi, pflx, vflx, qflx, & g1(:, :, iz, it, ivmu) = g1(:, :, iz, it, ivmu) + g0k ! perpendicular component +! g0k = -g(:, :, iz, it, ivmu) * zi * spread(aky, 2, nakx) * vperp2(ia, iz, imu) * geo_surf%rhoc & +! * (gds21(ia, iz) + theta0 * gds22(ia, iz)) * spec(is)%smz & +! / (geo_surf%qinp * geo_surf%shat * bmag(ia, iz)**2) g0k = -g(:, :, iz, it, ivmu) * zi * spread(aky, 2, nakx) * vperp2(ia, iz, imu) * geo_surf%rhoc & - * (gds21(ia, iz) + theta0 * gds22(ia, iz)) * spec(is)%smz & - / (geo_surf%qinp * geo_surf%shat * bmag(ia, iz)**2) + * (grad_x_dot_grad_y(ia, iz) + theta0 * geo_surf%shat * grad_x_dot_grad_x(ia, iz)) * spec(is)%smz & + / (geo_surf%qinp * bmag(ia, iz)**2) call gyro_average_j1(g0k, iz, ivmu, g2(:, :, iz, it, ivmu)) if (radial_variation) then @@ -871,9 +881,14 @@ subroutine get_fluxes_vmulo(g, phi, pflx, vflx, qflx, & * (dgds21dr(ia, iz) + theta0 * dgds22dr(ia, iz)) * aj1x(:, :, iz, ivmu) * spec(is)%smz & / (geo_surf%qinp * geo_surf%shat * bmag(ia, iz)**2) +! g0k = g0k - g(:, :, iz, it, ivmu) * zi * spread(aky, 2, nakx) * vperp2(ia, iz, imu) * geo_surf%rhoc & +! * (gds21(ia, iz) + theta0 * gds22(ia, iz)) * spec(is)%smz & +! / (geo_surf%qinp * geo_surf%shat * bmag(ia, iz)**2) & +! * (0.5 * aj0x(:, :, iz, ivmu) - aj1x(:, :, iz, ivmu)) & +! * (dkperp2dr(:, :, ia, iz) - dBdrho(iz) / bmag(ia, iz)) g0k = g0k - g(:, :, iz, it, ivmu) * zi * spread(aky, 2, nakx) * vperp2(ia, iz, imu) * geo_surf%rhoc & - * (gds21(ia, iz) + theta0 * gds22(ia, iz)) * spec(is)%smz & - / (geo_surf%qinp * geo_surf%shat * bmag(ia, iz)**2) & + * (grad_x_dot_grad_y(ia, iz) + theta0 * geo_surf%shat * grad_x_dot_grad_x(ia, iz)) * spec(is)%smz & + / (geo_surf%qinp * bmag(ia, iz)**2) & * (0.5 * aj0x(:, :, iz, ivmu) - aj1x(:, :, iz, ivmu)) & * (dkperp2dr(:, :, ia, iz) - dBdrho(iz) / bmag(ia, iz)) @@ -887,10 +902,14 @@ subroutine get_fluxes_vmulo(g, phi, pflx, vflx, qflx, & end if !subtract adiabatic contribution part of g +! g0k = -spec(is)%zt * fphi * phi(:, :, iz, it) * aj0x(:, :, iz, ivmu) * aj1x(:, :, iz, ivmu) & +! * zi * spread(aky, 2, nakx) * vperp2(ia, iz, imu) * geo_surf%rhoc & +! * (gds21(ia, iz) + theta0 * gds22(ia, iz)) * spec(is)%smz & +! / (geo_surf%qinp * geo_surf%shat * bmag(ia, iz)**2) g0k = -spec(is)%zt * fphi * phi(:, :, iz, it) * aj0x(:, :, iz, ivmu) * aj1x(:, :, iz, ivmu) & * zi * spread(aky, 2, nakx) * vperp2(ia, iz, imu) * geo_surf%rhoc & - * (gds21(ia, iz) + theta0 * gds22(ia, iz)) * spec(is)%smz & - / (geo_surf%qinp * geo_surf%shat * bmag(ia, iz)**2) + * (grad_x_dot_grad_y(ia, iz) + theta0 * geo_surf%shat * grad_x_dot_grad_x(ia, iz)) * spec(is)%smz & + / (geo_surf%qinp * bmag(ia, iz)**2) if (.not. maxwellian_normalization) then g0k = g0k * maxwell_vpa(iv, is) * maxwell_mu(ia, iz, imu, is) * maxwell_fac(is) end if @@ -902,11 +921,18 @@ subroutine get_fluxes_vmulo(g, phi, pflx, vflx, qflx, & * (dgds21dr(ia, iz) + theta0 * dgds22dr(ia, iz)) * spec(is)%smz & / (geo_surf%qinp * geo_surf%shat * bmag(ia, iz)**2) + ! g1k = g1k - spec(is)%zt * fphi * phi(:, :, iz, it) * aj0x(:, :, iz, ivmu) & + ! * maxwell_vpa(iv, is) * maxwell_mu(ia, iz, imu, is) * maxwell_fac(is) & + ! * zi * spread(aky, 2, nakx) * vperp2(ia, iz, imu) * geo_surf%rhoc & + ! * (gds21(ia, iz) + theta0 * gds22(ia, iz)) * spec(is)%smz & + ! / (geo_surf%qinp * geo_surf%shat * bmag(ia, iz)**2) & + ! * (0.5 * aj0x(:, :, iz, ivmu) - aj1x(:, :, iz, ivmu)) & + ! * (dkperp2dr(:, :, ia, iz) - dBdrho(iz) / bmag(ia, iz)) g1k = g1k - spec(is)%zt * fphi * phi(:, :, iz, it) * aj0x(:, :, iz, ivmu) & * maxwell_vpa(iv, is) * maxwell_mu(ia, iz, imu, is) * maxwell_fac(is) & * zi * spread(aky, 2, nakx) * vperp2(ia, iz, imu) * geo_surf%rhoc & - * (gds21(ia, iz) + theta0 * gds22(ia, iz)) * spec(is)%smz & - / (geo_surf%qinp * geo_surf%shat * bmag(ia, iz)**2) & + * (grad_x_dot_grad_y(ia, iz) + theta0 * geo_surf%shat * grad_x_dot_grad_x(ia, iz)) * spec(is)%smz & + / (geo_surf%qinp * bmag(ia, iz)**2) & * (0.5 * aj0x(:, :, iz, ivmu) - aj1x(:, :, iz, ivmu)) & * (dkperp2dr(:, :, ia, iz) - dBdrho(iz) / bmag(ia, iz)) diff --git a/stella_io.fpp b/stella_io.fpp index a3fa26dffc..0fa1677a27 100644 --- a/stella_io.fpp +++ b/stella_io.fpp @@ -586,7 +586,7 @@ contains # ifdef NETCDF use neasyf, only: neasyf_write use stella_geometry, only: bmag, gradpar, gbdrift, gbdrift0, & - cvdrift, cvdrift0, gds2, gds21, gds22, grho, jacob, & + cvdrift, cvdrift0, grad_y_dot_grad_y, grad_x_dot_grad_y, grad_x_dot_grad_x, grho, jacob, & drhodpsi, djacdrho, b_dot_grad_z use stella_geometry, only: geo_surf use physics_parameters, only: beta @@ -611,9 +611,9 @@ contains call neasyf_write(file_id, "cvdrift", cvdrift, dim_names=flux_surface_dim) call neasyf_write(file_id, "cvdrift0", cvdrift0, dim_names=flux_surface_dim) call neasyf_write(file_id, "kperp2", kperp2, dim_names=[character(len=5)::"ky", "kx", "alpha", "zed"]) - call neasyf_write(file_id, "gds2", gds2, dim_names=flux_surface_dim) - call neasyf_write(file_id, "gds21", gds21, dim_names=flux_surface_dim) - call neasyf_write(file_id, "gds22", gds22, dim_names=flux_surface_dim) + call neasyf_write(file_id, "grad_y_dot_grad_y", grad_y_dot_grad_y, dim_names=flux_surface_dim) + call neasyf_write(file_id, "grad_x_dot_grad_y", grad_x_dot_grad_y, dim_names=flux_surface_dim) + call neasyf_write(file_id, "grad_x_dot_grad_x", grad_x_dot_grad_x, dim_names=flux_surface_dim) call neasyf_write(file_id, "grho", grho, dim_names=flux_surface_dim) call neasyf_write(file_id, "jacob", jacob, dim_names=flux_surface_dim) call neasyf_write(file_id, "djacdrho", djacdrho, dim_names=flux_surface_dim)