Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Adds zed_scalefac to the output netcdf file. #109

Open
wants to merge 3 commits into
base: master
Choose a base branch
from
Open

Adds zed_scalefac to the output netcdf file. #109

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
a309e08
Added zed_scalefac to the output netcdf file. This lets us calculate …
May 10, 2023
0a530c5
zed_scalefac now set properly for miller.
May 11, 2023
3f5c36d
updated
Jun 28, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions flow_shear.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -94,7 +94,7 @@ subroutine init_flow_shear
prl_shear(ia, iz, ivmu) = -omprimfac * g_exb * code_dt * vpa(iv) * spec(is)%stm_psi0 &
* dydalpha * drhodpsi &
* (geo_surf%qinp_psi0 / geo_surf%rhoc_psi0) &
* (btor(iz) * rmajor(iz) / bmag(ia, iz)) * (spec(is)%mass / spec(is)%temp)
* (btor(ia, iz) * rmajor(iz) / bmag(ia, iz)) * (spec(is)%mass / spec(is)%temp)
end do
if (.not. maxwellian_normalization) then
do iz = -nzgrid, nzgrid
Expand All @@ -104,7 +104,7 @@ subroutine init_flow_shear
end if
if (radial_variation) then
energy = (vpa(iv)**2 + vperp2(:, :, imu)) * (spec(is)%temp_psi0 / spec(is)%temp)
prl_shear_p(:, :, ivmu) = prl_shear(:, :, ivmu) * (dIdrho / spread(rmajor * btor, 1, nalpha) &
prl_shear_p(:, :, ivmu) = prl_shear(:, :, ivmu) * (dIdrho / (spread(rmajor, 1, nalpha) * btor) &
- spread(dBdrho, 1, nalpha) / bmag &
- spec(is)%fprim - spec(is)%tprim * (energy - 2.5) &
- 2.*mu(imu) * spread(dBdrho, 1, nalpha))
Expand Down
38 changes: 25 additions & 13 deletions geo/stella_geometry.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,8 @@ module stella_geometry
public :: theta_vmec
public :: zeta
public :: zed_scalefac
public :: vmec_file
public :: geo_file
public :: dxdXcoord, dydalpha
public :: aref, bref
public :: twist_and_shift_geo_fac
Expand Down Expand Up @@ -67,7 +69,8 @@ module stella_geometry
real, dimension(:, :, :), allocatable :: dVolume
real, dimension(:, :), allocatable :: x_displacement_fac
real, dimension(:), allocatable :: dBdrho, d2Bdrdth, dgradpardrho
real, dimension(:), allocatable :: btor, Rmajor
real, dimension(:, :), allocatable :: btor
real, dimension(:), allocatable :: Rmajor
real, dimension(:), allocatable :: alpha
real, dimension(:, :), allocatable :: zeta

Expand All @@ -84,6 +87,7 @@ module stella_geometry
logical :: overwrite_gbdrift, overwrite_cvdrift, overwrite_gbdrift0
logical :: q_as_x
character(100) :: geo_file
character(200) :: vmec_file

logical :: vmec_chosen = .false.
logical :: geoinit = .false.
Expand Down Expand Up @@ -140,7 +144,9 @@ subroutine init_geometry(nalpha, naky)

! default is no re-scaling of zed
zed_scalefac = 1.0

! default is no VMEC file
vmec_file = ""

if (proc0) then
call read_parameters
select case (geo_option_switch)
Expand All @@ -157,7 +163,7 @@ subroutine init_geometry(nalpha, naky)
gds2(1, :), gds21(1, :), gds22(1, :), &
gds23(1, :), gds24(1, :), gradpar, &
gbdrift0(1, :), gbdrift(1, :), cvdrift0(1, :), cvdrift(1, :), &
dBdrho, d2Bdrdth, dgradpardrho, btor, rmajor, &
dBdrho, d2Bdrdth, dgradpardrho, btor(1, :), rmajor, &
dcvdrift0drho(1, :), dcvdriftdrho(1, :), &
dgbdrift0drho(1, :), dgbdriftdrho(1, :), &
dgds2dr(1, :), dgds21dr(1, :), &
Expand Down Expand Up @@ -193,6 +199,8 @@ 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
zed_scalefac = 1/geo_surf%qinp

case (geo_option_multibox)
! read in Miller local parameters
call read_local_parameters(nzed, nzgrid, geo_surf)
Expand All @@ -210,7 +218,7 @@ subroutine init_geometry(nalpha, naky)
gds2(1, :), gds21(1, :), gds22(1, :), &
gds23(1, :), gds24(1, :), gradpar, &
gbdrift0(1, :), gbdrift(1, :), cvdrift0(1, :), cvdrift(1, :), &
dBdrho, d2Bdrdth, dgradpardrho, btor, rmajor, &
dBdrho, d2Bdrdth, dgradpardrho, btor(1, :), rmajor, &
dcvdrift0drho(1, :), dcvdriftdrho(1, :), &
dgbdrift0drho(1, :), dgbdriftdrho(1, :), &
dgds2dr(1, :), dgds21dr(1, :), &
Expand Down Expand Up @@ -246,6 +254,8 @@ 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
zed_scalefac = 1/geo_surf%qinp


case (geo_option_inputprof)
! first read in some local parameters
Expand All @@ -263,7 +273,7 @@ subroutine init_geometry(nalpha, naky)
gds2(1, :), gds21(1, :), gds22(1, :), &
gds23(1, :), gds24(1, :), gradpar, &
gbdrift0(1, :), gbdrift(1, :), cvdrift0(1, :), cvdrift(1, :), &
dBdrho, d2Bdrdth, dgradpardrho, btor, rmajor, &
dBdrho, d2Bdrdth, dgradpardrho, btor(1, :), rmajor, &
dcvdrift0drho(1, :), dcvdriftdrho(1, :), &
dgbdrift0drho(1, :), dgbdriftdrho(1, :), &
dgds2dr(1, :), dgds21dr(1, :), &
Expand All @@ -290,8 +300,9 @@ subroutine init_geometry(nalpha, naky)
twist_and_shift_geo_fac = 2.0 * pi * geo_surf%shat
! aref and bref should not be needed so set to 1
aref = 1.0; bref = 1.0

zeta(1, :) = zed * geo_surf%qinp
zed_scalefac = 1/geo_surf%qinp


case (geo_option_vmec)
vmec_chosen = .true.
Expand Down Expand Up @@ -320,7 +331,7 @@ subroutine init_geometry(nalpha, naky)
gds23, gds24, gds25, gds26, gbdrift_alpha, gbdrift0_psi, &
cvdrift_alpha, cvdrift0_psi, sign_torflux, &
theta_vmec, zed_scalefac, aref, bref, alpha, zeta, &
field_period_ratio, x_displacement_fac)
field_period_ratio, x_displacement_fac, vmec_file)

