From 0e252c37045acb0718b6b04c35e211128ced0841 Mon Sep 17 00:00:00 2001
From: Jonathan Maddock <78556175+jonmaddock@users.noreply.github.com>
Date: Fri, 17 May 2024 12:34:21 +0100
Subject: [PATCH 1/6] Ensure idempotence of objective function and constraints

---
 process/caller.py     | 76 +++++++++++++++++++++++++++++++++++++++++--
 process/evaluators.py | 24 ++++----------
 2 files changed, 81 insertions(+), 19 deletions(-)

diff --git a/process/caller.py b/process/caller.py
index cb8469ca41..397e4e1c95 100644
--- a/process/caller.py
+++ b/process/caller.py
@@ -1,4 +1,9 @@
 from process import fortran as ft
+import numpy as np
+import logging
+from typing import Union, Tuple
+
+logger = logging.getLogger(__name__)
 
 
 class Caller:
@@ -17,9 +22,76 @@ def __init__(self, models, x):
         :type x: np.ndarray
         """
         self.models = models
-        self.call_models(x)
 
-    def call_models(self, xc):
+    @staticmethod
+    def check_agreement(
+        previous: Union[float, np.ndarray], current: Union[float, np.ndarray]
+    ) -> bool:
+        """Compare previous and current arrays for agreement within a tolerance.
+
+        :param previous: value(s) from previous models evaluation
+        :type previous: float | np.ndarray
+        :param current: value(s) from current models evaluation
+        :type current: float | np.ndarray
+        :return: whether values agree or not
+        :rtype: bool
+        """
+        return np.allclose(previous, current, rtol=1.0e-6)
+
+    def call_models(self, xc: np.ndarray, m: int) -> Tuple[float, np.ndarray]:
+        """Evalutate models until results are idempotent.
+
+        Ensure objective function and constraints are idempotent before returning.
+
+        :param xc: optimisation parameters
+        :type xc: np.ndarray
+        :param m: number of constraints
+        :type m: int
+        :raises RuntimeError: if values are non-idempotent after successive
+        evaluations
+        :return: objective function and constraints
+        :rtype: Tuple[float, np.ndarray]
+        """
+        objf_prev = None
+        conf_prev = None
+
+        # Evaluate models up to 10 times; any more implies non-converging values
+        for i in range(10):
+            self._call_models_once(xc)
+            # Evaluate objective function and constraints
+            objf = ft.function_evaluator.funfom()
+            conf, _, _, _, _ = ft.constraints.constraint_eqns(m, -1)
+
+            if objf_prev is None and conf_prev is None:
+                # First run: run again to check idempotence
+                logger.debug("New optimisation parameter vector being evaluated")
+                objf_prev = objf
+                conf_prev = conf
+                continue
+
+            # Check for idempotence
+            if self.check_agreement(objf_prev, objf) and self.check_agreement(
+                conf_prev, conf
+            ):
+                # Idempotent: no longer changing, so return
+                logger.debug(
+                    "Model evaluations idempotent, returning objective "
+                    "function and constraints"
+                )
+                return objf, conf
+            else:
+                # Not idempotent: still changing, so evaluate models again
+                logger.debug("Model evaluations not idempotent: evaluating again")
+                objf_prev = objf
+                conf_prev = conf
+
+        raise RuntimeError(
+            "Model evaluations at the current optimisation parameter vector "
+            "don't produce idempotent values for the objective function and "
+            "constraints."
+        )
+
+    def _call_models_once(self, xc):
         """Call the physics and engineering models.
 
         This method is the principal caller of all the physics and
diff --git a/process/evaluators.py b/process/evaluators.py
index a4d2db18b5..4b3c9d1ff1 100644
--- a/process/evaluators.py
+++ b/process/evaluators.py
@@ -1,12 +1,10 @@
 from process.caller import Caller
 from process.fortran import global_variables as gv
-from process.fortran import constraints
 from process.fortran import cost_variables as cv
 from process.fortran import numerics
 from process.fortran import physics_variables as pv
 from process.fortran import stellarator_variables as sv
 from process.fortran import times_variables as tv
-from process.fortran import function_evaluator
 import numpy as np
 import math
 import logging
@@ -53,15 +51,17 @@ def fcnvmc1(self, n, m, xv, ifail):
         conf = np.zeros(m, dtype=np.float64, order="F")
 
         # Evaluate machine parameters at xv
-        self.caller.call_models(xv)
+        objf, conf = self.caller.call_models(xv, m)
 
         # Convergence loop to ensure burn time consistency
+        # TODO This can be removed once model evaluations become effectively
+        # idempotent
         if sv.istell == 0:
             for _ in range(10):
                 if abs((tv.tburn - tv.tburn0) / max(tv.tburn, 0.01)) <= 0.001:
                     break
 
-                self.caller.call_models(xv)
+                objf, conf = self.caller.call_models(xv, m)
                 if gv.verbose == 1:
                     print("Internal tburn consistency check: ", tv.tburn, tv.tburn0)
             else:
