Evaluating a single PROCESS model¶
When understanding or investigating an individual model within Process, it can be useful to run the model in isolation and plot some responses. This is done here to investigate the effect of tungsten impurity concentration on radiated power and power incident on the divertor.
In [1]:
Copied!
import matplotlib.pyplot as plt
import numpy as np
import process
from process.main import SingleRun, setup_loggers
# setup the loggers so that the output is not spammed with model errors/warnings
setup_loggers()
import matplotlib.pyplot as plt
import numpy as np
import process
from process.main import SingleRun, setup_loggers
# setup the loggers so that the output is not spammed with model errors/warnings
setup_loggers()
Set up¶
First, inspect a variable to check its uninitialised value:
In [2]:
Copied!
print(
f"p_plasma_separatrix_mw = {process.data_structure.physics_variables.p_plasma_separatrix_mw}"
)
print(
f"p_plasma_separatrix_mw = {process.data_structure.physics_variables.p_plasma_separatrix_mw}"
)
p_plasma_separatrix_mw = None
In order to initialise all variables in Process with their values at a given point (design parameter vector), run an evaluation input file (one with no optimisation) to initialise values in all models. The "large tokamak" regression test solution is used here.
In [3]:
Copied!
single_run = SingleRun("data/large_tokamak_eval_IN.DAT")
single_run.run()
single_run = SingleRun("data/large_tokamak_eval_IN.DAT")
single_run.run()
The IN.DAT file does not contain any obsolete variables. ************************************************************************************************************** ************************************************** PROCESS *************************************************** ************************************** Power Reactor Optimisation Code *************************************** ************************************************************************************************************** Version : 3.2.1 Git Tag : v3.2.2-14-g1448faf3 Git Branch : main Date : 05/12/2025 UTC Time : 18:10 User : runner Computer : runnervmoqczp Directory : /home/runner/work/PROCESS/PROCESS/examples Input : /home/runner/work/PROCESS/PROCESS/examples/data/large_tokamak_eval_IN.DAT Run title : generic large tokamak Run type : Reactor concept design: Pulsed tokamak model, (c) UK Atomic Energy Authority ************************************************************************************************************** Equality constraints : 2 Inequality constraints : 23 Total constraints : 25 Iteration variables : 2 ************************************************************************************************************** Solving equality constraints using fsolve
/home/runner/work/PROCESS/PROCESS/process/init.py:92: UserWarning: temp_cs_superconductor_margin_min and tmargmin should not both be specified in IN.DAT temp_cs_superconductor_margin_min has been ignored check_process(inputs)
/home/runner/work/PROCESS/PROCESS/process/physics.py:7764: RuntimeWarning: divide by zero encountered in scalar divide * (nd_plasma_pedestal_electron / n_greenwald) ** -0.174
/home/runner/work/PROCESS/PROCESS/process/costs.py:3033: RuntimeWarning: invalid value encountered in sqrt annoam = cost_variables.ucoam[cost_variables.lsa - 1] * np.sqrt( /home/runner/work/PROCESS/PROCESS/process/costs.py:3104: RuntimeWarning: invalid value encountered in sqrt annwst = cost_variables.ucwst[cost_variables.lsa - 1] * np.sqrt(
************************************ PROCESS found a consistent solution *************************************
******************************************** Errors and Warnings ********************************************* (/home/runner/work/PROCESS/PROCESS/process/physics.py:6725) Diamagnetic fraction is more than 1%, but not calculated. Consider using i_diamagnetic_current=2 and i_pfirsch_schluter_current=1 (/home/runner/work/PROCESS/PROCESS/process/tf_coil.py:246) dr_tf_plasma_case to small to accommodate the WP, forced to minimum value (/home/runner/work/PROCESS/PROCESS/process/pfcoil.py:2385) CS not using max current density: further optimisation may be possible ******************************************* End of PROCESS Output ********************************************
In [4]:
Copied!
# Kernel can crash when running physics without correctly initialised values
# Doesn't crash after running a once-through
# Print initial values of interest
def print_values():
print(
f"W frac = {process.data_structure.impurity_radiation_module.f_nd_impurity_electron_array[13]:.3e}"
)
print(
f"p_plasma_rad_mw = {process.data_structure.physics_variables.p_plasma_rad_mw:.3e}"
)
print(
f"p_plasma_separatrix_mw = {process.data_structure.physics_variables.p_plasma_separatrix_mw:.3e}"
)
print_values()
# Kernel can crash when running physics without correctly initialised values
# Doesn't crash after running a once-through
# Print initial values of interest
def print_values():
print(
f"W frac = {process.data_structure.impurity_radiation_module.f_nd_impurity_electron_array[13]:.3e}"
)
print(
f"p_plasma_rad_mw = {process.data_structure.physics_variables.p_plasma_rad_mw:.3e}"
)
print(
f"p_plasma_separatrix_mw = {process.data_structure.physics_variables.p_plasma_separatrix_mw:.3e}"
)
print_values()
W frac = 5.000e-06 p_plasma_rad_mw = 2.150e+02 p_plasma_separatrix_mw = 1.734e+02
Now try increasing the tungsten impurity fraction to see if there's a change in the divertor power.
In [5]:
Copied!
process.data_structure.impurity_radiation_module.f_nd_impurity_electron_array[13] = (
5.0e-5
)
single_run.models.physics.physics()
print_values()
process.data_structure.impurity_radiation_module.f_nd_impurity_electron_array[13] = (
5.0e-5
)
single_run.models.physics.physics()
print_values()
W frac = 5.000e-05 p_plasma_rad_mw = 2.729e+02 p_plasma_separatrix_mw = 1.144e+02
With a higher W impurity fraction, the radiated power has increased and the power incident on the divertor has decreased.
Parameter study of W impurity¶
Now investigate effect of varying W impurity on impurity radiation power, divertor power and constraint 15 (L-H threshold constraint).
In [6]:
Copied!
from process.constraints import ConstraintManager
def run_impurities(w_imp_fracs):
"""Calculate responses to W impurities."""
n = w_imp_fracs.shape[0]
p_plasma_rad_mw = np.empty(n)
p_plasma_separatrix_mw = np.empty(n)
p_l_h_threshold_mw = np.empty(n)
con15 = np.empty(n)
# Loop over W impurity values, evaluate model and store responses at each point
for i, imp_frac in enumerate(w_imp_fracs):
# Set W impurity fraction, then run physics model
process.data_structure.impurity_radiation_module.f_nd_impurity_electron_array[
13
] = imp_frac
single_run.models.physics.physics()
# Evaluate constraint equation 15 (L-H threshold constraint)
con15_value = ConstraintManager.evaluate_constraint(15).normalised_residual
# Need to copy values
p_plasma_rad_mw[i] = (
process.data_structure.physics_variables.p_plasma_rad_mw.item()
)
p_plasma_separatrix_mw[i] = (
process.data_structure.physics_variables.p_plasma_separatrix_mw.item()
)
p_l_h_threshold_mw[i] = (
process.data_structure.physics_variables.p_l_h_threshold_mw.item()
)
# Need to flip sign of constraint so negative means violated
con15[i] = -con15_value
return p_plasma_rad_mw, p_plasma_separatrix_mw, p_l_h_threshold_mw, con15
from process.constraints import ConstraintManager
def run_impurities(w_imp_fracs):
"""Calculate responses to W impurities."""
n = w_imp_fracs.shape[0]
p_plasma_rad_mw = np.empty(n)
p_plasma_separatrix_mw = np.empty(n)
p_l_h_threshold_mw = np.empty(n)
con15 = np.empty(n)
# Loop over W impurity values, evaluate model and store responses at each point
for i, imp_frac in enumerate(w_imp_fracs):
# Set W impurity fraction, then run physics model
process.data_structure.impurity_radiation_module.f_nd_impurity_electron_array[
13
] = imp_frac
single_run.models.physics.physics()
# Evaluate constraint equation 15 (L-H threshold constraint)
con15_value = ConstraintManager.evaluate_constraint(15).normalised_residual
# Need to copy values
p_plasma_rad_mw[i] = (
process.data_structure.physics_variables.p_plasma_rad_mw.item()
)
p_plasma_separatrix_mw[i] = (
process.data_structure.physics_variables.p_plasma_separatrix_mw.item()
)
p_l_h_threshold_mw[i] = (
process.data_structure.physics_variables.p_l_h_threshold_mw.item()
)
# Need to flip sign of constraint so negative means violated
con15[i] = -con15_value
return p_plasma_rad_mw, p_plasma_separatrix_mw, p_l_h_threshold_mw, con15
In [7]:
Copied!
%matplotlib inline
# Run W impurity parameter study
w_imp_fracs = np.linspace(1.0e-6, 1.0e-4, 50)
p_plasma_rad_mw, p_plasma_separatrix_mw, p_l_h_threshold_mw, con15 = run_impurities(
w_imp_fracs
)
fig, ax = plt.subplots()
ax.scatter(w_imp_fracs, p_plasma_rad_mw, label="p_plasma_rad_mw")
ax.scatter(w_imp_fracs, p_plasma_separatrix_mw, label="p_plasma_separatrix_mw")
ax.plot(w_imp_fracs, p_l_h_threshold_mw, "r", label="p_l_h_threshold_mw")
ax.ticklabel_format(style="sci", axis="x", scilimits=(0, 0))
ax.set_title("W impurity fraction against radiated and divertor power")
ax.set_xlabel("W impurity fraction")
ax.set_ylabel("Power (MW)")
ax.legend()
%matplotlib inline
# Run W impurity parameter study
w_imp_fracs = np.linspace(1.0e-6, 1.0e-4, 50)
p_plasma_rad_mw, p_plasma_separatrix_mw, p_l_h_threshold_mw, con15 = run_impurities(
w_imp_fracs
)
fig, ax = plt.subplots()
ax.scatter(w_imp_fracs, p_plasma_rad_mw, label="p_plasma_rad_mw")
ax.scatter(w_imp_fracs, p_plasma_separatrix_mw, label="p_plasma_separatrix_mw")
ax.plot(w_imp_fracs, p_l_h_threshold_mw, "r", label="p_l_h_threshold_mw")
ax.ticklabel_format(style="sci", axis="x", scilimits=(0, 0))
ax.set_title("W impurity fraction against radiated and divertor power")
ax.set_xlabel("W impurity fraction")
ax.set_ylabel("Power (MW)")
ax.legend()
Out[7]:
<matplotlib.legend.Legend at 0x7f2f786c0790>
How does the L-H threshold constraint vary?
In [8]:
Copied!
%matplotlib inline
fig, ax = plt.subplots()
ax.scatter(w_imp_fracs, con15, label="con15")
ax.ticklabel_format(style="sci", axis="x", scilimits=(0, 0))
ax.set_title("W impurity fraction against L-H threshold constraint (15)")
ax.set_xlabel("W impurity fraction")
ax.set_ylabel("L-H threshold constraint value")
ax.hlines(0.0, xmin=ax.get_xlim()[0], xmax=ax.get_xlim()[1], colors="r")
ax.annotate("Satisfied", (0.0, 0.1))
ax.annotate("Violated", (0.0, -0.15))
%matplotlib inline
fig, ax = plt.subplots()
ax.scatter(w_imp_fracs, con15, label="con15")
ax.ticklabel_format(style="sci", axis="x", scilimits=(0, 0))
ax.set_title("W impurity fraction against L-H threshold constraint (15)")
ax.set_xlabel("W impurity fraction")
ax.set_ylabel("L-H threshold constraint value")
ax.hlines(0.0, xmin=ax.get_xlim()[0], xmax=ax.get_xlim()[1], colors="r")
ax.annotate("Satisfied", (0.0, 0.1))
ax.annotate("Violated", (0.0, -0.15))
Out[8]:
Text(0.0, -0.15, 'Violated')
The constraint becomes violated for W fraction values $> 6\times10^{-5}$.
This can easily be modified to investigate behaviour of any model in Process in isolation, without running other models or optimising.