write (*, '(A)') "############################################################"
write (*, '(A)') " BOUNDARY CONDITIONS"
Expand All @@ -344,7 +355,7 @@ subroutine init_geometry(nalpha, naky)
gds23, gds24, gds25, gds26, gbdrift_alpha, gbdrift0_psi, &
cvdrift_alpha, cvdrift0_psi, sign_torflux, &
theta_vmec, zed_scalefac, aref, bref, alpha, zeta, &
field_period_ratio, x_displacement_fac)
field_period_ratio, x_displacement_fac, vmec_file)
! Restart the variable twist_and_shift_geo_fac_full
twist_and_shift_geo_fac_full = 0
end if
Expand Down Expand Up @@ -599,8 +610,7 @@ subroutine overwrite_selected_geometric_coefficients(nalpha)
gds2_file, gds21_file, gds22_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_gradpar) b_dot_grad_z(ia, iz) = gradpar_file
if (overwrite_gds2) gds2(ia, iz) = gds2_file
if (overwrite_gds21) gds21(ia, iz) = gds21_file
if (overwrite_gds22) gds22(ia, iz) = gds22_file
Expand All @@ -612,6 +622,8 @@ subroutine overwrite_selected_geometric_coefficients(nalpha)
end do
end do
cvdrift0 = gbdrift0
if (overwrite_gradpar) gradpar = sum(b_dot_grad_z, dim=1)