@@ -71,12 +71,6 @@ def fcnvmc1(self, n, m, xv, ifail):
                 )
                 print("tburn, tburn0: ", tv.tburn, tv.tburn0)
 
-        # Evaluate figure of merit (objective function)
-        objf = function_evaluator.funfom()
-
-        # Evaluate constraint equations
-        conf, _, _, _, _ = constraints.constraint_eqns(m, -1)
-
         # Verbose diagnostics
         if gv.verbose == 1:
             summ = 0.0
@@ -141,14 +135,10 @@ def fcnvmc2(self, n, m, xv, lcnorm):
                     xbac[i] = xv[j] * (1.0 - numerics.epsfcn)
 
             # Evaluate at (x+dx)
-            self.caller.call_models(xfor)
-            ffor = function_evaluator.funfom()
-            cfor, _, _, _, _ = constraints.constraint_eqns(m, -1)
+            ffor, cfor = self.caller.call_models(xfor, m)
 
             # Evaluate at (x-dx)
-            self.caller.call_models(xbac)
-            fbac = function_evaluator.funfom()
-            cbac, _, _, _, _ = constraints.constraint_eqns(m, -1)
+            fbac, cbac = self.caller.call_models(xbac, m)
 
             # Calculate finite difference gradients
             fgrd[i] = (ffor - fbac) / (xfor[i] - xbac[i])
@@ -162,6 +152,6 @@ def fcnvmc2(self, n, m, xv, lcnorm):
         # variable in the solution vector is inconsistent with its value
         # shown elsewhere in the output file, which is a factor (1-epsfcn)
         # smaller (i.e. its xbac value above).
-        self.caller.call_models(xv)
+        self.caller.call_models(xv, m)
 
         return fgrd, cnorm

From c2023009cdcd307371afe9b27ead3183776257e4 Mon Sep 17 00:00:00 2001
From: Jonathan Maddock <78556175+jonmaddock@users.noreply.github.com>
Date: Mon, 20 May 2024 10:21:41 +0100
Subject: [PATCH 2/6] Ensure idempotence of final result

Ensure values in mfiles converge before finally writing them.
---
 process/caller.py               | 98 +++++++++++++++++++++++++++++++--
 process/evaluators.py           |  2 +-
 process/final.py                |  5 +-
 process/main.py                 | 31 ++++-------
 process/scan.py                 |  8 +--
 source/fortran/final_module.f90 |  6 +-
 source/fortran/init_module.f90  | 32 +++++++++++
 7 files changed, 152 insertions(+), 30 deletions(-)

diff --git a/process/caller.py b/process/caller.py
index 397e4e1c95..2f3d503201 100644
--- a/process/caller.py
+++ b/process/caller.py
@@ -1,7 +1,12 @@
+from __future__ import annotations
 from process import fortran as ft
 import numpy as np
 import logging
 from typing import Union, Tuple
+import process.main
+from process.final import finalise
+from process.io.mfile import MFile
+from process.utilities.f2py_string_patch import f2py_compatible_to_string
 
 logger = logging.getLogger(__name__)
 
@@ -9,7 +14,7 @@
 class Caller:
     """Calls physics and engineering models."""
 
-    def __init__(self, models, x):
+    def __init__(self, models):
         """Initialise all physics and engineering models.
 
         To ensure that, at the start of a run, all physics/engineering
@@ -18,8 +23,6 @@ def __init__(self, models, x):
 
         :param models: physics and engineering model objects
         :type models: process.main.Models
-        :param x: optimisation parameters
-        :type x: np.ndarray
         """
         self.models = models
 
@@ -56,7 +59,7 @@ def call_models(self, xc: np.ndarray, m: int) -> Tuple[float, np.ndarray]:
         conf_prev = None
 
         # Evaluate models up to 10 times; any more implies non-converging values
-        for i in range(10):
+        for _ in range(10):
             self._call_models_once(xc)
             # Evaluate objective function and constraints
             objf = ft.function_evaluator.funfom()
@@ -91,6 +94,75 @@ def call_models(self, xc: np.ndarray, m: int) -> Tuple[float, np.ndarray]:
             "constraints."
         )
 
