Variable not found when using squashoutput

[georgeholt1]

I sometimes find that calling squashoutput.squashoutput fails when I try to call it on SD1D simulation output. This seems to only occur for simulations that ran in parallel. All simulations I've tested on exit gracefully with the creation of a BOUT.stop file.

The (truncated) error message is always:

    squashoutput(
  File "/marconi_work/FUA36_UKAEA_ML/gholt/miniconda3/envs/sd1d/lib/python3.9/site-packages/boutdata/squashoutput.py", line 207, in squashoutput
    var = outputs[varname]
  File "/marconi_work/FUA36_UKAEA_ML/gholt/miniconda3/envs/sd1d/lib/python3.9/site-packages/boutdata/data.py", line 1683, in __getitem__
    data = self._collect(name)
  File "/marconi_work/FUA36_UKAEA_ML/gholt/miniconda3/envs/sd1d/lib/python3.9/site-packages/boutdata/data.py", line 1484, in _collect
    return collect(
  File "/marconi_work/FUA36_UKAEA_ML/gholt/miniconda3/envs/sd1d/lib/python3.9/site-packages/boutdata/collect.py", line 237, in collect
    varname = findVar(varname, f.list())
  File "/marconi_work/FUA36_UKAEA_ML/gholt/miniconda3/envs/sd1d/lib/python3.9/site-packages/boutdata/collect.py", line 57, in findVar
    raise ValueError("Variable '" + varname + "' not found")
ValueError: Variable 'flux_ion' not found

Grepping for flux_ion in the dmp.nc files shows it is only present in one of them. However, this also seems to be the case for other similar simulations that don't throw the above error when running squashoutput on them.

My exact call to squashoutput is:

from boutdata.squashoutput import squashoutput
datadir = <path_to_directory>
squashoutput(datadir=datadir, tind=-1)

The BOUT.inp file for reproducibility is below. Changing neutral density from 0.0007 to 0.0001 results in a simulation that squashoutput works on.

#
#
# 

NOUT = 1600     # number of output time-steps
TIMESTEP = 5000.  # time between outputs

MZ = 1     # number of points in z direction (2^n + 1)
MXG = 0    # No guard cells needed in X

[mesh]

ny = 50    # Resolution along field-line

length = 25        # Length of the domain in meters
length_xpt = 12.5  # Length from midplane to X-point [m]
area_expansion = 1 # Expansion factor Area(target) / Area(midplane)

dy = length / ny   # Parallel grid spacing [m]

ypos = y * length / (2*pi) # Y position [m]

nx = 1
dx = 1
ixseps1 = -1   # Branch-cut indices, specifying that
ixseps2 = -1   # the grid is in the SOL

# The following make the field-aligned
# metric tensor an identity metric
Rxy = 1
Bpxy = 1
Btxy = 0
Bxy = 1
hthe = 1
sinty = 0

symmetricGlobalY = true

##################################################
# derivative methods

[ddy]

first = C2
second = C2
upwind = W3

[solver]

type = cvode

mxstep = 100000  # Maximum number of internal steps per output

atol = 1e-10
rtol = 1e-5

use_precon=true

[SD1D]

diagnose = true  # Output additional diagnostics

# Normalisation factors
Nnorm = 1e20  # Reference density [m^-3]
Tnorm = 100   # Reference temperature [eV]
Bnorm = 1.0   # Reference magnetic field [T]
AA = 2.0      # Ion atomic number

Eionize = 13.6     # Energy lost per ionisation [eV]
excitation = true  # Include electron-neutral excitation

volume_source = true   # Sources spread over a volume
density_upstream = 4e+19

density_controller_p = 1e-2 # Density controller 'p' parameter
density_controller_i = 1e-3 # Density controller 'i' parameter

# Model parameters
vwall = 0.0        # Velocity of neutrals at the wall, as fraction of Franck-Condon energy

frecycle = 0.95            # Recycling fraction
fredistribute = 0.0        # Fraction of recycled neutrals redistributed evenly along length
redist_weight = h(y - pi)  # Weighting for redistribution

gaspuff = 0        # NOTE: In normalised units 
dneut    = 10.0     # Scale neutral gas diffusion rate
nloss    = 1e3     # Neutral gas loss rate [1/s]
fimp = 0.05         # Impurity fraction
impurity_adas = true  # Use Atomicpp to get ADAS data
impurity_species = c  # Species. currently "c" (carbon) or "n" (Nitrogen)

sheath_gamma = 6 # Sheath heat transmission
neutral_gamma = 0.  # Neutral gas heat transmission
density_sheath = 1  # 0 = free, 1 = Neumann, 2 = constant nV
pressure_sheath = 1  # 0 = free, 1 = Neumann, 2 = constant (5/2)Pv + (1/2)nv^3

atomic = true      # Include atomic processes (CX/iz/rc)

# Set flux tube area as function of parallel grid index
# using normalised y coordinate from 0 to 2pi
area = 1 #+ (mesh:area_expansion - 1) * h(mesh:ypos - mesh:length_xpt)*(mesh:ypos - mesh:length_xpt)/(mesh:length - mesh:length_xpt)

hyper = -1 # Numerical diffusion parameter on all terms
ADpar = -1  # 4th-order numerical dissipation
viscos = -0.0001 # Parallel viscosity
ion_viscosity = false  # Braginskii parallel ion viscosity (ions and neutrals)

heat_conduction = true # Heat conduction

density_form = 4
momentum_form = 6
energy_form = 8

[All]
scale = 0.0

bndry_all = neumann_o2  # Default boundary condition
                        # Note: Sheath boundary applied in code

[Ne] # Electron density 
scale = 1

# Initial conditions
function = 0.1

flux = 4e23  # Particles per m^2 per second input
source = (flux/(mesh:length_xpt))*h(mesh:length_xpt - mesh:ypos)  # Particle input source
                           # as function of normalised y coordinate

[NVi]  # Parallel ion momentum
scale = 1
vtarg = 0.3
function = vtarg * Ne:scale * Ne:function * y / (2*pi)  # Linear from 0 to 0.03 in y
bndry_target = dirichlet_o2

[P]    # Plasma pressure P = 2 * Ne * T
scale = 1
function=0.1   # Initial constant pressure

powerflux = 2e7  # Input power flux in W/m^2

source = (powerflux*2/3 / (mesh:length_xpt))*h(mesh:length_xpt - mesh:ypos)  # Input power as function of y

# 1e7 W/m^2 / (L/2) with  L = 100 m , factor of 2 because power only from y = 0 to y=pi
# * 2/3 to get from energy to Pe

[Nn]
# Neutral density
scale = 1
function = 0.0007

[NVn]
evolve = true # Evolve neutral momentum?

[Pn]
evolve = true # Evolve neutral pressure? Otherwise Tn = Te model

Tstart = 3.5 # Starting temperature in eV

scale = 1.0
function = Nn:scale * Nn:function * Pn:Tstart / SD1D:Tnorm

[johnomotani]

It looks like flux_ion is a scalar, but is only written on the process that has the 'upper-y' boundary.

[BTW I have no idea how changing the neutral density could affect this error!]

The assumption in squashoutput is that scalars are things that are the same on all processes, and would be written on all processes, so scalars will be read from (I think) BOUT.dmp.0.nc. This usually makes sense - e.g. a 'flux to the target' in a 3d tokamak simulation would be needed at up to 4 separate targets so there would be no way to collect those into a single scalar variable.

I suspect a solution requires a change to some code, either boutdata.squashoutput or SD1D (or both):

• the 'best' solution is probably to

  1. make ion_flux a FieldPerp object in SD1D
  2. call SAVE_REPEAT(ion_flux); on all processes, not just the last one
  3. on the process where ion_flux is actually needed, call ion_flux.set_index() to give it a sensible y-index (i.e. a local index in range for that process, I think either a grid cell or boundary cell, but not a guard cell).

  That should let squashoutput() (via BoutOutputs) handle ion_flux correctly, although the FieldPerp handling is a tricky special case and not as well tested as scalar/Field2D/Field3D.

• a workaround (and possibly a useful feature anyway) would be to give squash_output() an argument like drop_variables, which would be a list of names that it should just ignore if they are found. I'll make a PR for this - it should be simple...

[georgeholt1]

Thanks for the thorough explanation and quickfix!

I tried to run squashoutput on the simulation output that it previously worked for and this time it threw the above error. It could be intermittent, or I could be making a mistake somewhere.