close (geofile_unit)

Expand Down Expand Up @@ -697,10 +709,10 @@ subroutine allocate_arrays(nalpha, nzgrid)
if (.not. allocated(dl_over_b)) allocate (dl_over_b(nalpha, -nzgrid:nzgrid))
if (.not. allocated(d_dl_over_b_drho)) allocate (d_dl_over_b_drho(nalpha, -nzgrid:nzgrid))
if (.not. allocated(b_dot_grad_z)) allocate (b_dot_grad_z(nalpha, -nzgrid:nzgrid))

if (.not. allocated(btor)) allocate (btor(nalpha, -nzgrid:nzgrid))

if (.not. allocated(gradpar)) allocate (gradpar(-nzgrid: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))
Expand Down Expand Up @@ -1015,7 +1027,7 @@ subroutine write_geometric_coefficients(nalpha)
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), &
gds24(ia, iz), gbdrift(ia, iz), cvdrift(ia, iz), gbdrift0(ia, iz), &
bmag_psi0(ia, iz), btor(iz)
bmag_psi0(ia, iz), btor(ia,iz)
end do
write (geometry_unit, *)
end do
Expand Down
7 changes: 5 additions & 2 deletions geo/vmec_geo.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@ module vmec_geo
real :: nfield_periods
real :: zeta_center, torflux
logical :: verbose
character(2000) :: vmec_filename
character(200) :: vmec_filename
integer :: n_tolerated_test_arrays_inconsistencies

contains
Expand Down Expand Up @@ -97,7 +97,7 @@ subroutine get_vmec_geo(new_zeta_min, stellarator_symmetric_BC, &
gds23, gds24, gds25, gds26, gbdrift_alpha, gbdrift0_psi, cvdrift_alpha, &
cvdrift0_psi, sign_torflux, &
theta_vmec, zed_scalefac, L_reference, B_reference, alpha, zeta, &
field_period_ratio, x_displacement_fac)
field_period_ratio, x_displacement_fac, vmec_file)