+    def call_models_full_idempotence(self, xc: np.ndarray, ifail: int) -> None:
+        """Evaluate models until all results are idempotent.
+
+        Ensure all outputs in mfile are idempotent before returning, by
+        evaluating models multiple times. Typically used at the end of an
+        optimisation, or in a non-optimising evaluation.
+
+        :param xc: optimisation parameter
+        :type xc: np.ndarray
+        :param ifail: return code of solver
+        :type ifail: int
+        :raises RuntimeError: if values are non-idempotent after successive
+        evaluations
+        """
+        # TODO The only way to ensure idempotence in all outputs is by comparing
+        # mfiles at this stage
+        previous_mfile_arr = None
+        file_prefix = ft.global_variables.fileprefix
+        # Evaluate models up to 10 times; any more implies non-converging values
+        for _ in range(10):
+            # Divert OUT.DAT and MFILE.DAT output to scratch files for
+            # idempotence checking
+            ft.init_module.open_idempotence_files()
+            self._call_models_once(xc)
+            # Write mfile
+            finalise(self.models, ifail)
+
+            # Extract mfile data
+            mfile_path = (
+                f2py_compatible_to_string(file_prefix).split("IN.DAT")[0]
+                + "IDEM_MFILE.DAT"
+            )
+            mfile = MFile(mfile_path)
+            mfile_data = {}
+            for var in mfile.data.keys():
+                mfile_data[var] = mfile.data[var].get_scan(-1)
+
+            # Extract floats from mfile dict into array for straightforward
+            # comparison: only compare floats
+            current_mfile_arr = np.array(
+                [val for val in mfile_data.values() if type(val) is float]
+            )
+            if previous_mfile_arr is None:
+                # First run: need another run to compare with
+                logger.debug(
+                    "New mfile created: evaluating models again to check idempotence"
+                )
+                previous_mfile_arr = np.copy(current_mfile_arr)
+                continue
+
+            if self.check_agreement(previous_mfile_arr, current_mfile_arr):
+                # Previous and current mfiles agree (idempotent)
+                logger.debug("Mfiles idempotent, returning")
+                # Divert OUT.DAT and MFILE.DAT output back to original files
+                # now idempotence checking complete
+                ft.init_module.close_idempotence_files()
+                # Write final output file and mfile
+                finalise(self.models, ifail)
+                return
+            else:
+                # Mfiles not yet idempotent: re-evaluate models
+                logger.debug("Mfiles not idempotent, evaluating models again")
+                previous_mfile_arr = np.copy(current_mfile_arr)
+
+        raise RuntimeError(
+            "Model evaluations at the current optimisation parameter vector "
+            "don't produce idempotent values in the final output."
+        )
+
     def _call_models_once(self, xc):
         """Call the physics and engineering models.
 
@@ -227,3 +299,21 @@ def _call_models_once(self, xc):
         self.models.costs.run()
 
         # FISPACT and LOCA model (not used)- removed
+
+
+def write_idempotent_result(models: process.main.Models, ifail: int) -> None:
+    """Ensure result idempotence, then write output files.
+
+    :param models: physics and engineering models
+    :type models: process.main.Models
+    :param ifail: solver return code
+    :type ifail: int
+    """
+    n = ft.numerics.nvar
+    x = ft.numerics.xcm[:n]
+    # Call models, ensuring output mfiles are fully idempotent
+    caller = Caller(models)
+    caller.call_models_full_idempotence(
+        xc=x,
+        ifail=ifail,
+    )
diff --git a/process/evaluators.py b/process/evaluators.py
index 4b3c9d1ff1..d1bfd3532d 100644
--- a/process/evaluators.py
+++ b/process/evaluators.py
@@ -23,7 +23,7 @@ def __init__(self, models, x):
         :param x: optimisation parameters
         :type x: np.ndarray
         """
-        self.caller = Caller(models, x)
+        self.caller = Caller(models)
 
     def fcnvmc1(self, n, m, xv, ifail):
         """Function evaluator for VMCON.
diff --git a/process/final.py b/process/final.py
index 94a61e6502..06a31142fb 100644
--- a/process/final.py
+++ b/process/final.py
@@ -1,4 +1,5 @@
 """Final output at the end of a scan."""
+
 from process import fortran as ft
 from process.fortran import final_module as fm
 from process import output as op
@@ -7,6 +8,8 @@
 def finalise(models, ifail):
     """Routine to print out the final point in the scan.
 
+    Writes to OUT.DAT and MFILE.DAT.
+
     :param models: physics and engineering model objects
     :type models: process.main.Models
     :param ifail: error flag
@@ -18,7 +21,7 @@ def finalise(models, ifail):
     if ft.numerics.ioptimz == -2:
         fm.no_optimisation()
 
-    # Write output to OUT.DAT
+    # Write output to OUT.DAT and MFILE.DAT
     op.write(models, ft.constants.nout)
 
     fm.final_output()
diff --git a/process/main.py b/process/main.py
index 3d6859988d..351e0dc6f8 100644
--- a/process/main.py
+++ b/process/main.py
@@ -49,7 +49,6 @@
 from process.plasma_geometry import PlasmaGeom
 from process.pulse import Pulse
 from process.scan import Scan
-from process import final
 from process.stellarator import Stellarator
 from process.structure import Structure
 from process.build import Build
@@ -61,7 +60,6 @@
 from process.availability import Availability
 from process.ife import IFE
 from process.costs_2015 import Costs2015
-from process.caller import Caller
 from process.power import Power
 from process.cs_fatigue import CsFatigue
 from process.physics import Physics
@@ -73,6 +71,8 @@
 from process.fw import Fw
 from process.current_drive import CurrentDrive
 from process.impurity_radiation import initialise_imprad
