Merge branch 'main' into 3388-update-the-beta-limit-section-of-the-do…

…cs-to-be-up-to-date-and-include-all-models
ukaea · Jan 6, 2025 · 1023e2d · 1023e2d
2 parents 906dddb + bc0dc68
commit 1023e2d
Show file tree

Hide file tree

Showing 11 changed files with 109 additions and 792 deletions.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -90,7 +90,6 @@ LIST(APPEND PROCESS_SRCS
     buildings_variables.f90
     maths_library.f90
     iteration_variables.f90
-    evaluators.f90
     water_usage_variables.f90
     constraint_equations.f90
     constants.f90

diff --git a/README.md b/README.md
@@ -6,7 +6,9 @@
 
 PROCESS is the reactor systems code at the [UK Atomic Energy Authority](https://ccfe.ukaea.uk/). More information on PROCESS can be found on the PROCESS [webpage](https://ccfe.ukaea.uk/resources/process/).
 
-PROCESS was originally a Fortran code, but is currently a mixture of Python and Python-wrapped Fortran; the eventual aim is to have an entirely Python code base. In order to use PROCESS, the Fortran must be compiled and a Python-Fortran interface generated for the Python to import. Once built, it can be installed and run as a Python package.
+PROCESS was originally a Fortran code, but is currently a mixture of Python and Python-wrapped Fortran; the eventual aim is to have an entirely Python code base. In order to use PROCESS, the Fortran must be compiled and a Python-Fortran interface generated for the Python to import. Once built, it can be installed and run as a Python package. 
+
+Due to this ongoing conversion work, **PROCESS version 3 is unstable and does not guarantee backward compatibility**. PROCESS version 4 will be the first major version to enforce backward-compatible API changes and will be released following confirmation of variable names and the Python data structure. 
 
 
 

diff --git a/process/caller.py b/process/caller.py
@@ -3,6 +3,7 @@
 import numpy as np
 import logging
 from process.final import finalise
+from process.objectives import objective_function
 from process.io.mfile import MFile
 from process.utilities.f2py_string_patch import f2py_compatible_to_string
 from typing import Union, Tuple, TYPE_CHECKING
@@ -70,7 +71,7 @@ def call_models(self, xc: np.ndarray, m: int) -> Tuple[float, np.ndarray]:
         for _ in range(10):
             self._call_models_once(xc)
             # Evaluate objective function and constraints
-            objf = ft.function_evaluator.funfom()
+            objf = objective_function(ft.numerics.minmax)
             conf, _, _, _, _ = ft.constraints.constraint_eqns(m, -1)
 
             if objf_prev is None and conf_prev is None:

diff --git a/process/final.py b/process/final.py
@@ -5,11 +5,11 @@
 from process.fortran import (
     process_output as po,
     constants,
-    function_evaluator,
     numerics,
     constraints,
 )
 from process import output as op
+from process.objectives import objective_function
 from process.utilities.f2py_string_patch import f2py_compatible_to_string
 
 
@@ -50,7 +50,7 @@ def output_once_through():
     po.oblnkl(constants.nout)
 
     # Evaluate objective function
-    norm_objf = function_evaluator.funfom()
+    norm_objf = objective_function(numerics.minmax)
     po.ovarre(
         constants.mfile, "Normalised objective function", "(norm_objf)", norm_objf
     )

diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py
@@ -462,13 +462,13 @@ def color_key(axis, mfile_data, scan, colour_scheme):
                 [0.7, -0.3], 1, 0.4, lw=0, facecolor=CRYOSTAT_COLOUR[colour_scheme - 1]
             )
         )
