diff --git a/CMakeLists.txt b/CMakeLists.txt
index 46979e38c..a2be73910 100755
--- 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
index 2e6bb2efd..925a12452 100644
--- 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
index 8bbe91018..b80d03a9b 100644
--- 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
index df764b977..bd0775d01 100644
--- 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
index 7f9420ebb..ee33c92f9 100755
--- 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
new file mode 100644
index 000000000..e27039cd5
--- /dev/null
+++ 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
diff --git a/process/process_script_advanced.py b/process/process_script_advanced.py
deleted file mode 100644
index 3512ec856..000000000
--- a/process/process_script_advanced.py
+++ /dev/null
@@ -1,242 +0,0 @@
-#!/usr/bin/env python
-"""Script to build and run PROCESS simply.
-
-This compiles, runs and displays output from PROCESS. The aim is to provide a
-common entry point to the code in Python, and to automate a standard workflow.
-"""
-# TODO This requires converting to run the Python-wrapped version of Process
-import subprocess
-import argparse
-from shutil import copy
-from pathlib import Path
-from difflib import unified_diff
-from utilities.process_io_lib import input_validator
-
-# Set required paths outside the Python package directory
-ROOT_DIR = Path(__file__).parent.parent
-# The root project directory (for running cmake)
-PROCESS_EXE_PATH = ROOT_DIR.joinpath("bin/process.exe")
-# The Process executable (for running Process)
-
-
-class Process(object):
- """A Process workflow based on command line arguments."""
-
- def __init__(self):
- """Parse command line arguments and run accordingly."""
- # File paths
- self.run_dir = None
- self.input = None
- self.ref_input = None
-
- # Run actions
- self.parse_args()
- if self.args.build:
- self.build_process()
- self.create_dicts()
-
- # Find the input file and look for changes to it
- self.set_input()
- self.set_ref_input()
- self.input_file_diff()
-
- self.validate_input()
- self.run_process()
- if self.args.util:
- # Only create dicts if necessary for utilities
- self.create_dicts()
- self.run_utils()
-
- def parse_args(self):
- """Parse the command line arguments."""
- parser = argparse.ArgumentParser(
- description="Script to automate " "running PROCESS"
- )
-
- # Optional arguments
- parser.add_argument(
- "--input",
- "-i",
- metavar="input_file_path",
- help="The path to the input file that Process runs on",
- )
- parser.add_argument(
- "--ref_input",
- "-r",
- metavar="reference_input",
- help="Reference input file to record changes against",
- )
- parser.add_argument(
- "--build", "-b", action="store_true", help="Rebuilds Process if present"
- )
- parser.add_argument(
- "--util",
- "-u",
- metavar="util_name",
- help="Utility to run after Process runs",
- )
-
- self.args = parser.parse_args()
- # Store namespace object of the args
-
- def build_process(self):
- """Build Process"""
- # cmake must be run from the project root dir
- subprocess.run(["cmake", "-H.", "-Bbuild"], cwd=ROOT_DIR)
- subprocess.run(["cmake", "--build", "build"], cwd=ROOT_DIR)
-
- def set_input(self):
- """Try to find or validate the input file path, then store it.
-
- Also set the run directory according to the input file path.
- """
- if self.args.input:
- # Input file path specified as CLI argument
- input_path = Path(self.args.input)
- if input_path.exists():
- self.input = input_path
- # Set input as Path object
- else:
- raise FileNotFoundError("Input file not found on this path.")
- else:
- # No input file specified; try to find one in the current dir
- input_files = [
- path
- for path in Path.cwd().iterdir()
- if "IN.DAT" in path.name and "REF_IN.DAT" not in path.name
- ]
-
- if len(input_files) == 0:
- # No input file found
- raise FileNotFoundError("Input file not found in this dir.")
- elif len(input_files) == 1:
- # Only one input file in this dir; use it
- self.input = input_files[0]
- else:
- # More than one input file; ambiguous which one to use
- raise Exception(
- "More than one IN.DAT found in this dir. "
- 'Specifiy which one to use with the "-i" option, or '
- "remove the other ones."
- )
-
- # self.input is now defined
- self.run_dir = self.input.parent
- # Set run directory as dir that contains the input file
-
- def set_ref_input(self):
- """Find an input file to use as a reference for changes.
-
- Find or create a reference input file to compare with the current
- input file, to allow the changes to be tracked.
- """
- # self.args.ref_input: reference input file path; can already be set as
- # command line arg
- if self.args.ref_input:
- # Ref input specified as command line arg. Check it exists
- path = Path(self.args.ref_input)
- if path.exists():
- # Store path object
- self.ref_input = path
- else:
- raise FileNotFoundError(
- "Reference input file not found; " "check the path."
- )
-
- # Input file exists
- if "REF_IN.DAT" in self.ref_input.name:
- # Ref file specified; got a ref file
- pass
- elif "IN.DAT" in self.ref_input.name:
- # Regular input file specified. Make a reference copy as
- # REF_IN.DAT and update self.ref_input
- self.create_ref_input()
- else:
- # File isn't an IN.DAT or a REF_IN.DAT
- raise ValueError(
- "Reference input file should end in IN.DAT or " "REF_IN.DAT"
- )
- else:
- # Ref input not specified: try to find one in the input file's dir
- ref_inputs = list(self.run_dir.glob("*REF_IN.DAT"))
- if len(ref_inputs) == 0:
- # No ref input files: create one from input file
- self.create_ref_input()
- elif len(ref_inputs) == 1:
- # One pre-existing ref input file: set as ref input
- self.ref_input = ref_inputs[0]
- else:
- # More than one ref file: error
- raise Exception(
- "More than one REF_IN.DAT found in this dir. "
- 'Specifiy which one to use with the "-r" option, or '
- "remove the other ones."
- )
-
- # Now either have raised an exception or have a ref input file set,
- # with suffix "REF_IN.DAT"
-
- def create_ref_input(self):
- """Create a reference input file from a pre-existing input file."""
- # First, ensure that self.ref_input is set as an IN.DAT to copy
- if not self.ref_input:
- # No IN.DAT specified to use as ref: use input file instead
- self.ref_input = self.input
-
- # Copy the ref input and save it with the REF_IN.DAT extension
- # Then set this as the new ref input
- new_ref_input_name = self.ref_input.name.replace("IN.DAT", "REF_IN.DAT")
- new_ref_input = Path(self.run_dir).joinpath(new_ref_input_name)
- copy(self.ref_input, new_ref_input)
- self.ref_input = new_ref_input
-
- def input_file_diff(self):
- """Perform a diff between the reference input and input files.
-
- This compares the REF_IN.DAT and IN.DAT files to highlight the
- modifications.
- """
- # Run a diff on the files
- before_file = open(self.ref_input)
- after_file = open(self.input)
- before = before_file.readlines()
- after = after_file.readlines()
- before_file.close()
- after_file.close()
- diff = unified_diff(
- before, after, fromfile=self.ref_input.name, tofile=self.input.name
- )
-
- # Write diff to output file
- diff_path = self.run_dir.joinpath("input.diff")
- with open(diff_path, "w") as diff_file:
- diff_file.writelines(diff)
-
- def create_dicts(self):
- """Create Python dictionaries"""
- # cmake must be run from the project root directory
- subprocess.run(["cmake", "--build", "build", "--target", "dicts"], cwd=ROOT_DIR)
-
- def validate_input(self):
- """Validate the input file using the input_validator module."""
- input_validator.validate(self.input)
-
- def run_process(self):
- """Run Process using the input file path."""
- subprocess.run([PROCESS_EXE_PATH, self.input])
-
- def run_utils(self):
- """Run a utility if specified on the command line."""
- # TODO allow multiple utils to run
- # TODO allow options to be passed to utils
- subprocess.run(["python", "./utilities/" + self.args.util + ".py"])
-
-
-def main():
- """Run Process."""
- print("Running process.py")
- Process()
-
-
-if __name__ == "__main__":
- main()
diff --git a/source/fortran/evaluators.f90 b/source/fortran/evaluators.f90
deleted file mode 100644
index 7a79c03de..000000000
--- a/source/fortran/evaluators.f90
+++ /dev/null
@@ -1,214 +0,0 @@
-! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-module function_evaluator
-
- !! Module containing function evaluators for HYBRD and VMCON
- !! solvers
- !! author: P J Knight, CCFE, Culham Science Centre
- !! N/A
- !! This module contains the function evaluators required
- !! by the two equation solvers in the code.
-
- !! None
- !
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-#ifndef dp
- use, intrinsic :: iso_fortran_env, only: dp=>real64
-#endif
- implicit none
-
- public
-
-contains
-
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-! fcnhyb() is commented out temporarily. It calls the caller() subroutine, which
-! is being moved to Python. fcnhyb() is passed into maths_library.hybrd() as an
-! external subroutine argument, which calls fcnhyb() within a goto block. This
-! is difficult to unravel when converting to Python (so that caller() is called
-! from Python), so it has been decided to comment out hybrd() and fcnhyb()
-! temporarily, disabling the non-optimising solver, to allow Python conversion
-! work using the optimising solver (vmcon()) to continue.
-
-! subroutine fcnhyb(n,xc,rc,iflag)
-
-! !! Function evaluator for EQSOLV
-! !! author: P J Knight, CCFE, Culham Science Centre
-! !! n : input integer : Number of equations and unknowns
-! !! xc(n) : input/output real array : On input XC must contain
-! !! an initial estimate of the solution vector. On output XC
-! !! contains the final estimate of the solution vector.
-! !! rc(n) : output real array : Functions evaluated at the output XC
-! !! iflag : input/output integer : Terminate execution of EQSOLV
-! !! by setting IFLAG to a negative integer.
-! !! This subroutine is the function evaluator for
-! !! EQSOLV (q.v.).
-! !! None
-! !
-! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-! use constraints, only: constraint_eqns
-! use numerics, only: neqns
-! use caller_module, only: caller
-! implicit none
-
-! ! Arguments
-
-! integer, intent(in) :: n
-! real(dp), dimension(n), intent(inout) :: xc
-! real(dp), dimension(n), intent(out) :: rc
-! integer, intent(inout) :: iflag
-
-! ! Local variables
-
-! integer :: ncon, nvars
-
-! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-! nvars = neqns
-! ncon = neqns
-
-! call caller(xc,nvars)
-! call constraint_eqns(ncon,rc,-1)
-
-! ! Set iflag < 0 if program is to be terminated here.
-
-! iflag = 1 * iflag
-
-! end subroutine fcnhyb
-
- subroutine funfom(fc)
-
- !! Objective function evaluator for VMCON
- !! author: P J Knight, CCFE, Culham Science Centre
- !! fc : output real : value of objective function at the output point
- !! This routine evaluates the value of the objective function
- !! i.e. the (normalised) figure-of-merit, at the nvar-dimensional
- !! point of interest.
- !! Each equation for fc gives a value of the
- !! order of unity for the sake of the numerics.
- !!
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- use global_variables, only: xlabel, iscan_global
- use constants, only: nout, iotty, mfile
- use constraints, only: constraint_eqns
- use cost_variables, only: concost, cfactr, cdirt, ireactor, iavail, coe
- use current_drive_variables, only: bigq, porbitlossmw, pinjmw
- use divertor_variables, only: hldiv
- use error_handling, only: idiags, fdiags, errors_on, report_error
- use heat_transport_variables, only: pnetelmw
- use numerics, only: minmax
- use physics_variables, only: fusion_power, bt, rmajor, wallmw, aspect, pohmmw
- use pf_power_variables, only: srcktpm
- use process_output, only: int_to_string3
- use tfcoil_variables, only: tfcmw
- use times_variables, only: t_burn
- implicit none
-
- ! Arguments
-
- real(dp), intent(out) :: fc
-! real(c_double), intent(out) :: fc
- ! Local variables
-
- integer :: iab
- real(dp) :: sgn
-
-! write(*,*) 'Figure of merit 2 (fusion power / input power) is not used.'
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- iab = abs(minmax)
- sgn = sign(1.0D0, real(minmax, kind(1.0D0)))
-
- ! If sgn is -1 the value of fc will be maximised
- ! If sgn is +1 the value of fc will be minimised
-
- select case (iab)
-
- case (1) ! major radius
- fc = sgn * 0.2D0 * rmajor
-
- case (2) ! fusion power / input power
- write(*,*) 'Figure of merit 2 (fusion power / input power) is not used.'
- write(*,*) 'Figure of merit 5 (fusion gain Q) is available.'
- stop 1
- ! fc = sgn * fusion_power / (pinjmw + porbitlossmw + tfcpmw + ppump/1.0D6)
-
- case (3) ! neutron wall load
- fc = sgn * wallmw
-
- case (4) ! TF coil + PF coil power
- fc = sgn * (tfcmw + 1.0D-3*srcktpm)/10.0D0
-
- case (5) ! Q = fusion gain Issue #540
- fc = sgn * fusion_power / (pinjmw + porbitlossmw + pohmmw)
- !fc = sgn * fusion_power / pinjmw
-
- case (6) ! cost of electricity
- fc = sgn * coe/100.0D0
-
- case (7) ! direct/constructed/capital cost
- if (ireactor == 0) then
- fc = sgn * cdirt/1.0D3
- else
- fc = sgn * concost/1.0D4
- end if
-
- case (8) ! aspect ratio
- fc = sgn * aspect
-
- case (9) ! divertor heat load
- fc = sgn * hldiv
-
- case (10) ! toroidal field on axis
- fc = sgn * bt
-
- case (11) ! injection power
- fc = sgn * pinjmw
-
- case (12) ! hydrogen production capital cost
- ! #506 OBSOLETE
- write(*,*) 'Figure of Merit 13 (Hydrogen production) is no longer supported.'
- stop 1
- case (13) ! hydrogen production rate
- ! #506 OBSOLETE
- write(*,*) 'Figure of Merit 13 (Hydrogen production) is no longer supported.'
- stop 1
-
- case (14) ! pulse length
- fc = sgn * t_burn / 2.0D4
-
- case (15) ! plant availability factor (N.B. requires iavail = 1)
-
- if (iavail /= 1) call report_error(23)
-
- fc = sgn * cfactr
-
- case (16) ! major radius/burn time
- fc = sgn * ( 0.95d0 * (rmajor/9.0d0) - 0.05d0 * (t_burn/7200.d0) )
-
- case (17) ! net electrical output
- fc = sgn * pnetelmw / 500.0d0
-
- case (18) ! Null figure of merit
- fc = 1d0
-
- case (19) ! major radius/burn time
- fc = sgn * ( -0.5d0 * (bigq/20.0D0) - 0.5d0 * (t_burn/7200.d0) )
-
- case default
- idiags(1) = iab ; call report_error(24)
-
- end select
-
- ! Crude method of catching NaN errors
-
- if ((abs(fc) > 9.99D99).or.(fc /= fc)) then
- idiags(1) = iab ; call report_error(25)
- end if
-
- end subroutine funfom
-
-end module function_evaluator
diff --git a/source/fortran/scan.f90 b/source/fortran/scan.f90
index ced511e14..397b168e1 100644
--- a/source/fortran/scan.f90
+++ b/source/fortran/scan.f90
@@ -869,7 +869,6 @@ subroutine post_optimise(ifail)
!!
use constraints
use error_handling
- use function_evaluator
use numerics
use process_output
use utilities, only:upper_case
diff --git a/tests/integration/data/test_solver.conf b/tests/integration/data/test_solver.conf
deleted file mode 100644
index a3acd009c..000000000
--- a/tests/integration/data/test_solver.conf
+++ /dev/null
@@ -1,42 +0,0 @@
-* CONFIG FILE FOR SOLVER TEST
-*
-
-* No iterations
-NITER = 1
-
-* sets flag of same name in IN.DAT:
-* None - keeps original value
-* -1 - Non-optimisation only
-* 0 - one non-optimisation step followed by optimised iteration
-* 1 - optimisation only
-* 2 - NAG solver
-IOPTIMZ = 1
-
-* sets the error tolerance for VMCON (default = 1.0D-3)
-* None - keeps the original value (IN.DAT or default)
-EPSVMC = None
-
-* sets the step length for HYBRD (default = 1.0D-3)
-* None - keeps the original value (IN.DAT or default)
-EPSFCN = None
-
-* sets the feasibility tolerance of the NAG solver
-* None - keeps the original value (IN.DAT or default)
-FEASTOL = None
-
-*sets the function precision of the NAG solver
-* None - keeps the original value (IN.DAT or default)
-FUNCPREC = None
-
-*sets the figure of merit
-*positive value - minimise, negative value - maximise
-* None - keeps the original value
-* 1 - plasma major radius
-* 6 - cost of electricity
-MINMAX = None
-
-* integer seed for random number generator; use None for random seed
-SEED = 2
-
-* factor within which the iteration variables are changed
-FACTOR = 1.1
diff --git a/tests/integration/test_solver.py b/tests/integration/test_solver.py
deleted file mode 100755
index a2ccf0de4..000000000
--- a/tests/integration/test_solver.py
+++ /dev/null
@@ -1,207 +0,0 @@
-"""Tests for the Process solver."""
-
-import time
-from numpy import histogram
-from process.io.mfile import MFile
-from process.io.python_fortran_dicts import get_dicts
-from process.io.process_config import TestProcessConfig
-from process.io.process_funcs import (
- get_neqns_itervars,
- update_ixc_bounds,
- get_variable_range,
- vary_iteration_variables,
- process_stopped,
- get_solution_from_mfile,
- process_warnings,
-)
-from process import fortran
-
-# Required to prevent pytest collecting this as a test
-TestProcessConfig.__test__ = False
-
-# Load dicts from dicts JSON file
-process_dicts = get_dicts()
-IFAIL_SUCCESS = process_dicts["IFAIL_SUCCESS"]
-
-
-# TODO: I'm not sure how useful this test is anymore? TN
-def test_solver(temp_data):
- """Code to test PROCESS Solver by choosing different starting values
- for the iteration parameters around the initial values in the INPUT file
- Code modified by Sarah Medley in April 2015 to also calculate Q (fusion power/
- injected power) and output this to SolverTest.out
-
- Author: H. Lux (Hanni.Lux@ccfe.ac.uk)
-
- Date: November 2013 (Initial version)
- March 2014 (Initial release version)
- April 2015 - Code modified by Sarah Medley to also calculate Q - see above
-
- - Input files -
- test_solver.conf: in integration/data dir
-
- - Output files -
- Saved to temporary test dir
- OUT.DAT - PROCESS output of last run
- MFILE.DAT - PROCESS output of last run
- process.log - logfile of PROCESS output to stdout
- time.info - stores the run time of the inner loop
- SolverTest.out -
-
- Ifail values:
- -1: Error in logfile, no solution vector in OUT.DAT
- 0: Error in VMCON: improper input parameters
- 1: normal run
- 2: too many function calls in VMCON
- 3/4: either constraint/objective function + gradients
- are not calculated well or data too noisy
- 5: either no feasible solution to the problem,
- or identity matrix is a bad approximation ot the Hessian
- 6: restriction by an artificial bound or singular matrix
- in quadratic problem.
-
- :param temp_data: temp test dir containing integration test data
- :type temp_data: Path
- """
- # Path to config file and varied input file
- conf_path = temp_data / "test_solver.conf"
- input_path = temp_data / "large_tokamak_IN.DAT"
-
- CONFIG = TestProcessConfig(conf_path)
- CONFIG.setup()
-
- NEQNS, ITERVARS = get_neqns_itervars()
-
- update_ixc_bounds()
-
- fortran.init_module.init_all_module_vars()
- fortran.init_module.init()
-
- LBS, UBS = get_variable_range(ITERVARS, CONFIG.factor)
-
- TABLE_ERR = []
- TABLE_OBJ = []
- TABLE_CON = []
- TABLE_IN = []
- TABLE_OUT = []
- TABLE_POWF = []
- TABLE_PINJ = []
- TABLE_RATIO_POWF_PINJ = []
-
- WARNING_CNT = 0
-
- START_TIME = time.time()
-
- for i in range(CONFIG.niter):
- # modify IN.DAT each time
- input_values = vary_iteration_variables(ITERVARS, LBS, UBS)
-
- TABLE_IN += [input_values]
-
- print(i, end=" ")
-
- # Actually perform the Process run for this varied input file
- CONFIG.run_process(input_path)
-
- if process_stopped():
- ifail, objective_function, constraints, table_sol, table_res = (
- -1,
- "0",
- "0",
- ["0"] * len(ITERVARS),
- ["0"] * NEQNS,
- )
- fusion_power = "0"
- injected_power = "0"
- ratio_fus_inj_power = 0
- else:
- (
- ifail,
- objective_function,
- constraints,
- table_sol,
- table_res,
- ) = get_solution_from_mfile(NEQNS, len(ITERVARS))
-
- # Access MFILE.DAT
- m_file = MFile(filename=str(temp_data / "large_tokamak_MFILE.DAT"))
- ind = -1 # last scan point
-
- # Obtain fusion power and injected power from MFILE.DAT
- fusion_power = m_file.data["fusion_power"].get_scan(ind)
- injected_power = m_file.data["pinjmw"].get_scan(ind)
- ratio_fus_inj_power = fusion_power / injected_power
-
- if ifail == IFAIL_SUCCESS and process_warnings():
- WARNING_CNT += 1
-
- TABLE_OUT += [table_sol + table_res]
- TABLE_ERR += [ifail]
- TABLE_OBJ += [objective_function]
- TABLE_CON += [constraints]
- TABLE_POWF += [fusion_power]
- TABLE_PINJ += [injected_power]
- TABLE_RATIO_POWF_PINJ += [ratio_fus_inj_power]
-
- i += 1
-
- END_TIME = time.time()
-
- ##########
- TIMEFILE = open("time.info", "w")
- TIMEFILE.write("wall clock time of run %f s" % (END_TIME - START_TIME))
- TIMEFILE.close()
-
- #########
- print("\n Error codes:")
- ERRHIST, ERRBINS = histogram(
- TABLE_ERR, bins=[-1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5]
- )
-
- for i in range(len(ERRHIST)):
- print(
- "code %2i: %4i out of %4i = %3.1f "
- % (
- int((ERRBINS[i + 1] + ERRBINS[i]) / 2.0),
- ERRHIST[i],
- CONFIG.niter,
- ERRHIST[i] * 100.0 / CONFIG.niter,
- )
- + "%"
- )
-
- if WARNING_CNT > 0:
- # TODO: calculate fraction of warnings in sucessful runs.
- print("\nIn %f%% of the successful runs warnings occurred." % WARNING_CNT)
-
- ##########
- OUTFILE = open("SolverTest.out", "w")
-
- OUTHEADER = "#Err\t ObjectiveF\t Constraints\t"
- for inp_str in ITERVARS: # input variables
- OUTHEADER += "%s_in\t" % inp_str
- for var_no in range(len(ITERVARS)): # output variables
- OUTHEADER += "itvar%03i\t" % (var_no + 1)
- for con_no in range(NEQNS): # constraints
- OUTHEADER += "constr%03i\t" % (con_no + 1)
- # OUTHEADER += 'Fusion_Power_(MW)\t'
- # OUTHEADER += 'Injected_Power_(MW)\t'
- OUTHEADER += "Ratio_fus_inj_power_(MW)"
- OUTHEADER += "\n"
- OUTFILE.write(OUTHEADER)
-
- for i in range(CONFIG.niter):
- outstr = "%d\t %s\t %s\t" % (TABLE_ERR[i], TABLE_OBJ[i], TABLE_CON[i])
- for valuein in TABLE_IN[i]:
- outstr += "%s\t" % valuein
- for valueout in TABLE_OUT[i]:
- outstr += "%s\t" % valueout
- # outstr += '%d\t' % TABLE_POWF[i]
- # outstr += '%d\t' % TABLE_PINJ[i]
- outstr += "%.1f\t" % TABLE_RATIO_POWF_PINJ[i]
- outstr += "\n"
- OUTFILE.write(outstr)
-
- OUTFILE.close()
-
- # TODO This might need to assert something, like ifail == 1?