+from process.caller import write_idempotent_result
+
 
 from pathlib import Path
 import os
@@ -459,25 +459,18 @@ def run_scan(self, solver):
         :param solver: which solver to use, as specified in solver.py
         :type solver: str
         """
-        if fortran.numerics.ioptimz >= 0:
+        if fortran.numerics.ioptimz == 1:
+            # VMCON optimisation
             self.scan = Scan(self.models, solver)
-        else:
-            # If no optimisation will be done, compute the OP variables now
-            if fortran.numerics.ioptimz == -2:
-                # Get optimisation parameters x, perform first evaluation of models
-                fortran.define_iteration_variables.loadxc()
-                n = fortran.numerics.nvar
-                x = fortran.numerics.xcm[:n]
-                caller = Caller(self.models, x)
-
-                # To ensure that, at the start of a run, all physics/engineering
-                # variables are fully initialised with consistent values, we perform
-                # a second evaluation call here
-                caller.call_models(x)
-                self.ifail = 6
-
-            final.finalise(self.models, self.ifail)
+        elif fortran.numerics.ioptimz == -2:
+            # No optimisation: compute the output variables now
+            # Get optimisation parameters x, evaluate models
+            fortran.define_iteration_variables.loadxc()
+            self.ifail = 6
+            write_idempotent_result(models=self.models, ifail=self.ifail)
             self.show_errors()
+        else:
+            raise ValueError(f"Invalid ioptimz value: {fortran.numerics.ioptimz}")
 
     def show_errors(self):
         """Report all informational/error messages encountered."""
diff --git a/process/scan.py b/process/scan.py
index 664c4eff48..a40eb533ce 100644
--- a/process/scan.py
+++ b/process/scan.py
@@ -2,8 +2,8 @@
 from process.fortran import scan_module
 from process.fortran import numerics
 from process.optimiser import Optimiser
-from process import final
 import numpy as np
+from process.caller import write_idempotent_result
 
 
 class Scan:
@@ -34,7 +34,7 @@ def run_scan(self):
 
         if scan_module.isweep == 0:
             ifail = self.doopt()
-            final.finalise(self.models, ifail)
+            write_idempotent_result(models=self.models, ifail=ifail)
             error_handling.show_errors()
             return
 
@@ -76,7 +76,7 @@ def scan_1d(self):
             scan_module.scan_1d_write_point_header(iscan)
             ifail = self.doopt()
             scan_1d_ifail_dict[iscan] = ifail
-            final.finalise(self.models, ifail)
+            write_idempotent_result(models=self.models, ifail=ifail)
 
             # outvar is an intent(out) of scan_1d_store_output()
             outvar = scan_module.scan_1d_store_output(
@@ -146,7 +146,7 @@ def scan_2d(self):
                 )
                 ifail = self.doopt()
 
-                final.finalise(self.models, ifail)
+                write_idempotent_result(models=self.models, ifail=ifail)
 
                 outvar, sweep_1_vals, sweep_2_vals = scan_module.scan_2d_store_output(
                     ifail,
diff --git a/source/fortran/final_module.f90 b/source/fortran/final_module.f90
index b67fc728e4..06ae0c2680 100644
--- a/source/fortran/final_module.f90
+++ b/source/fortran/final_module.f90
@@ -73,7 +73,7 @@ subroutine no_optimisation()
   end subroutine no_optimisation
 
   subroutine final_output()
-    use constants, only: iotty
+    use constants, only: iotty, mfile
     use numerics, only: nfev1, ncalls, xcm, ioptimz, nviter, &
       nvar
     use define_iteration_variables, only: loadxc
@@ -117,5 +117,9 @@ subroutine final_output()
           t2,'The HYBRD point required ',i5,' iterations',/, &
           t2,'The optimisation required ',i5,' iterations',/, &
           t2,'There were ',i6,' function calls')
+
+    ! Required to ensure mfile fully written before any model evaluation
+    ! idempotence checks
+    flush(mfile)
   end subroutine final_output
 end module final_module
diff --git a/source/fortran/init_module.f90 b/source/fortran/init_module.f90
index 008332c9ca..2e1f6e1a66 100644
--- a/source/fortran/init_module.f90
+++ b/source/fortran/init_module.f90
@@ -168,6 +168,38 @@ subroutine init
 
    end subroutine init
 
+   subroutine open_idempotence_files
+      ! Open new output file and mfile to write output to
+      ! This is used when checking model evaluation idempotence, to avoid
+      ! polluting the final output file and mfile with intermediate result checks
+      use global_variables, only: output_prefix
+      use constants, only: nout, mfile
+      implicit none
+
+      ! Close existing output file and mfile (could be original out and mfiles
+      ! or idem scratch files)
+      close(unit = nout)
+      close(unit = mfile)
+      ! Open scratch files with same units
+      open(unit=nout, file=trim(output_prefix)//'IDEM_OUT.DAT', action='write', status='replace')
+      open(unit=mfile, file=trim(output_prefix)//'IDEM_MFILE.DAT', action='write', status='replace')
+   end subroutine open_idempotence_files
+
+   subroutine close_idempotence_files
+      ! Revert back to writing to original OUT.DAT and MFILE.DAT, after
+      ! checking model evaluation idempotence
+      use global_variables, only: output_prefix
+      use constants, only: nout, mfile
+      implicit none
+
+      ! Close idempotence files, deleting them in the process
+      close(unit = nout, status="delete")
+      close(unit = mfile, status="delete")
+      ! Re-open original output file and mfile, appending future output to them
+      open(unit=nout, file=trim(output_prefix)//'OUT.DAT', action='write', position='append')
+      open(unit=mfile, file=trim(output_prefix)//'MFILE.DAT', action='write', position='append')
+   end subroutine close_idempotence_files
+
    subroutine finish
       ! Originally at the end of the "program", this subroutine writes some final
       ! lines via the output module and then closes any open files. This is

From f34bc617272c900b274ceb4be80bed3a3eba62e5 Mon Sep 17 00:00:00 2001
From: Jonathan Maddock <78556175+jonmaddock@users.noreply.github.com>
Date: Thu, 23 May 2024 09:17:01 +0100
Subject: [PATCH 3/6] Improve exception messages

Remove superfluous else statements

Fix type hints in caller.py

Prevent circular import of process.main.

Use output file prefix for idempotence mfile checks

Improve naming and docstrings
---
 process/caller.py | 61 +++++++++++++++++++++++++----------------------
 process/main.py   |  9 ++++---
 process/scan.py   |  8 +++----
 3 files changed, 42 insertions(+), 36 deletions(-)

diff --git a/process/caller.py b/process/caller.py
index 2f3d503201..621e53791f 100644
--- a/process/caller.py
+++ b/process/caller.py
@@ -2,11 +2,13 @@
 from process import fortran as ft
 import numpy as np
 import logging
-from typing import Union, Tuple
-import process.main
 from process.final import finalise
 from process.io.mfile import MFile
 from process.utilities.f2py_string_patch import f2py_compatible_to_string
+from typing import Union, Tuple, TYPE_CHECKING
+
+if TYPE_CHECKING:
+    from process.main import Models
 
 logger = logging.getLogger(__name__)
 
@@ -14,7 +16,7 @@
 class Caller:
     """Calls physics and engineering models."""
 
-    def __init__(self, models):
+    def __init__(self, models: Models) -> None:
         """Initialise all physics and engineering models.
 
         To ensure that, at the start of a run, all physics/engineering