-
-    axis.text(-5, 1, "Cryostat", ha="left", va="top", size="medium")
-    axis.add_patch(
-        patches.Rectangle(
-            [0.7, 0.7], 1, 0.1, lw=0, facecolor=CRYOSTAT_COLOUR[colour_scheme - 1]
+    else:
+        axis.text(-5, 1, "Cryostat", ha="left", va="top", size="medium")
+        axis.add_patch(
+            patches.Rectangle(
+                [0.7, 0.7], 1, 0.1, lw=0, facecolor=CRYOSTAT_COLOUR[colour_scheme - 1]
+            )
         )
-    )
 
 
 def toroidal_cross_section(axis, mfile_data, scan, demo_ranges, colour_scheme):
@@ -883,44 +883,6 @@ def read_imprad_data(skiprows, data_path):
     return impdata
 
 
-def synchrotron_rad():
-    """Function for Synchrotron radiation power calculation from Albajar, Nuclear Fusion 41 (2001) 665
-      Fidone, Giruzzi, Granata, Nuclear Fusion 41 (2001) 1755
-
-    Arguments:
-    """
-    # tbet is betaT in Albajar, not to be confused with plasma beta
-
-    tbet = 2.0
-    # rpow is the(1-Rsyn) power dependence based on plasma shape
-    # (see Fidone)
-    rpow = 0.62
-    kap = plasma_volume / (2.0 * 3.1415**2 * rmajor * rminor**2)
-
-    # No account is taken of pedestal profiles here, other than use of
-    # the correct ne0 and te0...
-    de2o = 1.0e-20 * ne0
-    pao = 6.04e3 * (rminor * de2o) / bt
-    gfun = 0.93 * (1.0 + 0.85 * np.exp(-0.82 * rmajor / rminor))
-    kfun = (alphan + 3.87e0 * alphat + 1.46) ** (-0.79)
-    kfun = kfun * (1.98 + alphat) ** 1.36 * tbet**2.14
-    kfun = kfun * (tbet**1.53 + 1.87 * alphat - 0.16) ** (-1.33)
-    dum = 1.0 + 0.12 * (te0 / (pao**0.41)) * (1.0 - ssync) ** 0.41
-    # Very high T modification, from Fidone
-    dum = dum ** (-1.51)
-
-    psync = 3.84e-8 * (1.0e0 - ssync) ** rpow * rmajor * rminor**1.38
-    psync = psync * kap**0.79 * bt**2.62 * de2o**0.38
-    psync = psync * te0 * (16.0 + te0) ** 2.61 * dum * gfun * kfun
-
-    # psyncpv should be per unit volume
-    # Albajar gives it as total
-    psyncpv = psync / plasma_volume
-    print("psyncpv = ", psyncpv * plasma_volume)  # matches the out.dat file
-
-    return psyncpv
-
-
 def plot_radprofile(prof, mfile_data, scan, impp, demo_ranges) -> float:
     """Function to plot radiation profile, formula taken from ???.
 
@@ -1001,22 +963,11 @@ def plot_radprofile(prof, mfile_data, scan, impp, demo_ranges) -> float:
                     1 - min(0.9999, rhopedt)
                 )
 
-        # ncore = neped + (ne0-neped) * (1-rhocore**2/rhopedn**2)**alphan
-        # nsep = nesep + (neped-nesep) * (1-rhosep)/(1-min(0.9999, rhopedn))
-        # ne = np.append(ncore, nsep)
-
-        # The temperatue profile
-        # tcore = teped + (te0-teped) * (1-(rhocore/rhopedt)**tbeta)**alphat
-        # tsep = tesep + (teped-tesep)* (1-rhosep)/(1-min(0.9999,rhopedt))
-        # te = np.append(tcore,tsep)
-
     # Intailise the radiation profile arrays
     pimpden = np.zeros([imp_data.shape[0], te.shape[0]])
     lz = np.zeros([imp_data.shape[0], te.shape[0]])
     prad = np.zeros(te.shape[0])
 
-    # psyncpv = synchrotron_rad()
-
     # Intailise the impurity radiation profile
     for k in range(te.shape[0]):
         for i in range(imp_data.shape[0]):
@@ -1041,18 +992,6 @@ def plot_radprofile(prof, mfile_data, scan, impp, demo_ranges) -> float:
         for l in range(imp_data.shape[0]):  # noqa: E741
             prad[k] = prad[k] + pimpden[l][k] * 2.0e-6
 
-    # benchmark prad again outfile so mod prad
-    # pbremint = (rho[1:] * pbrem[1:]) @ drho
-    # pradint = prad[1:] @ drho * 2.0e-5
-    # pbremint = pbrem[1:] @ drho * 2.0e-5
-
-    # print('prad = ',prad)
-    # print('pbrem = ',pbrem)
-    # print(1.0e32*lz[12])
-    # print('pradpv = ',pradint)
-    # print('pbremmw = ',pbremint*plasma_volume)
-    # print('pradmw = ', pradint*plasma_volume, 'MW') # pimp = pline + pbrem
-
     prof.plot(rho, prad, label="Total")
     prof.plot(rho, pimpden[0] * 2.0e-6, label="H")
     prof.plot(rho, pimpden[1] * 2.0e-6, label="He")
@@ -1438,19 +1377,6 @@ def plot_firstwall(axis, mfile_data, scan, colour_scheme):
         )
 
 
-def angle_check(angle1, angle2):
-    """Function to perform TF coil angle check"""
-    if angle1 > 1:
-        angle1 = 1
-    if angle1 < -1:
-        angle1 = -1
-    if angle2 > 1:
-        angle2 = 1
-    if angle2 < -1:
-        angle2 = -1
-    return angle1, angle2
-
-
 def plot_tf_coils(axis, mfile_data, scan, colour_scheme):
     """Function to plot TF coils
 
@@ -2791,7 +2717,7 @@ def plot_current_drive_info(axis, mfile_data, scan):
             (flh, r"$\frac{P_{\mathrm{div}}}{P_{\mathrm{LH}}}$", ""),
             (hstar, "H* (non-rad. corr.)", ""),
         ]
