Creating your first POPCON

Welcome to cfspopcon!

This notebook will work you through how to set up a POPCON run and plot the results.

To start with, execute the cell below to import some additional libraries that we’ll frequently use, as well as the cfspopcon library itself. If you can execute this cell without error, that’s a great start. If not, make sure that the cfspopcon library is installed correctly.

[1]:

%reload_ext autoreload
%autoreload 2

import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path

import cfspopcon
from cfspopcon.unit_handling import ureg

# Change to the top-level directory. Required to find radas_dir in its default location.
%cd {Path(cfspopcon.__file__).parents[1]}
# As a sanity check, print out the current working directory
print(f"Running in {Path('').absolute()}")

/Users/tbody/Projects/cfspopcon
Running in /Users/tbody/Projects/cfspopcon

The simplest way to configure a POPCON analysis is via a configuration file. For this, we use YAML files (which map neatly to python dictionaries).

In the next cell, we read in the SPARC_PRD/input.yaml file, applying a few conversions on the data. This includes

  1. The algorithms entry is converted into a cfspopcon.CompositeAlgorithm which we’ll talk about later. This basically gives the list of operations that we want to perform on the input data.

  2. The points entry is stored in a separate dictionary. This gives a set of key-value pairs of ‘optimal’ points (for instance, giving the point with the maximum fusion power gain).

  3. The grids entry is converted into an xr.DataArray storing a np.linspace or np.logspace of values which we scan over. We usually scan over average_electron_density and average_electron_temp, but there’s nothing preventing you from scanning over other numerical input variables or having more than 2 dimensions which you scan over (n.b. this can get expensive!).

  4. Each input variable is checked to see if its name matches one of the enumerators in cfspopcon.named_options. These are used to store switch values, such as cfspopcon.named_options.ReactionType.DT which indicates that we’re interested in the DT fusion reaction.

  5. Each input variable is converted into its default units. Default units are retrieved via the cfspopcon.unit_handling.default_unit function. This will set, for instance, the average_electron_temp values to have units of keV.

[2]:

input_parameters, algorithm, points, plots  = cfspopcon.read_case("example_cases/SPARC_PRD")

We first want to make sure that our analysis won’t crash due to missing input variables. For this, we use the validate_inputs method of the cfspopcon.CompositeAlgorithm object, which makes sure that all required input variables are defined by the input.yaml file. It will also tell you if the algorithms you’re requesting are out of order.

[3]:

algorithm.validate_inputs(input_parameters);

Since we’re looking at the example, it’s unsurprisingly well behaved. Let’s intentionally drop an input variable to show that it works.

[4]:

incomplete_input_parameters = input_parameters.copy()
del incomplete_input_parameters["major_radius"]

algorithm.validate_inputs(incomplete_input_parameters, raise_error_on_missing_inputs=False);

/var/folders/x2/fhfghwm566d83ws71kgzj86c0000gp/T/ipykernel_90477/170706882.py:4: UserWarning: Missing input parameters [major_radius].
  algorithm.validate_inputs(incomplete_input_parameters, raise_error_on_missing_inputs=False);

Next, we pack all of the (complete) input parameters into an xarray.Dataset.

If you’re not familiar with xarray, a Dataset functions like a dictionary which stores labelled xr.DataArrays, which in turn function like np.arrays but with additional annotation to describe the coordinates and units of the array. We’ll use a lot of xarray features, so it’s worth reading the docs, including their useful How do I … section.

Our starting dataset isn’t super interesting: it’s exactly what we defined in the input.yaml file. In the next cell, we construct a dataset from out input parameters, and then print into the notebook a representation of the notebook (this is the somewhat odd-looking second line).

[5]:

dataset = xr.Dataset(input_parameters)

dataset

[5]:

<xarray.Dataset> Size: 2kB
Dimensions:                               (dim_species: 3,
                                           dim_average_electron_density: 40,
                                           dim_average_electron_temp: 30)
Coordinates:
  * dim_species                           (dim_species) object 24B AtomicSpec...
  * dim_average_electron_density          (dim_average_electron_density) float64 320B ...
  * dim_average_electron_temp             (dim_average_electron_temp) float64 240B ...
Data variables: (12/48)
    major_radius                          float64 8B [m] 1.85
    magnetic_field_on_axis                float64 8B [T] 12.2
    inverse_aspect_ratio                  float64 8B [] 0.3081
    areal_elongation                      float64 8B [] 1.75
    elongation_ratio_sep_to_areal         float64 8B [] 1.125
    elongation_ratio_areal_to_psi95       float64 8B [] 1.025
    ...                                    ...
    total_flux_available_from_CS          int64 8B [Wb] 35
    ejima_coefficient                     float64 8B [] 0.6
    safety_factor_on_axis                 int64 8B [] 1
    vertical_magnetic_field_equation      <U4 16B 'Barr'
    average_electron_density              (dim_average_electron_density) float64 320B [1e19 m^-3] ...
    average_electron_temp                 (dim_average_electron_temp) float64 240B [keV] ...