@@ -22,7 +24,7 @@ def __init__(self, models):
         called with the initial optimisation paramters, x.
 
         :param models: physics and engineering model objects
-        :type models: process.main.Models
+        :type models: Models
         """
         self.models = models
 
@@ -82,24 +84,25 @@ def call_models(self, xc: np.ndarray, m: int) -> Tuple[float, np.ndarray]:
                     "function and constraints"
                 )
                 return objf, conf
-            else:
-                # Not idempotent: still changing, so evaluate models again
-                logger.debug("Model evaluations not idempotent: evaluating again")
-                objf_prev = objf
-                conf_prev = conf
+
+            # Not idempotent: still changing, so evaluate models again
+            logger.debug("Model evaluations not idempotent: evaluating again")
+            objf_prev = objf
+            conf_prev = conf
 
         raise RuntimeError(
-            "Model evaluations at the current optimisation parameter vector "
-            "don't produce idempotent values for the objective function and "
-            "constraints."
+            "After 10 model evaluations at the current optimisation parameter "
+            "vector, values for the objective function and constraints haven't "
+            "converged (don't produce idempotent values)."
         )
 
-    def call_models_full_idempotence(self, xc: np.ndarray, ifail: int) -> None:
-        """Evaluate models until all results are idempotent.
+    def call_models_and_write_output(self, xc: np.ndarray, ifail: int) -> None:
+        """Evaluate models until results are idempotent, then write output files.
 
         Ensure all outputs in mfile are idempotent before returning, by
         evaluating models multiple times. Typically used at the end of an
-        optimisation, or in a non-optimising evaluation.
+        optimisation, or in a non-optimising evaluation. Writes OUT.DAT and
+        MFILE.DAT with final results.
 
         :param xc: optimisation parameter
         :type xc: np.ndarray
@@ -111,7 +114,10 @@ def call_models_full_idempotence(self, xc: np.ndarray, ifail: int) -> None:
         # TODO The only way to ensure idempotence in all outputs is by comparing
         # mfiles at this stage
         previous_mfile_arr = None
-        file_prefix = ft.global_variables.fileprefix
+        # Path of output file prefix required for writing intermediate idempotence-
+        # checking mfiles in the same directory as input file
+        output_prefix = ft.global_variables.output_prefix
+
         # Evaluate models up to 10 times; any more implies non-converging values
         for _ in range(10):
             # Divert OUT.DAT and MFILE.DAT output to scratch files for
@@ -122,10 +128,7 @@ def call_models_full_idempotence(self, xc: np.ndarray, ifail: int) -> None:
             finalise(self.models, ifail)
 
             # Extract mfile data
-            mfile_path = (
-                f2py_compatible_to_string(file_prefix).split("IN.DAT")[0]
-                + "IDEM_MFILE.DAT"
-            )
+            mfile_path = f2py_compatible_to_string(output_prefix) + "IDEM_MFILE.DAT"
             mfile = MFile(mfile_path)
             mfile_data = {}
             for var in mfile.data.keys():
@@ -153,17 +156,17 @@ def call_models_full_idempotence(self, xc: np.ndarray, ifail: int) -> None:
                 # Write final output file and mfile
                 finalise(self.models, ifail)
                 return
-            else:
-                # Mfiles not yet idempotent: re-evaluate models
-                logger.debug("Mfiles not idempotent, evaluating models again")
-                previous_mfile_arr = np.copy(current_mfile_arr)
+
+            # Mfiles not yet idempotent: re-evaluate models
+            logger.debug("Mfiles not idempotent, evaluating models again")
+            previous_mfile_arr = np.copy(current_mfile_arr)
 
         raise RuntimeError(
             "Model evaluations at the current optimisation parameter vector "
             "don't produce idempotent values in the final output."
         )
 
-    def _call_models_once(self, xc):
+    def _call_models_once(self, xc: np.ndarray) -> None:
         """Call the physics and engineering models.
 
         This method is the principal caller of all the physics and
