gkvp_f0.63

New feature: velocity-space coordinate (vl,vp)

Previous GKV code employs (vl,mu) velocity-space coordinates, namely, the parallel velocity and magnetic moment. The gkvp_f0.63 added a new choice of velocity-space coordinate (vl,vp), replacing the perpendicular velocity-space coordinate from the magnetic moment to the perpendicular velocity.

This choice is beneficial when the magnetic field strength significantly varies, like the dipole geometry.

Switching the velocity-space coordinates is done in run/gkvp_namelist:

    vp_coord = 0, # For (vl,mu) coordinates

or

    vp_coord = 1, # For (vl,vp) coordinates

In the equation, the parallel advection term is modified:

$$\frac{\partial}{\partial z}_{(z,v_\parallel,\mu)} \rightarrow \frac{\partial}{\partial z}_{(z,v_\parallel,v_\perp)} + \frac{v_\perp}{2B}\frac{\partial B}{\partial z}\frac{\partial}{\partial v_\perp}_{(z,v_\parallel,v_\perp)}$$

Note that the recent rotating flux-tube implementation, the rotating operator is unchanged,

$$\frac{\partial}{\partial t''}_{(z,v_\parallel,\mu)} - \frac{\gamma_E}{\hat{s}} \frac{\partial}{\partial z''}_{(z,v_\parallel,\mu)} \rightarrow \frac{\partial}{\partial t''}_{(z,v_\parallel,v_\perp)} - \frac{\gamma_E}{\hat{s}} \frac{\partial}{\partial z''}_{(z,v_\parallel,v_\perp)}$$

because the operators of $\partial_{v_{\perp}}$ from the first ($\partial_{t''}$) and second (\partial_{z''}) terms cancel out.

What's Changed

• Choice of (vl,vp) velocity-space coordinates by @smaeyama in #14

Full Changelog: f0.62...f0.63

gkvp_f0.62

New feature: Ring dipole geometry

equib_type = "ring" is added for the ring dipole geometry.

In run/gkvp_namelist,

&ring  ring_a = 0.5d0,
       kxmin  = 0.05d0, &end

specify the field line.

Consider flux coordinates (Psi,Theta,phi), where the magnetic poloidal flux Psi<0, the geometrical poloidal angle Theta = arctan(Z/(R-a)), the azimuthal angle of the cylindrical coordinate phi. There is a ring current in direction of phi at R=a. The field line passing through (R,Z)=(R0,0) is picked up as a flux-tube domain. GKV coordinates (x,y,z) are (right-handed system):
x = cx*(Psi0 - Psi)/Psi0
y = cy*phi
z = Theta
with cx=R0, cy=R0. Note that Psi0 is the magnetic poloidal flux at the center of the considered flux-tube domain. In these definitions, the factor on the magnetic field B = cb * \nabla x \times \nabla y is cb = Psi0/(R0*R0) = B0, where B0 is the magnetic field strength at (R,Z)=(R0,0). Normalized omg = B(z)/B0 and cb = 1 in the B0 unit. The reference length is set to be R0 (not the ring current at R=a). The normalized parameter to specify the flux-tube is
ring_a = a / R0

New feature: Rotating flux tube model for background flow shear

In gkvp_f0.62, the novel formulation for treating background flow shear is implemented and named the rotating flux-tube model. It is presented in ICNSP 2022 (Abstract_ICNSP_maeyama_v01.pdf). The full paper is under preparation.

Until gkvp_f0.61, the effect of background flow shear (gamma_e in run/gkvp_namelist) is implemented by the so-called wave-vector remap method with nearest grid-point approximation [Hammett, APS (2006)], which is discontinuous in time. The rotating flux-tube model in gkvp_f0.62 has lots of advantages:

• Time-continuous implementation.
• A linear mode is described by a single (kx,ky) in ballooning space.
• Mathematical consistency with Floquet theorem for periodic time-dependent operator.
• Nonlinear computation cost is reduced than gkvp_f0.61 which requires large box size in radial direction x to reduce the error in time-discontinuous remap for small gamma_e. Properly resolved wave-number space with keeping ballooning angle constant in gkvp_f0.61 allows the same box size and resolution even for small gamma_e.
• Since the remap procedure is not required in the rotating flux-tube model, the analysis of the time evolution of the mode (e.g., nonlinear mode-to-mode interactions) is available even in the presence of shear flow. This is a distinguished feature from another time-continuous implementation with remap [e.g. by Christen, JPP (2021)].

Followings are source-code-level explanations for developer.

• The new module GKV_geom (gkvp_geom.f90) is added. This module extract geometry settings from gkvp_set.f90.
• GKV_geom also stores metric information in -global_nz:global_nz-1 and Fourier coefficients of omg, domgdr, domgdt, domgdq, ggrr, ggrt, ..., defiend in the flux coordinates (r,t,q). Using the Fourier coefficients, metric at arbitrary position zz is calculated. Relation between flux coordinates (r,t,q) and field aligned coordinates (x,y,z) defines the metric domgdx, domgdy, domgdz, ggxx, ggxy, ..., which are used to calculate various operators, like vdx, vdy, ksq, fmx, j0, ....
• The argument of "iz" of the subroutine vmecbzx_boozx_coeff in gkvp_vmecbzx.f90 is changed from local iz for each rankz to global giz.
• All other modifications in source code is small; only in gkvp_advnc.f90, gkvp_header.f90, gkvp_main.f90. These are marked as !%%% For shearflow rotating flux tube model %%%.