use constants, only: pi
use common_types, only: flux_surface_type
Expand Down Expand Up @@ -127,6 +127,8 @@ subroutine get_vmec_geo(new_zeta_min, stellarator_symmetric_BC, &
integer, intent(out) :: sign_torflux
real, intent(out) :: field_period_ratio

character(200), intent(out) :: vmec_file

logical, parameter :: debug = .false.

integer :: tmpunit
Expand Down Expand Up @@ -156,6 +158,7 @@ subroutine get_vmec_geo(new_zeta_min, stellarator_symmetric_BC, &

integer :: ierr

vmec_file = vmec_filename
!> To avoid writting twice in the output file when recomputing zeta.
if (stellarator_symmetric_BC) verbose = .false.
!> first read in equilibrium information from vmec file
Expand Down
12 changes: 6 additions & 6 deletions stella_diagnostics.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -576,7 +576,7 @@ subroutine get_fluxes(g, pflx, vflx, qflx, &

! get momentum flux
! parallel component
gtmp1 = g(:, :, ikxkyz) * spread(vpa, 2, nmu) * geo_surf%rmaj * btor(iz) / bmag(ia, iz)
gtmp1 = g(:, :, ikxkyz) * spread(vpa, 2, nmu) * geo_surf%rmaj * btor(ia, 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 &
Expand Down Expand Up @@ -610,7 +610,7 @@ subroutine get_fluxes(g, pflx, vflx, qflx, &
! Apar contribution to momentum flux
! parallel component
gtmp1 = -spread(vpa**2, 2, nmu) * spec(is)%stm * g(:, :, ikxkyz) &
* geo_surf%rmaj * btor(iz) / bmag(1, iz)
* geo_surf%rmaj * btor(ia, iz) / bmag(ia, iz)
call gyro_average(gtmp1, ikxkyz, gtmp2)
! perp component
gtmp1 = spread(vpa, 2, nmu) * spec(is)%stm * g(:, :, ikxkyz) &
Expand Down Expand Up @@ -825,15 +825,15 @@ subroutine get_fluxes_vmulo(g, phi, pflx, vflx, qflx, &
do it = 1, ntubes
do iz = -nzgrid, nzgrid
! parallel component
g0k = g(:, :, iz, it, ivmu) * vpa(iv) * geo_surf%rmaj * btor(iz) / bmag(ia, iz)
g0k = g(:, :, iz, it, ivmu) * vpa(iv) * geo_surf%rmaj * btor(ia, iz) / bmag(ia, iz)
call gyro_average(g0k, iz, ivmu, g1(:, :, iz, it, ivmu))

if (radial_variation) then
g0k = g1(:, :, iz, it, ivmu) &
* (-0.5 * aj1x(:, :, iz, ivmu) / aj0x(:, :, iz, ivmu) * (spec(is)%smz)**2 &
* (kperp2(:, :, ia, iz) * vperp2(ia, iz, imu) / bmag(ia, iz)**2) &
* (dkperp2dr(:, :, ia, iz) - dBdrho(iz) / bmag(ia, iz)) &
+ dIdrho / (geo_surf%rmaj * btor(iz)))
+ dIdrho / (geo_surf%rmaj * btor(ia,iz)))

call multiply_by_rho(g0k)

Expand All @@ -842,7 +842,7 @@ 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)**2 &
* vpa(iv) * geo_surf%rmaj * btor(iz) / bmag(ia, iz)
* vpa(iv) * geo_surf%rmaj * btor(ia, iz) / bmag(ia, iz)
if (.not. maxwellian_normalization) then
g0k = g0k * maxwell_vpa(iv, is) * maxwell_mu(ia, iz, imu, is) * maxwell_fac(is)
end if
Expand All @@ -852,7 +852,7 @@ subroutine get_fluxes_vmulo(g, phi, pflx, vflx, qflx, &
- aj1x(:, :, iz, ivmu) / aj0x(:, :, iz, ivmu) * (spec(is)%smz)**2 &
* (kperp2(:, :, ia, iz) * vperp2(ia, iz, imu) / bmag(ia, iz)**2) &
* (dkperp2dr(:, :, ia, iz) - dBdrho(iz) / bmag(ia, iz)) &
+ dIdrho / (geo_surf%rmaj * btor(iz)))
+ dIdrho / (geo_surf%rmaj * btor(ia, iz)))

call multiply_by_rho(g1k)

Expand Down
62 changes: 59 additions & 3 deletions stella_io.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,8 @@ module stella_io
# ifdef NETCDF
integer(kind_nf) :: ncid
integer(kind_nf) :: char10_dim
integer(kind_nf) :: char200_dim


integer :: code_id

Expand Down Expand Up @@ -122,7 +124,7 @@ contains

! Dimensions for various string variables
call neasyf_dim(file_id, "char10", dim_size=10, dimid=char10_dim)
call neasyf_dim(file_id, "char200", dim_size=200)
call neasyf_dim(file_id, "char200", dim_size=200, dimid=char200_dim)
call neasyf_dim(file_id, "nlines", unlimited=.true., long_name="Input file line number")

call neasyf_dim(file_id, "ri", dim_size=2, long_name="Complex components", units="(real, imaginary)")
Expand Down Expand Up @@ -193,7 +195,7 @@ contains
integer, parameter :: line_length = 200
character(line_length), dimension(:), allocatable :: input_file_array
integer :: n, unit, status, dim_id, previous_nlines

! Don't attempt to write zero-sized arrays
if (num_input_lines <= 0) return

Expand Down Expand Up @@ -587,16 +589,23 @@ contains
use neasyf, only: neasyf_write
use stella_geometry, only: bmag, gradpar, gbdrift, gbdrift0, &
cvdrift, cvdrift0, gds2, gds21, gds22, grho, jacob, &
drhodpsi, djacdrho, b_dot_grad_z
drhodpsi, djacdrho, b_dot_grad_z, zed_scalefac, alpha, &
vmec_file, geo_file, gds23, gds24, btor, bmag_psi0, &
aref, bref, exb_nonlin_fac, exb_nonlin_fac_p
use stella_geometry, only: geo_surf
use physics_parameters, only: beta
use dist_fn_arrays, only: kperp2
use kt_grids, only: naky, nakx, jtwist

use netcdf!, only: nf90_inq_varid, nf90_def_var, nf90_enddef, nf90_put_var, NF90_NOERR, NF90_EBADDIM
#endif
implicit none
!> NetCDF ID of the file to write to
integer, intent(in) :: file_id

integer :: status, code_id


# ifdef NETCDF
character(len=*), dimension(*), parameter :: flux_surface_dim = ["alpha", "zed "]

Expand All @@ -614,9 +623,19 @@ contains
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, "gds23", gds23, dim_names=flux_surface_dim)
call neasyf_write(file_id, "gds24", gds24, dim_names=flux_surface_dim)
call neasyf_write(file_id, "btor", btor, dim_names=flux_surface_dim)
call neasyf_write(file_id, "bmag_psi0", bmag_psi0, 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)
call neasyf_write(file_id, "zed_scalefac", zed_scalefac, long_name="Factor to divide zed with to get zeta.")
!call neasyf_write(file_id, "vmec_file", vmec_file, long_name="VMEC filename used to calculate geometry input.", dim_names=["char200"])
!call neasyf_write(file_id, "geo_file", geo_file, long_name="Geometry file potentially used for overriding geometric coefficients.", dim_names=["char200"])

call neasyf_write(file_id, "alpha_grid", alpha, dim_names=[character(len=5)::"alpha"], long_name="Field line labels in simulation.")
call neasyf_write(file_id, "beta", beta, &
long_name="Reference beta", units="2.mu0.nref.Tref/B_a**2")
call neasyf_write(file_id, "q", geo_surf%qinp, &
Expand All @@ -629,6 +648,43 @@ contains
call neasyf_write(file_id, "d2psidr2", geo_surf%d2psidr2)
call neasyf_write(file_id, "jtwist", jtwist, &
long_name="2*pi*shat*dky/dkx")

call neasyf_write(file_id, "aref", aref, &
long_name="Reference length")
call neasyf_write(file_id, "bref", bref, &
long_name="Reference magnetic field")
call neasyf_write(file_id, "rhotor", geo_surf%rhotor, &
long_name="Radial location in sqrt(normalized toroidal flux)")
call neasyf_write(file_id, "rhoc", geo_surf%rhoc, &
long_name="Radial location in sqrt(normalized toroidal flux)")
call neasyf_write(file_id, "exb_nonlin", exb_nonlin_fac, &
long_name="Radial location in sqrt(normalized toroidal flux)")
call neasyf_write(file_id, "exb_nonlin_p", exb_nonlin_fac_p, &
long_name="Radial location in sqrt(normalized toroidal flux)")


! Could not get writing string to work with neasy-f, so using the normal interface.
status = nf90_inq_varid(file_id, 'vmec_file', code_id)
if (status /= nf90_noerr) then
status = nf90_def_var(file_id, 'vmec_file', nf90_char, char200_dim, code_id)
if (status /= nf90_noerr) call netcdf_error(status, file_id, code_id, att="vmec_file")
status = nf90_enddef(file_id)
if (status /= nf90_noerr) call netcdf_error(status, file_id, code_id, att="vmec_file")
end if
status = nf90_put_var(file_id, code_id, trim(vmec_file))
if(status /= nf90_noerr) call netcdf_error(status, file_id, code_id, att="vmec_file")

status = nf90_inq_varid(file_id, 'geo_file', code_id)
if (status /= nf90_noerr) then
status = nf90_def_var(file_id, 'geo_file', nf90_char, char200_dim, code_id)
if (status /= nf90_noerr) call netcdf_error(status, file_id, code_id, att="geo_file")
status = nf90_enddef(file_id)
if (status /= nf90_noerr) call netcdf_error(status, file_id, code_id, att="geo_file")
end if
status = nf90_put_var(file_id, code_id, trim(geo_file))
if(status /= nf90_noerr) call netcdf_error(status, file_id, code_id, att="geo_file")


# endif
end subroutine nc_geo

Expand Down