stellaGK
Clean-up stella's input file / namelists / input variables.
@GeorgiaActon and I have written out our ideal stella input file, so that we can clean up the stella namelists and variables. The goal is that the namelists will be more intuitive to understand and easier to use. Before merging this into the master a Python script will be written to convert the old input files to these new input files. We copy our ideal input file here below so other users can make comments / suggestions on the input file while we work on this clean up.
!-------------------------------
!---------- Geometry ----------
!-------------------------------
! Define geometry parameters
! Depending on <geo_option>, different extra namelists will be read:
! <geo_option> = ['local', 'default', 'miller'] uses the <millergeo_parameters> namelist
! <geo_option> = 'vmec' uses the <vmec_parameters> namelist
! <geo_option> = 'zpinch' does not use any extra namelists
! <geo_option> = 'input.profiles' does not use any extra namelists
! Note that q_as_x = <radial_variation>, from the <parameters_physics> namelist
&geometry_option
geo_option = "local"
q_as_x = .false.
/
&overwrite_geometry
geo_file = "input.geometry"
overwrite_bmag = .false.
overwrite_b_dot_grad_zeta = .false.
overwrite_gds2 = .false.
overwrite_gds21 = .false.
overwrite_gds22 = .false.
overwrite_gds23 = .false.
overwrite_gds24 = .false.
overwrite_gbdrift = .false.
overwrite_cvdrift = .false.
overwrite_gbdrift0 = .false.
set_bmag_const = .false.
/
&geometry_vmec (renamed from vmec_parameters)
alpha0 = 0.0
zeta_center = 0.0
rectangular_cross_section = .false.
nfield_periods = -1.0
torflux = 0.6354167
zgrid_refinement_factor = 1.0
surface_option = 0.0
radial_coordinate = "sgn(psi_t)psi_t"
verbose = .true.
vmec_filename = "equilibria/wout_w7x_standardconfig.nc"
n_tolerated_test_arrays_inconsistencies = 0.0
/
&geometry_miller (renamed from millergeo_parameters)
rhoc = 0.5
rmaj = 3.0
shift = 0.0
qinp = 1.4
shat = 0.8
kappa = 0.0
kapprim = 0.0
tri = 0.0
triprim = 0.0
rgeo = 3.0
betaprim = 0.0
betadbprim = 0.0
d2qdr2 = 0.0
d2psidr2 = 0.0
nzed_local = 128.0
read_profile_variation = .false.
write_profile_variation = .false.
/
!-------------------------------
!---------- Physics ----------
!-------------------------------
&gyrokinetic_terms
simulation_domain = 'fluxtube' !! options are{ft, ffs, multibox}
include_parallel_streaming = .true.
include_mirror = .true.
include_xdrift = .true. !! xdriftknob = 1.0 or 0.0
include_ydrift = .true. !! ydriftknob = 1.0 or 0.0
include_drive = .true. !! wstarknob = 1.0 or 0.0
include_nonlinear = .false.
include_parallel_nonlinearity = .false.
include_electromagnetic = .false.
/
&scale_gyrokinetic_terms
xdriftknob = 1.0
ydriftknob = 1.0
wstarknob = 1.0
fphi = 1.0
suppress_zonal_interaction = .false.
/
&adiabatic_electron_response
! Abort if kinetic electrons are included
adiabatic_option = "field-line-average-term" ! Remove this flag if adiabatic_ion_response namelist exists
tite = 1.0
nine = 1.0
/
&adiabatic_ion_response
! The <adiabatic_option> flag should not be set by the user
! Stella should be automatically choosing the option without the flux surface
! average for adiabatic ions and with the flux surface average for electrons
/
&electromagnetic
!! set the below to false if include_electromagnetic = .false.
include_apar = .true.
include_bpar = .true.
beta = 0.0
/
&full_flux_surface
! Abort stella if y0 and rhostar are both set
rhostar = -1.0
irhostar = -1.0 !! Delete irhostar
/
&extra_physics
prp_shear_enabled = .false.
hammett_flow_shear = .true.
include_pressure_variation = .false.
g_exb = 0.0
g_exbfac = 1.0
omprimfac = 1.0
/
!-------------------------------
!----------- Grids -----------
!-------------------------------
&vpamu_grid (renamed from vpamu_grids_parameters)
nvgrid = 24.0
nmu = 12.0
vpa_max = 3.0
vperp_max = 3.0
equally_spaced_mu_grid = .false.
conservative_wgts_vpa = .false.
/
&z_grid (renamed from zgrid_parameters)
nzed = 24.0
nperiod = 1.0
ntubes = 1.0
zed_equal_arc = .false.
/
&z_boundary_condition
boundary_option = "default"
shat_zero = 1e-05
grad_x_grad_y_zero = 1e-05
dkx_over_dky = -1.0
/
&species_knobs
nspec = 2.0
species_option = "stella"
read_profile_variation = .false.
write_profile_variation = .false.
/
&species_parameters_1
z = 1.0
mass = 1.0
dens = 1.0
temp = 1.0
tprim = -999.9
fprim = -999.9
d2ndr2 = 0.0
d2tdr2 = 0.0
bess_fac = 1.0
type = "default"
/
&species_parameters_2
z = 1.0
mass = 1.0
dens = 1.0
temp = 1.0
tprim = -999.9
fprim = -999.9
d2ndr2 = 0.0
d2tdr2 = 0.0
bess_fac = 1.0
type = "default"
/
&kxky_grid_option
grid_option = "default"
/
&kxky_grid_range
naky = 1.0
nakx = 1.0
aky_min = 0.0
aky_max = 0.0
theta0_min = 0.0
theta0_max = -1.0
akx_min = 0.0
akx_max = -1.0
kyspacing_option = "default"
/
&kxky_grid_box
nx = 1.0
ny = 1.0
jtwist = -1.0
jtwistfac = 1.0
x0 = -1.0
y0 = -1.0
centered_in_rho = .true.
periodic_variation = .false.
randomize_phase_shift = .false.
phase_shift_angle = 0.0
/
!--------------------------------
!--------- Diagnostics --------
!--------------------------------
&diagnostics
nwrite = 50.0
navg = 50.0
nsave = -1.0
nc_mult = 1.0
save_for_restart = .false.
write_all = .false.
write_all_time_traces = .true.
write_all_spectra_kxkyz = .false.
write_all_spectra_kxky = .false.
write_all_velocity_space = .false.
write_all_potential = .false.
write_all_omega = .false.
write_all_distribution = .false.
write_all_fluxes = .false.
write_all_moments = .false.
/
&diagnostics_potential
write_phi2_vs_time = .true.
write_apar2_vs_time = .true.
write_bpar2_vs_time = .true.
write_phi_vs_kxkyz = .false.
write_phi2_vs_kxky = .false.
/
&diagnostics_omega
write_omega_vs_kxky = .true.
write_omega_avg_vs_kxky = .false.
/
&diagnostics_distribution
write_g2_vs_vpamus = .false.
write_g2_vs_zvpas = .false.
write_g2_vs_zmus = .false.
write_g2_vs_kxkyzs = .false.
write_g2_vs_zvpamus = .false.
write_distribution_g = .true.
write_distribution_h = .false.
write_distribution_f = .false.
/
&diagnostics_fluxes
flux_norm = .true.
write_fluxes_vs_time = .true.
write_radial_fluxes = .false.
write_fluxes_kxkyz = .false.
write_fluxes_kxky = .false.
/
&diagnostics_moments
write_moments = .false.
write_radial_moments = .false.
/
!--------------------------------
!-------- Init fields ---------
!--------------------------------
! Choose the option, then read the namelist of this option
&init_distribution (renamed from init_g_knobs)
ginit_option = "default"
phiinit = 1.0
chop_side = .false.
left = .true.
scale_to_phiinit = .false.
/
&init_distribution_default
width0 = -3.5
den0 = 1.0
upar0 = 0.0
oddparity = .false.
/
&init_distribution_noise
zf_init = 1.0
rng_seed = -1
/
&init_distribution_kpar
width0 = -3.5
den0 = 1.0
upar0 = 0.0
tpar0 = 0.0
tperp0 = 0.0
den1 = 0.0
upar1 = 0.0
tpar1 = 0.0
tperp1 = 0.0
den2 = 0.0
upar2 = 0.0
tpar2 = 0.0
tperp2 = 0.0
/
&init_distribution_rh
kxmax = 1e+100
kxmin = 0.0
imfac = 0.0
refac = 1.0
/
&restart_options
scale = 1.0
tstart = 0.0
restart_file = "dummy.nc"
restart_dir = "./"
read_many = .true.
/
!-------------------------------
!-------- Dissipation --------
!-------------------------------
&dissipation_and_collisions_options (Renamed from &dissipation)
include_collisions = .false.
collisions_implicit = .true.
collision_model = "dougherty"
hyper_dissipation = .false.
ecoll_zeff = .false.
zeff = 1.0
vnew_ref = -1.0
/
&collisions_dougherty
momentum_conservation = .true.
energy_conservation = .true.
vpa_operator = .true.
mu_operator = .true.
/
&collisions_fokker_planck (Renamed from &collisions_fp)
testpart = .true.
fieldpart = .false.
lmax = 1.0
jmax = 1.0
nvel_local = 512.0
interspec = .true.
intraspec = .true.
iiknob = 1.0
ieknob = 1.0
eeknob = 1.0
eiknob = 1.0
eiediffknob = 1.0
eideflknob = 1.0
deflknob = 1.0
eimassr_approx = .false.
advfield_coll = .true.
spitzer_problem = .false.
density_conservation = .false.
density_conservation_field = .false.
density_conservation_tp = .false.
exact_conservation = .false.
exact_conservation_tp = .false.
vpa_operator = .true.
mu_operator = .true.
cfac = 1.0
cfac2 = 1.0
nuxfac = 1.0
i1fac = 1.0
i2fac = 0.0
no_j1l1 = .true.
no_j1l2 = .false.
no_j0l2 = .false.
/
&hyper_dissipation
d_hyper = 0.05
d_zed = 0.05
d_vpa = 0.05
hyp_zed = .false.
hyp_vpa = .false.
use_physical_ksqr = .true.
scale_to_outboard = .false.
/
!-------------------------------
!---------- Numerics ---------
!-------------------------------
&time_trace_options
nstep = -1.0
tend = -1.0
autostop = .true.
avail_cpu_time = 10000000000.0
/
&time_step
delt = 0.03
delt_option = "default"
delt_max = -1.0
delt_min = 1e-10
cfl_cushion_upper = 0.5
cfl_cushion_middle = 0.25
cfl_cushion_lower = 1e-05
/
&numerical_algorithms
explicit_algorithm = "rk3" (renamed from explicit_option)
flip_flop = .false.
stream_implicit = .true.
stream_iterative_implicit = .false.
stream_matrix_inversion = .false.
driftkinetic_implicit = .false.
mirror_implicit = .true.
mirror_semi_lagrange = .true.
mirror_linear_interp = .false.
drifts_implicit = .false.
fully_implicit = .false.
fully_explicit = .false.
maxwellian_inside_zed_derivative = .false.
use_deltaphi_for_response_matrix = .false.
maxwellian_normalization = .false. !!!! DELETE THIS
/
&numerical_upwinding_for_derivatives
zed_upwind = 0.02
vpa_upwind = 0.02
time_upwind = 0.02
/
!-------------------------------
!-------- Neoclassics --------
!-------------------------------
&neoclassical_input
include_neoclassical_terms = .false.
neo_option = "sfincs"
nradii = 5.0
drho = 0.01
/
&euterpe_parameters
nradii = 1000.0
data_file = "euterpe.dat"
/
&sources
source_option = "none"
nu_krook = 0.05
tcorr_source = 0.02
ikxmax_source = 1.0
krook_odd = .true.
exclude_boundary_regions = .false.
tcorr_source_qn = 0.0
exclude_boundary_regions_qn = .false.
from_zero = .true.
conserve_momentum = .false.
conserve_density = .false.
/
!--------------------------------
!------ Radial variation ------
!--------------------------------
&radial_variation !! merged with &multibox_parameters
radial_variation = .false.
include_geometric_variation = .true.
ky_solve_radial = 0.0
ky_solve_real = .false.
boundary_size = 4.0
krook_size = 0.0
smooth_zfs = .false.
zf_option = "default"
lr_debug_option = "default"
krook_option = "default"
rk_step = .false.
nu_krook_mb = 0.0
mb_debug_step = -1.0
krook_exponent = 0.0
comm_at_init = .false.
phi_bound = 0.0
phi_pow = 0.0
krook_efold = 3.0
use_dirichlet_bc = .false.
/
!--------------------------------
!------- Parallelisation ------
!--------------------------------
¶llelisation (renamed from &layout_knobs)
xyzs_layout = "yxzs"
vms_layout = "vms"
!! Double check
mat_gen = .false.
mat_read = .false.
lu_option = "default"
fields_kxkyz = .false.
/
&debug_flags
print_extra_info_to_terminal = .true.
debug_all = .false.
stella_debug = .false.
ffs_solve_debug = .false.
fields_all_debug = .false.
fields_debug = .false.
fields_fluxtube_debug = .false.
fields_electromagnetic_debug = .false.
fields_ffs_debug = .false.
implicit_solve_debug = .false.
parallel_streaming_debug = .false.
mirror_terms_debug = .false.
neoclassical_terms_debug = .false.
response_matrix_debug = .false.
time_advance_debug = .false.
extended_grid_debug = .false.
diagnostics_all_debug = .false.
diagnostics_parameters = .false.
diagnostics_fluxes_fluxtube_debug = .false.
diagnostics_omega_debug = .false.
diagnostics_debug = .false.
dist_fn_debug = .false.
gyro_averages_debug = .false.
fluxes_debug = .false.
geometry_debug = .false.
const_alpha_geo = .false.
/
I am happy with this, with a few minor exceptions. Here are my comments/questions:
-
How does include_electromagnetic work with the &electromagnetic namelist? E.g., if include_electromagnetic is false but include_apar is true (or vice versa), which one overrides the other?
-
I do not understand the comment about &adiabatic_electron_response ! Abort if kinetic electrons are included Will this still work in the same way as before, where adiabatic response is only used if no kinetic electron species is included? It will not actually abort if an adiabatic option is set, right?
-
What is the point of this namelist if the code automatically will choose the adiabatic option for you? &adiabatic_ion_response ! The <adiabatic_option> flag should not be set by the user ! Stella should be automatically choosing the option without the flux surface ! average for adiabatic ions and with the flux surface average for electrons /
-
How do you decide what has its own dedicated namelist? E.g., electromagnetic gets its own (despite really only having 2 unique flags/parameters), while flow shear gets lumped into 'extra_physics' despite having several associated flags/parameters. It might be worth including flow shear into the gyrokinetic_terms namelist and giving it a separate namelist for the associated flags/parameters; happy to chat this through with you, but basically the terms that are turned on/off are the perpendicular flow shear terms and the PVG term.
-
&euterpe_parameters is not related to neoclassical physics at all
-
I think the &sources namelist is related to radially global (not neoclassics)
I am happy with this, with a few minor exceptions. Here are my comments/questions:
- How does include_electromagnetic work with the &electromagnetic namelist? E.g., if include_electromagnetic is false but include_apar is true (or vice versa), which one overrides the other?
include_electromagnetic will override include_apar and include_bpar. If include_electromagnetic=.False. then the other two will be automatically set to .False. . If include_electromagnetic=.True. then the default is that the other two are true, but they are now able to be toggled individually to look at one or the other. We can put a message to screen about which physics is contained.
- I do not understand the comment about> &adiabatic_electron_response ! Abort if kinetic electrons are included Will this still work in the same way as before, where adiabatic response is only used if no kinetic electron species is included? It will not actually abort if an adiabatic option is set, right?
There was discussion about this. I think it would be nice to keep it as it is, where you can set the adiabatic response but if nspec=2 this is irrelevant. I think it might just be a case of adding a message to the screen about what is happening.
- What is the point of this namelist if the code automatically will choose the adiabatic option for you? &adiabatic_ion_response ! The <adiabatic_option> flag should not be set by the user ! Stella should be automatically choosing the option without the flux surface ! average for adiabatic ions and with the flux surface average for electrons
I believe this is added for our benefit only and will not be an available name list. This is because Cookie and others at CIEMAT are interested in multi-scale simulations in the future.
- How do you decide what has its own dedicated namelist? E.g., electromagnetic gets its own (despite really only having 2 unique flags/parameters), while flow shear gets lumped into 'extra_physics' despite having several associated flags/parameters. It might be worth including flow shear into the gyrokinetic_terms namelist and giving it a separate namelist for the associated flags/parameters; happy to chat this through with you, but basically the terms that are turned on/off are the perpendicular flow shear terms and the PVG term.
I think this is because Cookie and I are not 100% sure which flags correspond to what physics. We have currently put them under 'extra_physics' but we could make a separate name list for these, which would also help users in making it clear what these toggles do. We can chat and decide where to add them.
- &euterpe_parameters is not related to neoclassical physics at all
Will add.
- I think the &sources namelist is related to radially global (not neoclassics)
Will look into and change.
Potential feature:
- have maximal input file printed upon make with all default values in the code.
&extra_physics prp_shear_enabled = .false. hammett_flow_shear = .true. include_pressure_variation = .false. g_exb = 0.0 g_exbfac = 1.0 omprimfac = 1.0
changed to:
& flow_shear g_exb = 0.0 --- EXB shearing rate prp_shear_enabled = .false. --- Add shear in the EXB flow (ignore in example input with flow shear - should be true if g_exb and g_exbfac are non-zero) hammett_flow_shear = .true. --- Use wavenumber remapping algorithm for the EXB flow shear (ignore in example input with flow shear - should be the default algorithmic choice if flow shear is included) g_exbfac = 1.0 --- Set to 1.0 to include shear in perp. component of the toroidal flow (ignore in example input with flow shear - set default to 1.0) omprimfac = 1.0 --- set to 1.0 to include shear in the parallel component of the toroidal flow (ignore in example input with flow shear - set default to 1.0)
include_pressure_variation = .false. --- should be in radial variation.
To comment on the adiabatic species namelist, I think ideally one should be able to set:
adiabatic_electron_response = .true.
adiabatic_ion_reponse = .true.
In one of the namelists (obviously not both at the same time), where both are set to false by default. Georgia made a comment that it should be false if nspec=2, but be careful that Jose often runs with nspec=2 which are 2 impurities and then adding adiabatic electrons, so this assumption can not be made.
What I do not like about stella, is that it is not straightforward to add adiabatic ions or electrons. If we add adiabatic electrons, I think we should only set:
adiabatic_electron_response = .true.
tite = 1
nine = 1
And the adiabatic option in this case should be set to "field-line-average-term" by stella, it should not be set manually by the user. Similarly if we want to simulate turbulence on electron scales using adiabatic ions, we should be able to set:
adiabatic_ion_response = .true.
tite = 1
nine = 1
And the adiabatic option in this case should be set to "no-field-line-average-term".
Instead of the original suggestion for the namelist, perhaps we need a namelist:
&adiabatic_species
adiabatic_electron_response = .false.
adiabatic_ion_response = .false.
tite = 1
nine = 1
/
And if there are more input parameters which control the adiabatic species, they should go into this namelist. Right now I still find it very confusing on how to perform simulations with adiabatic ions.
In the current version of stella, to control the adiabatic species we require that none of the &species_parameters_X namelists have type = 'electron' and we need to set:
¶meters_physics
adiabatic_option = "field-line-average-term"
/
¶meters
nine = 1.0
tite = 1.0
/
Which is extremely confusing, also it is not clear that "nine" and "tite" are only used if we have adiabatic species. And to add adiabatic ions I think you need to use some hacks to go around this requirement,
These changes have been implemented in the upcoming stella release, I will close this issue.