gkvp_f0.61

• equib_type="s-alpha-shift" is added.
  Previous equib_type="s-alpha" is large-aspect tokamak model without Shafranov shift, which is valid in low beta limit.
  New equib_type="s-alpha-shift" includes Shafranov shift, that is used in A. Ishizawa, et al., PRL 123, 025003 (2019).

• random_seed is added in src/gkvp_math_portable.f90 for reproducibility of initial random phases (, but not ensured between different architectures).

• Modify initial random phase rr in src/gkvp_set.f90, so as to be specified by global indices of (kx,ky). Then, the initial random phase is independent of wave-number-space domain decomposition by MPI.

• Special treatment of collision for tracer particles (fcs=0).
  Set the field-particle collision operator to be zero for the tracer particles, namely, gy_fld(ia,:)=0 in src/gkvp_colliimp.f90. Physically, bulk plasma species (electrons and main ions) collide only with bulk plasma species. Tracer particles (fcs=0) only has test-particle collision with bulk plasma.
  Setting of fcs=0 nullifies most of output about the tracer particle (e.g., particle flux hst/gkvp.ges.X.001 or fluid momentum data phi/gkvp*.mom.001 become zero) in the present GKV normalization. Only availables are cnt/ and fxv/.

gkvp_f0.60

NetCDF4+Parallel HDF5 is added for GKV binary output. This is the parallel file I/O of NetCDF4 interface with Parallel HDF5 engine(, and is NOT parallel-NetCDF (pnetcdf).)

src/gkvp_fileio_fortran.f90 is a module for Fortran binary output used so far,
src/gkvp_fileio_netcdf.f90 is a new module for NetCDF4+Parallel HDF5.

In Makefile, these modules are switch by choosing appropriate source file during compilation,
FILEIO=gkvp_fileio_fortran
or
FILEIO=gkvp_fileio_netcdf

For now, verification tests of NetCDF4+parallel HDF5 are finished on

Ubuntu (WSL2) on laptop (FFTW and NetCDF by SPACK)
p-srv at P Lab. (build by ifort, hdf5-1.10.7, netcdf-c-4.7.4, netcdf-fortran-4.5.3)
JFRS-1 at QST (Cray, XC50. module load cray-fftw cray-hdf5-parallel/1.10.5.1 cray-netcdf-hdf5parallel/4.6.3.1)
Plasma Simulator at NIFS (NEC, SX-Aurora TSUBASA. module load NECNLC-sx netcdf-parallelIO-fortran-sx/4.5.2)
ITO at Kyushu Univ. (Fujitsu PRIMERGY CX2550/CX2560 M4. module load netcdf)

There are still problem of NetCDF4+parallel HDF5 on some systems (I am now asking to helpdesk):

FLOW at Nagoya Univ. (FUJITSU PRIMEHPC FX1000. Parallel HDF5 does not work?)
Fugaku at RIKEN (FUJITSU PRIMEHPC FX1000. Parallel HDF5 does not work?)

gkvp_f0.59

Updates from gkvp_f0.55 (available from GKV homepage) are summarized below.

• Version number f0.** is removed from filename.
  *** Update history should be clarified in each module files too. ***

• Makefile, shoot, sub.q for Flow supercomputer (Nagoya Univ, 2020) are added.

• Makefile, shoot, sub.q for Fugaku supercomputer (RIKEN R-CCS, 2020) are added.

• Makefile, shoot, sub.q for SX-Aurora TSUBASA (NIFS, 2020) are added. Optimized source codes by NEC and NIFS are enclosed as src/gkvp_f0.56_tune_***.f90, which is switched by Makefile.

• equib_type = "slab" is added to gkvp_namelist for shearless slab geometry.

• "init_random" is added to gkvp_namelist to turn on/off initial random numbers.

• benchmarks/ contains some benchmark data of GKV (ITG with adiabatic/kinetic electrons). One may reproduce the standard outputs in hst/ by setting gkvp_namelist and gkvp_header.f90.

• A bug at bias factor of collision, nu, is modified in colliimp module.

• (s-alpha model with Shafranov shift [Ishizawa, PRL (2019)] is supported, but it is still in source code level. One needs to uncomment lines with "with Shafranov shift" in gkvp_set.f90.)

• Minor debugs:
  Padding region iend_y < my is carefully initalized.
  menergy(:,:)=0 for beta=0
  freq module is modified to support nx0=0.

• Extensions of diag.
  out_mominxml.f90 has been maintained to draw 3D full torus visualization by Paraview.
  out_netcdf.f90 is added to convert GKV binary files, phi, Al, mom, fxv, cnt, trn, tri to NetCDF files.

• Python version of post-processing program diag_python is added. For now, diag_python can read NetCDF files and plot some types of time slices.