Let’s run a POPCON!

We’ll use the update_dataset method of our algorithm and pass in our dataset of input parameters.

The CompositeAlgorithm calls the update_dataset for each of its Algorithms. These calls check that all of the required inputs are defined in the dataset, and then write the results back into the dataset.

[6]:

dataset = algorithm.update_dataset(dataset)

dataset

[6]:

<xarray.Dataset> Size: 3MB
Dimensions:                                (dim_average_electron_temp: 30,
                                            dim_average_electron_density: 40,
                                            dim_species: 5, dim_rho: 50)
Coordinates:
  * dim_average_electron_temp              (dim_average_electron_temp) float64 240B ...
  * dim_average_electron_density           (dim_average_electron_density) float64 320B ...
  * dim_species                            (dim_species) object 40B AtomicSpe...
Dimensions without coordinates: dim_rho
Data variables: (12/148)
    current_relaxation_time                (dim_average_electron_temp, dim_average_electron_density) float64 10kB [s] ...
    peak_pressure                          (dim_average_electron_temp, dim_average_electron_density) float64 10kB [Pa] ...
    fusion_triple_product                  (dim_average_electron_density, dim_average_electron_temp) float64 10kB [1e10 m^-3·keV·s] ...
    rho_star                               (dim_average_electron_temp) float64 240B [] ...
    nu_star                                (dim_average_electron_density, dim_average_electron_temp) float64 10kB [] ...
    core_radiated_power_fraction           (dim_average_electron_density, dim_average_electron_temp) float64 10kB [] ...
    ...                                     ...
    total_flux_available_from_CS           int64 8B [Wb] 35
    ejima_coefficient                      float64 8B [] 0.6
    safety_factor_on_axis                  int64 8B [] 1
    vertical_magnetic_field_equation       <U4 16B 'Barr'
    average_electron_density               (dim_average_electron_density) float64 320B [1e19 m^-3] ...
    average_electron_temp                  (dim_average_electron_temp) float64 240B [keV] ...

Plotting and interrogating the results

That’s it: the only thing left is to look at our results. We can use the built-in plotting functionality, which is configured by dictionaries which we read in from YAML files using read_plot_style.

Nothing fancy is happening here: it’s just reading in a dictionary from a file and not doing any conversions. You can modify the YAML file, or directly make a dictionary yourself. The structure of the dictionary is

  • type: what sort of plot do you want to make. Currently we only support popcon, but room for growth.

  • figsize: the size of the figure in inches, defining both the aspect ratio and text size of the resulting plot. If you want larger overall labels, try reducing figsize

  • show_dpi: the pixels-per-inch of the resulting figure. Increase this if your figure is blurry. This will also make the figure larger in the Jupyter notebook.

  • save_as: a name for the figure when saving it as a .png

  • coords or new_coords: what to use as the x-axis and y-axis. coords must be already in the coordinates of the dataset, while new_coords can be output variables (see the example below)

  • fill: a block for a variable to be plotted as a colormesh. Inside this, we define a section where which we use to build a mask (discussed below)

  • points: a block for points to be highlighted on the plot. These must correspond to points defined in the input.yaml file

  • contour: a block for variables which we should plot contours for. We pick a single color per contour and then label the contour lines themselves at specific values, which lets us show multiple variables at once

[7]:

plot_style = cfspopcon.read_plot_style("example_cases/SPARC_PRD/plot_popcon.yaml")

plot_style

[7]:

{'type': 'popcon',
 'figsize': [8, 6],
 'show_dpi': 150,
 'coords': {'x': {'dimension': 'average_electron_temp',
   'label': '$<T_e>$ [$keV$]',
   'units': 'keV'},
  'y': {'dimension': 'average_electron_density',
   'label': '$<n_e>$ [$10^{20} m^{-3}$]',
   'units': 'n20'}},
 'fill': {'variable': 'Q',
  'where': {'Q': {'min': 1.0},
   'P_auxiliary_launched': {'min': 0.0, 'max': 25.0, 'units': 'MW'},
   'greenwald_fraction': {'max': 0.9},
   'ratio_of_P_SOL_to_P_LH': {'min': 1.0},
   'max_flattop_duration': {'min': 0.0, 'units': 'seconds'}}},
 'points': {'PRD': {'label': 'PRD',
   'marker': 'x',
   'color': 'red',
   'size': 50.0}},
 'contour': {'Q': {'label': '$Q$',
   'levels': [0.1, 1.0, 2.0, 5.0, 10.0, 50.0],
   'color': 'tab:red',
   'format': '1.2g'},
  'ratio_of_P_SOL_to_P_LH': {'label': '$P_{SOL}/P_{LH}$',
   'color': 'tab:blue',
   'levels': [1.0],
   'format': '1.2g'},
  'P_auxiliary_launched': {'label': '$P_{aux,launched}$',
   'levels': [1.0, 5.0, 10.0, 25.0, 50.0],
   'color': 'tab:gray',
   'format': '1.2g'},
  'P_fusion': {'label': '$P_{fusion}$',
   'color': 'tab:purple',
   'levels': [50.0, 100.0, 150.0, 200.0],
   'format': '1.2g'},
  'max_flattop_duration': {'label': '$t_{flattop}$',
   'color': 'tab:orange',
   'levels': [0.0, 1.0, 2.0, 4.0, 6.0, 8.0, 10.0],
   'format': '1.2g'}}}

Once we have our plot_style, we call cfspopcon.plotting.make_plot to make a figure.

N.b. we’ve set output_dir=None to stop make_plot from saving the figure to a .png file. If you’re using the command-line tool, you can delete the save_as key in plot_style.

It’s pretty easy to make this plot, but not entirely straightforward to interpret it.

Firstly, the axes are for the average density and for the average temperature. These are more often the outputs of a calculation, rather than the independent variables. This is a consequence of the particular ‘back to front’ algorithm that we’re using, which uses the average temperature and density to define the plasma stored energy and the resulting confinement time and input power required to stay at that point.

Next, it’s helpful to remember that all of the points shown are completely independent. Each represents a time-independent solution of the given algorithm.

Although we get a solution at all \((\bar n_e, \bar T_e)\) points, they don’t all make sense as operational points. There are a set of limits that we impose via the where block inside the fill block. For this example, these are

  1. We’re looking for high fusion gains, with at least \(Q>1\),

  2. Since we’re using the ITER98y2 H-mode energy confinement scaling, we require that \(P_{SOL} > P_{LH}\),

  3. We can’t switch off the Ohmic power, so we neglect points with \(P_{RF} <= 0\),

  4. To avoid density-limit disruptions, we require that the volume-averaged density is no more than 90% of the Greenwald density limit,

  5. To avoid damaging the magnets with neutrons, we impose a limit of \(P_{fusion} < 140MW\).

We mask points which don’t meet these requirements, and then plot a variable over the unmasked region (in this case, we’re plotting \(Q\)).

On top of this, we plot contours of fields of interest such as the divertor heat-flux metric \(P_{SOL} B_{pol} / R n^2\) so we can see the trends of these field as we go across the accessible parameter space.

Finally, within the accessible parameter, we select points which minimize or maximize the values of fields, such as the point giving the maximum value of \(Q\).

[8]:

cfspopcon.plotting.make_plot(
    dataset,
    plot_style,
    points=points,
    title="POPCON example",
    save_name=None
)
../_images/doc_sources_getting_started_15_0.png

There’s lots of functionality which you can explore (and extend!). For example, if you’d prefer to look at the results in terms of the auxiliary heating power instead of the average electron temperature, you can use cfspopcon.transform to map the results onto new axes.

Some of these features are a bit experimental. For instance, you can see here that by remapping we’ve shrunk the range of \(Q\) in the accessible parameter space, since we apply a transformation on the masked array. Increasing the resolution of the analysis grid should reduce the difference when remapping.

[9]:

cfspopcon.plotting.make_plot(dataset,
                             cfspopcon.read_plot_style("example_cases/SPARC_PRD/plot_remapped.yaml"),
                             points=points,
                             title="POPCON example",
                             save_name=None,
                            )
../_images/doc_sources_getting_started_17_0.png

We often want to study individual points in detail.

The simplest way to do this is to use file_io.write_point_to_file which writes the values at the points specified in the points dictionary to a JSON file.

[10]:

for point, point_params in points.items():
    cfspopcon.file_io.write_point_to_file(dataset, point, point_params, output_dir=Path("example_cases/SPARC_PRD/output"))

You can dig into the output file SPARC_PRD/outputs/max_Q.json, or alternatively we can peek under the hood to understand what is going on. Let’s pick the max_Q point which we defined in input.yaml.

[11]:

point_params = points["PRD"]

point_params

[11]:

{'maximize': 'Q',
 'where': {'P_auxiliary_launched': {'min': 0.0, 'max': 25.0, 'units': 'MW'},
  'greenwald_fraction': {'max': 0.9},
  'ratio_of_P_SOL_to_P_LH': {'min': 1.0},
  'P_fusion': {'max': 140.0, 'units': 'MW'}}}

We first build a mask from the where key. This lets us hide parts of operational space which aren’t physical (for example, if they have \(P_{SOL}\) below the LH transition).

This can be used with the where method of an xarray.DataArray.

[12]:

mask = cfspopcon.shaping_and_selection.build_mask_from_dict(dataset, point_params)

dataset.Q.where(mask).plot()

[12]:

<matplotlib.collections.QuadMesh at 0x32fde29f0>
../_images/doc_sources_getting_started_23_1.png

Next, let’s find the point here with the highest value of Q. This gives us the indices of the corresponding point.

[13]:

point_indices = cfspopcon.shaping_and_selection.find_coords_of_maximum(dataset.Q, mask=mask)

point_indices

[13]:

{'dim_average_electron_density': <xarray.DataArray 'Q' ()> Size: 8B
 array(24),
 'dim_average_electron_temp': <xarray.DataArray 'Q' ()> Size: 8B
 array(8)}

We can use those indices to select a point of the array, see the coordinates and value at that point.

[14]:

dataset.Q.isel(point_indices)

[14]:

<xarray.DataArray 'Q' ()> Size: 8B
<Quantity(12.720160972996633, 'dimensionless')>
Coordinates:
    dim_average_electron_temp     float64 8B 9.138
    dim_average_electron_density  float64 8B 25.0

We can also select this point for all of the variables defined in our dataset, which lets us look at the values of all of the different fields at that point (you might need to click the arrow next to “data variables” to see these).

[15]:

point = dataset.isel(point_indices)

point

[15]:

<xarray.Dataset> Size: 3kB
Dimensions:                                (dim_species: 5, dim_rho: 50)
Coordinates:
    dim_average_electron_temp              float64 8B 9.138
    dim_average_electron_density           float64 8B 25.0
  * dim_species                            (dim_species) object 40B AtomicSpe...
Dimensions without coordinates: dim_rho
Data variables: (12/148)
    current_relaxation_time                float64 8B [s] 16.36
    peak_pressure                          float64 8B [Pa] 2.469e+06
    fusion_triple_product                  float64 8B [1e10 m^-3·keV·s] 45.28
    rho_star                               float64 8B [] 0.002213
    nu_star                                float64 8B [] 0.02375
    core_radiated_power_fraction           float64 8B [] 0.1872
    ...                                     ...
    total_flux_available_from_CS           int64 8B [Wb] 35
    ejima_coefficient                      float64 8B [] 0.6
    safety_factor_on_axis                  int64 8B [] 1
    vertical_magnetic_field_equation       <U4 16B 'Barr'
    average_electron_density               float64 8B [1e19 m^-3] 25.0
    average_electron_temp                  float64 8B [keV] 9.138

If we want to look at a particular value, we can pull the single-element DataArray out of the Dataset.

For instance, we can return how much power is crossing the separatrix. By default, this is stored in megawatts, but we can easily convert it to other compatible units such as watts.

[16]:

point["power_crossing_separatrix"].pint.to(ureg.watt)

[16]:

<xarray.DataArray 'power_crossing_separatrix' ()> Size: 8B
<Quantity(25574170.517921906, 'watt')>
Coordinates:
    dim_average_electron_temp     float64 8B 9.138
    dim_average_electron_density  float64 8B 25.0

When we selected the point, we selected a single average electron density and temperature, but left the rest of the structure of the data unchanged. If we’d defined more than 2 variables in our grid section of our input.yaml, we could look at how the point (in this case maximum Q) varies as a function of the third parameter (which could be something like the confinement time scalar \(H\)). Even for our simple 2D POPCON, there’s still some 1D data at our point: namely, the assumed temperature and density profiles.

[17]:

point.electron_temp_profile.pint.to(ureg.keV).plot()

[17]:

[<matplotlib.lines.Line2D at 0x32febccb0>]
../_images/doc_sources_getting_started_33_1.png

That’s the basics for using the POPCON library.

You can perform the analysis in this notebook using the command line tool. If you wanted to run this exact case from a terminal at the top-level of this repository, the command would be:

poetry run popcon example_cases/SPARC_PRD --show

where example_cases/SPARC_PRD is the path to the input.yaml file.

Feel free to try things out and change the input.yaml file. At some point, you’ll probably want to dig into the formulas and algorithms and start implementing your own. For this, see the other notebooks in this folder.