-        if "iefrffix" in mfile_data.data.keys():
+        if mfile_data.data["iefrffix"].get_scan(scan) != 0:
             data.insert(
                 1, ("pinjmwfix", f"{secondary_heating} secondary auxiliary power", "MW")
             )

diff --git a/process/objectives.py b/process/objectives.py
@@ -0,0 +1,95 @@
+import numpy as np
+
+from process.fortran import (
+    physics_variables,
+    tfcoil_variables,
+    pf_power_variables,
+    current_drive_variables,
+    cost_variables,
+    divertor_variables,
+    times_variables,
+    heat_transport_variables,
+)
+
+
+def objective_function(minmax: int) -> float:
+    """Calculate the specified objective function
+
+    :param minimax: the ID and sign of the figure of merit to evaluate.
+    A negative value indicates maximisation.
+    A positive value indicates minimisation.
+    * 1: Major radius
+    * 3: Neutron wall load
+    * 4: TF coil + PF coil power
+    * 5: Fusion gain
+    * 6: Cost of electricity
+    * 7: Direct/constructed/capital cost
+    * 8: Aspect ratio
+    * 9: Divertor heat load
+    * 10: Toroidal field on axis
+    * 11: Injected power
+    * 14: Pulse length
+    * 15: Plant availability
+    * 16: Major radius/burn time
+    * 17: Net electrical output
+    * 18: NULL, f(x) = 1
+    * 19: Major radius/burn time
+    :type minimax: int
+    """
+    figure_of_merit = abs(minmax)
+
+    # -1 = maximise
+    # +1 = minimise
+    objective_sign = np.sign(minmax)
+
+    match figure_of_merit:
+        case 1:
+            objective_metric = 0.2 * physics_variables.rmajor
+        case 3:
+            objective_metric = physics_variables.wallmw
+        case 4:
+            objective_metric = (
+                tfcoil_variables.tfcmw + 1e-3 * pf_power_variables.srcktpm
+            ) / 10.0
+        case 5:
+            objective_metric = physics_variables.fusion_power / (
+                current_drive_variables.pinjmw
+                + current_drive_variables.porbitlossmw
+                + physics_variables.pohmmw
+            )
+        case 6:
+            objective_metric = cost_variables.coe / 100.0
+        case 7:
+            objective_metric = (
+                cost_variables.cdirt / 1.0e3
+                if cost_variables.ireactor == 0
+                else cost_variables.concost / 1.0e4
+            )
+        case 8:
+            objective_metric = physics_variables.aspect
+        case 9:
+            objective_metric = divertor_variables.hldiv
+        case 10:
+            objective_metric = physics_variables.bt
+        case 11:
+            objective_metric = current_drive_variables.pinjmw
+        case 14:
+            objective_metric = times_variables.t_burn / 2.0e4
+        case 15:
+            if cost_variables.iavail != 1:
+                raise ValueError("minmax=15 requires iavail=1")
+            objective_metric = cost_variables.cfactr
+        case 16:
+            objective_metric = 0.95 * (physics_variables.rmajor / 9.0) - 0.05 * (
+                times_variables.t_burn / 7200.0
+            )
+        case 17:
+            objective_metric = heat_transport_variables.pnetelmw / 500.0
+        case 18:
+            objective_metric = 1.0
+        case 19:
+            objective_metric = -0.5 * (current_drive_variables.bigq / 20.0) - 0.5 * (
+                times_variables.t_burn / 7200.0
+            )
+
+    return objective_sign * objective_metric