@@ -301,11 +304,11 @@ def _call_models_once(self, xc):
         # FISPACT and LOCA model (not used)- removed
 
 
-def write_idempotent_result(models: process.main.Models, ifail: int) -> None:
-    """Ensure result idempotence, then write output files.
+def write_output_files(models: Models, ifail: int) -> None:
+    """Evaluate models and write output files (OUT.DAT and MFILE.DAT).
 
     :param models: physics and engineering models
-    :type models: process.main.Models
+    :type models: Models
     :param ifail: solver return code
     :type ifail: int
     """
@@ -313,7 +316,7 @@ def write_idempotent_result(models: process.main.Models, ifail: int) -> None:
     x = ft.numerics.xcm[:n]
     # Call models, ensuring output mfiles are fully idempotent
     caller = Caller(models)
-    caller.call_models_full_idempotence(
+    caller.call_models_and_write_output(
         xc=x,
         ifail=ifail,
     )
diff --git a/process/main.py b/process/main.py
index 351e0dc6f8..ef91d5e740 100644
--- a/process/main.py
+++ b/process/main.py
@@ -71,7 +71,7 @@
 from process.fw import Fw
 from process.current_drive import CurrentDrive
 from process.impurity_radiation import initialise_imprad
-from process.caller import write_idempotent_result
+from process.caller import write_output_files
 
 
 from pathlib import Path
@@ -467,10 +467,13 @@ def run_scan(self, solver):
             # Get optimisation parameters x, evaluate models
             fortran.define_iteration_variables.loadxc()
             self.ifail = 6
-            write_idempotent_result(models=self.models, ifail=self.ifail)
+            write_output_files(models=self.models, ifail=self.ifail)
             self.show_errors()
         else:
-            raise ValueError(f"Invalid ioptimz value: {fortran.numerics.ioptimz}")
+            raise ValueError(
+                f"Invalid ioptimz value: {fortran.numerics.ioptimz}. Please "
+                "select either 1 (optimise) or -2 (no optimisation)."
+            )
 
     def show_errors(self):
         """Report all informational/error messages encountered."""
diff --git a/process/scan.py b/process/scan.py
index a40eb533ce..56689fd737 100644
--- a/process/scan.py
+++ b/process/scan.py
@@ -3,7 +3,7 @@
 from process.fortran import numerics
 from process.optimiser import Optimiser
 import numpy as np
-from process.caller import write_idempotent_result
+from process.caller import write_output_files
 
 
 class Scan:
@@ -34,7 +34,7 @@ def run_scan(self):
 
         if scan_module.isweep == 0:
             ifail = self.doopt()
-            write_idempotent_result(models=self.models, ifail=ifail)
+            write_output_files(models=self.models, ifail=ifail)
             error_handling.show_errors()
             return
 
@@ -76,7 +76,7 @@ def scan_1d(self):
             scan_module.scan_1d_write_point_header(iscan)
             ifail = self.doopt()
             scan_1d_ifail_dict[iscan] = ifail
-            write_idempotent_result(models=self.models, ifail=ifail)
+            write_output_files(models=self.models, ifail=ifail)
 
             # outvar is an intent(out) of scan_1d_store_output()
             outvar = scan_module.scan_1d_store_output(
@@ -146,7 +146,7 @@ def scan_2d(self):
                 )
                 ifail = self.doopt()
 
-                write_idempotent_result(models=self.models, ifail=ifail)
+                write_output_files(models=self.models, ifail=ifail)
 
                 outvar, sweep_1_vals, sweep_2_vals = scan_module.scan_2d_store_output(
                     ifail,

From 35229159f941bb71cc71b4ddde1339c9e2a6dae8 Mon Sep 17 00:00:00 2001
From: Jonathan Maddock <78556175+jonmaddock@users.noreply.github.com>
Date: Thu, 23 May 2024 14:48:18 +0100
Subject: [PATCH 4/6] Remove burn time consistency check

---
 process/evaluators.py | 19 -------------------
 1 file changed, 19 deletions(-)

diff --git a/process/evaluators.py b/process/evaluators.py
index d1bfd3532d..3d3ad3f9fe 100644
--- a/process/evaluators.py
+++ b/process/evaluators.py
@@ -3,7 +3,6 @@
 from process.fortran import cost_variables as cv
 from process.fortran import numerics
 from process.fortran import physics_variables as pv
-from process.fortran import stellarator_variables as sv
 from process.fortran import times_variables as tv
 import numpy as np
 import math
@@ -53,24 +52,6 @@ def fcnvmc1(self, n, m, xv, ifail):
         # Evaluate machine parameters at xv
         objf, conf = self.caller.call_models(xv, m)
 
-        # Convergence loop to ensure burn time consistency
-        # TODO This can be removed once model evaluations become effectively
-        # idempotent
-        if sv.istell == 0:
-            for _ in range(10):
-                if abs((tv.tburn - tv.tburn0) / max(tv.tburn, 0.01)) <= 0.001:
-                    break
-
-                objf, conf = self.caller.call_models(xv, m)
-                if gv.verbose == 1:
-                    print("Internal tburn consistency check: ", tv.tburn, tv.tburn0)
-            else:
-                print(
-                    "Burn time values are not consistent in iteration: ",
-                    numerics.nviter,
-                )
-                print("tburn, tburn0: ", tv.tburn, tv.tburn0)
-
         # Verbose diagnostics
         if gv.verbose == 1:
             summ = 0.0

From 69ec995f5619c72f23f7f9aa75589a2b6341a994 Mon Sep 17 00:00:00 2001
From: Jonathan Maddock <78556175+jonmaddock@users.noreply.github.com>
Date: Fri, 31 May 2024 16:59:02 +0100
Subject: [PATCH 5/6] Remove unnecessary variable in caller.py

Improve docstring in close_idempotence_files()
---
 process/caller.py              | 10 +++++-----
 source/fortran/init_module.f90 |  5 +++--
 2 files changed, 8 insertions(+), 7 deletions(-)

diff --git a/process/caller.py b/process/caller.py
index 621e53791f..a28f5c0aa6 100644
--- a/process/caller.py
+++ b/process/caller.py
@@ -114,9 +114,6 @@ def call_models_and_write_output(self, xc: np.ndarray, ifail: int) -> None:
         # TODO The only way to ensure idempotence in all outputs is by comparing
         # mfiles at this stage
         previous_mfile_arr = None
-        # Path of output file prefix required for writing intermediate idempotence-
-        # checking mfiles in the same directory as input file
-        output_prefix = ft.global_variables.output_prefix
 
         # Evaluate models up to 10 times; any more implies non-converging values
         for _ in range(10):
@@ -127,8 +124,11 @@ def call_models_and_write_output(self, xc: np.ndarray, ifail: int) -> None:
             # Write mfile
             finalise(self.models, ifail)
 
-            # Extract mfile data
-            mfile_path = f2py_compatible_to_string(output_prefix) + "IDEM_MFILE.DAT"
+            # Extract data from intermediate idempotence-checking mfile
+            mfile_path = (
+                f2py_compatible_to_string(ft.global_variables.output_prefix)
+                + "IDEM_MFILE.DAT"
+            )
             mfile = MFile(mfile_path)
             mfile_data = {}
             for var in mfile.data.keys():
diff --git a/source/fortran/init_module.f90 b/source/fortran/init_module.f90
index 2e1f6e1a66..18d80364be 100644
--- a/source/fortran/init_module.f90
+++ b/source/fortran/init_module.f90
@@ -186,8 +186,9 @@ subroutine open_idempotence_files
    end subroutine open_idempotence_files
 
    subroutine close_idempotence_files
-      ! Revert back to writing to original OUT.DAT and MFILE.DAT, after
-      ! checking model evaluation idempotence
+      ! Close the intermediate idempotence-check files, deleting them in the process
+      ! Re-open the original OUT.DAT and MFILE.DAT output files, ready to write
+      ! the final data, now model evaluation idempotence has been checked
       use global_variables, only: output_prefix
       use constants, only: nout, mfile
       implicit none

From e4b35b89d7dc9b8adcc487f8bca8c7fea7ceef10 Mon Sep 17 00:00:00 2001
From: Jonathan Maddock <78556175+jonmaddock@users.noreply.github.com>
Date: Fri, 31 May 2024 17:10:37 +0100
Subject: [PATCH 6/6] Ensure idempotence files are deleted if exception raised

Improve type check in caller.py
---
 process/caller.py | 98 +++++++++++++++++++++++++----------------------
 1 file changed, 52 insertions(+), 46 deletions(-)

diff --git a/process/caller.py b/process/caller.py
index a28f5c0aa6..6f04225dee 100644
--- a/process/caller.py
+++ b/process/caller.py
@@ -115,56 +115,62 @@ def call_models_and_write_output(self, xc: np.ndarray, ifail: int) -> None:
         # mfiles at this stage
         previous_mfile_arr = None
 
-        # Evaluate models up to 10 times; any more implies non-converging values
-        for _ in range(10):
-            # Divert OUT.DAT and MFILE.DAT output to scratch files for
-            # idempotence checking
-            ft.init_module.open_idempotence_files()
-            self._call_models_once(xc)
-            # Write mfile
-            finalise(self.models, ifail)
+        try:
+            # Evaluate models up to 10 times; any more implies non-converging values
+            for _ in range(10):
+                # Divert OUT.DAT and MFILE.DAT output to scratch files for
+                # idempotence checking
+                ft.init_module.open_idempotence_files()
+                self._call_models_once(xc)
+                # Write mfile
+                finalise(self.models, ifail)
 
-            # Extract data from intermediate idempotence-checking mfile
-            mfile_path = (
-                f2py_compatible_to_string(ft.global_variables.output_prefix)
-                + "IDEM_MFILE.DAT"
-            )
-            mfile = MFile(mfile_path)
-            mfile_data = {}
-            for var in mfile.data.keys():
-                mfile_data[var] = mfile.data[var].get_scan(-1)
-
-            # Extract floats from mfile dict into array for straightforward
-            # comparison: only compare floats
-            current_mfile_arr = np.array(
-                [val for val in mfile_data.values() if type(val) is float]
-            )
-            if previous_mfile_arr is None:
-                # First run: need another run to compare with
-                logger.debug(
-                    "New mfile created: evaluating models again to check idempotence"
+                # Extract data from intermediate idempotence-checking mfile
+                mfile_path = (
+                    f2py_compatible_to_string(ft.global_variables.output_prefix)
+                    + "IDEM_MFILE.DAT"
+                )
+                mfile = MFile(mfile_path)
+                mfile_data = {}
+                for var in mfile.data.keys():
+                    mfile_data[var] = mfile.data[var].get_scan(-1)
+
+                # Extract floats from mfile dict into array for straightforward
+                # comparison: only compare floats
+                current_mfile_arr = np.array(
+                    [val for val in mfile_data.values() if isinstance(val, float)]
                 )
+                if previous_mfile_arr is None:
+                    # First run: need another run to compare with
+                    logger.debug(
+                        "New mfile created: evaluating models again to check idempotence"
+                    )
+                    previous_mfile_arr = np.copy(current_mfile_arr)
+                    continue
+
+                if self.check_agreement(previous_mfile_arr, current_mfile_arr):
+                    # Previous and current mfiles agree (idempotent)
+                    logger.debug("Mfiles idempotent, returning")
+                    # Divert OUT.DAT and MFILE.DAT output back to original files
+                    # now idempotence checking complete
+                    ft.init_module.close_idempotence_files()
+                    # Write final output file and mfile
+                    finalise(self.models, ifail)
+                    return
+
+                # Mfiles not yet idempotent: re-evaluate models
+                logger.debug("Mfiles not idempotent, evaluating models again")
                 previous_mfile_arr = np.copy(current_mfile_arr)
-                continue
-
-            if self.check_agreement(previous_mfile_arr, current_mfile_arr):
-                # Previous and current mfiles agree (idempotent)
-                logger.debug("Mfiles idempotent, returning")
-                # Divert OUT.DAT and MFILE.DAT output back to original files
-                # now idempotence checking complete
-                ft.init_module.close_idempotence_files()
-                # Write final output file and mfile
-                finalise(self.models, ifail)
-                return
-
-            # Mfiles not yet idempotent: re-evaluate models
-            logger.debug("Mfiles not idempotent, evaluating models again")
-            previous_mfile_arr = np.copy(current_mfile_arr)
 
-        raise RuntimeError(
-            "Model evaluations at the current optimisation parameter vector "
-            "don't produce idempotent values in the final output."
-        )
+            raise RuntimeError(
+                "Model evaluations at the current optimisation parameter vector "
+                "don't produce idempotent values in the final output."
+            )
+        except Exception:
+            # If exception in model evaluations or idempotence can't be
+            # achieved, delete intermediate idempotence files to clean up
+            ft.init_module.close_idempotence_files()
+            raise
 
     def _call_models_once(self, xc: np.ndarray) -> None:
         """Call the physics and engineering models.