Impurity radiation corrections

process/core/caller.py

@@ -206,9 +206,28 @@ def call_models_and_write_output(self, xc: np.ndarray, ifail: int):
                     OutputFileManager.close_idempotence_files(
                         self.data.globals.output_prefix
                     )
+
+                    # Now idempotent, return
+                    # Pass model caller and opt params for stability constraint evaluation
+                    con_residuals_normalised, _, con_residuals, _, _ = (
+                        constraints.constraint_eqns(
+                            self.data.numerics.neqns + self.data.numerics.nineqns,
+                            -1,
+                            self.data,
+                        )
+                    )
+                    # Evaluate constraints and store in numerics during solver iterations:
+                    # can be used in objective function 20. Hence evaluate before objective
+                    # calculated
+                    self.data.numerics.constraint_residuals_normalised = (
+                        con_residuals_normalised
+                    )
+                    self.data.numerics.constraint_residuals = con_residuals
+                    # Evaluate objective function and constraints
+                    objf = objective_function(self.data.numerics.minmax, self.data)
                     # Write final output file and mfile
                     finalise(self.models, self.data, ifail)
-                    return
+                    return objf, con_residuals_normalised
 
                 # Mfiles not yet idempotent: need to re-evaluate models
                 logger.debug("Mfiles not idempotent, evaluating models again")
@@ -419,7 +438,7 @@ def finalise(models, data, ifail: int, non_idempotent_msg: str | None = None):
         po.oheadr(constants.NOUT, "Final UNFEASIBLE Point")
 
     # Output relevant to no optimisation
-    if data.numerics.ioptimz == PROCESSRunMode.EVALUATION:
+    if data.numerics.ioptimz in {PROCESSRunMode.EVALUATION, PROCESSRunMode.SOLUTION}:
         output_evaluation(data)
 
     # Print non-idempotence warning to OUT.DAT only

process/data_structure/numerics.py

@@ -49,13 +49,15 @@ class PROCESSRunMode(IntEnum):
     """In this mode, the code will not perform any optimisation, and will instead
     simply evaluate the constraints for the given input parameters, which is useful
     for testing and for evaluating the performance of a given design point without
-    trying to optimise it. Internally, PROCESS uses `fsolve` (a Newton-Krylov/hybrd
+    trying to optimise it.
+    """
+    SOLUTION = (-1, "Solution mode (no optimisation)")
+    """Internally, PROCESS uses `fsolve` (a Newton-Krylov/hybrd
     root-finding method from `scipy.optimize`) to seek a *consistent* solution by
     varying a subset of the iteration variables until the consistency constraints
     (equality constraints whose residuals must be driven to zero) are simultaneously
     satisfied; no figure-of-merit is optimised, and the solver simply tries to find
-    a root of the constraint-residual vector.
-    """
+    a root of the constraint-residual vector."""
     OPTIMISATION = (1, "Optimisation mode (e.g. via VMCON)")
     """In this mode, the code will perform optimisation using the VMCON solver
     (or a custom solver if specified) to try to find a design point that optimises
@@ -182,6 +184,14 @@ class NumericsData:
     nvar: int = 0
     """number of iteration variables to use"""
 
+    # Constraint residuals, updated on every iteration
+    constraint_residuals_normalised: list[float] = field(
+        default_factory=lambda: np.array([0] * IPEQNS)
+    )
+    constraint_residuals: list[float] = field(
+        default_factory=lambda: np.array([0] * IPEQNS)
+    )
+
     nviter: int = 0
     """number of optimisation iterations performed"""
 

process/data_structure/physics_variables.py

@@ -1078,6 +1078,9 @@ class PhysicsData:
     pden_plasma_core_rad_mw: float = 0.0
     """total core radiation power per volume (MW/m3)"""
 
+    pden_plasma_core_rad_tauE_mw: float = 0.0
+    """reduced total core radiation power for tauE calculation (MW/m3)"""
+
     p_dd_total_mw: float = 0.0
     """deuterium-deuterium fusion power (MW)"""
 

process/main.py

@@ -119,6 +119,8 @@
 )
 from process.models.vacuum import Vacuum, VacuumVessel
 from process.models.water_use import WaterUse
+from process.core.caller import write_output_files
+from process.core.solver.iteration_variables import load_iteration_variables
 
 PACKAGE_LOGGING = True
 """Can be set False to disable package-level logging, e.g. in the test suite"""
@@ -449,10 +451,19 @@ def run_scan(self):
         # ioptimz == 1: optimisation
         if self.data.numerics.ioptimz == PROCESSRunMode.OPTIMISATION:
             pass
-        elif self.data.numerics.ioptimz == PROCESSRunMode.EVALUATION:
-            # No optimisation:
-            # solve equality (consistency) constraints only using fsolve (HYBRD)
+        # ioptimz == -1: solution
+        elif self.data.numerics.ioptimz == PROCESSRunMode.SOLUTION:
+            # Solve equality (consistency) constraints only using fsolve (HYBRD)
             self.solver = "fsolve"
+        # ioptimz == -2: evaluation
+        elif self.data.numerics.ioptimz == PROCESSRunMode.EVALUATION:
+            # Evalutation only: compute the output variables now
+            # Get optimisation parameters x, evaluate models
+            load_iteration_variables(self.data)
+            self.ifail = 6
+            write_output_files(data=self.data, models=self.models, ifail=self.ifail)
+            self.show_errors()
+            return
         else:
             raise ValueError(
                 f"Invalid ioptimz value: {self.data.numerics.ioptimz}. Please "

process/models/costs/costs.py

@@ -2854,7 +2854,8 @@ def coelc(self):
 
         #  Capital recovery factor
 
-        crfcdr = (fefcdr * self.data.costs.discount_rate) / (fefcdr - 1.0e0)
+        # crfcdr = (fefcdr * self.data.costs.discount_rate) / (fefcdr - 1.0e0)
+        crfcdr = 1.0
 
         #  Annual cost of replacements
 

process/models/physics/confinement_time.py

@@ -18,6 +18,7 @@
     PlasmaIgnitionModel,
 )
 from process.models.physics.plasma_geometry import PlasmaGeom
+from process.models.physics import impurity_radiation
 
 logger = logging.getLogger(__name__)
 
@@ -166,8 +167,16 @@ def calculate_confinement_time(
         try:
             model = ConfinementRadiationLossModel(int(self.data.physics.i_rad_loss))
 
+            # pden_plasma_core_rad_mw was reduced here, but if full radiation, not used!
+            # rad_reduction_for_tauE_only option allows full radiation model in PPB,
+            # but reduced radiation here in confinement time calculation
             if model == ConfinementRadiationLossModel.FULL_RADIATION:
-                p_plasma_loss_mw -= self.data.physics.pden_plasma_rad_mw * vol_plasma
+                if impurity_radiation.rad_reduction_for_tauE_only:
+                    # Reduced rad
+                    p_plasma_loss_mw -= pden_plasma_core_rad_mw * vol_plasma
+                else:
+                    # Not reduced rad
+                    p_plasma_loss_mw -= self.data.physics.pden_plasma_rad_mw * vol_plasma
             elif model == ConfinementRadiationLossModel.CORE_ONLY:
                 p_plasma_loss_mw -= pden_plasma_core_rad_mw * vol_plasma
             # NO_RADIATION: do not adjust p_plasma_loss_mw for radiation

process/models/physics/impurity_radiation.py

@@ -22,6 +22,10 @@
     from process.models.physics.plasma_profiles import PlasmaProfile
 
 logger = logging.getLogger(__name__)
+include_edge_radiation = False
+int_edge_rad = False
+rho_fix = False
+rad_reduction_for_tauE_only = False
 
 
 def initialise_imprad(data: DataStructure):
@@ -629,6 +633,12 @@ def element2index(element: str, data: DataStructure):
         ) from e
 
 
+# Globals for ease of extraction for investigation plotting only
+global pden_impurity_rad_profile
+global pden_impurity_core_rad_profile
+global pden_impurity_rad_edge_profile
+
+
 class ImpurityRadiation:
     """Calculates the impurity radiation losses for given temperature and
     density profiles. The considers the  total impurity radiation from the core
@@ -660,13 +670,17 @@ def __init__(self, plasma_profile: PlasmaProfile, data_structure: DataStructure)
         self.pden_impurity_core_rad_profile = np.zeros(
             self.data.physics.n_plasma_profile_elements
         )
+        self.pden_impurity_core_rad_profile_tauE = np.zeros(
+            self.data.physics.n_plasma_profile_elements
+        )
         self.pden_impurity_rad_edge_profile = np.zeros(
             self.data.physics.n_plasma_profile_elements
         )
 
         self.pden_impurity_rad_total_mw = 0.0
         self.pden_impurity_core_rad_total_mw = 0.0
         self.pden_impurity_rad_edge_total_mw = 0.0
+        self.pden_impurity_core_rad_total_tauE_mw = 0.0
 
     def run(self):
         """ImpurityRadiation model isn't run"""
@@ -705,25 +719,94 @@ def calculate_radiation_loss_profiles(self):
         radiation  (pden_impurity_rad_total_mw). Update the stored arrays with the
         values.
         """
-        pden_impurity_rad_total = (
-            self.pden_impurity_radiation_profile
-            * self.plasma_profile.neprofile.profile_x
+        # Core region radiation profile
+        # Multiplication by "profile_x" (formerly rho, normalised minor radius)
+        # wrong here: causes 0 power density at rho = 0. Should be performed in
+        # power integral instead
+        if rho_fix:
+            rho = np.ones_like(self.plasma_profile.neprofile.profile_x)
+        else:
+            rho = self.plasma_profile.neprofile.profile_x
+
+        # Optionally treat core radiation reduction separately for tauE calculation and
+        # plasma power balance
+        f_p_plasma_core_rad_reduction_tauE = (
+            self.data.impurity_radiation.f_p_plasma_core_rad_reduction
         )
+        if rad_reduction_for_tauE_only:
+            # Only reduce radiation for tauE calculation
+            f_p_plasma_core_rad_reduction = 1.0
+        else:
+            # Reduce radiation in PPB as well
+            f_p_plasma_core_rad_reduction = f_p_plasma_core_rad_reduction_tauE
+
         pden_impurity_core_rad_total = self.pden_impurity_radiation_profile * (
-            self.plasma_profile.neprofile.profile_x
+            rho
             * create_f_rad_core_profile(
                 rho=self.plasma_profile.neprofile.profile_x,
                 radius_plasma_core_norm=self.data.impurity_radiation.radius_plasma_core_norm,
-                f_p_plasma_core_rad_reduction=self.data.impurity_radiation.f_p_plasma_core_rad_reduction,
+                f_p_plasma_core_rad_reduction=f_p_plasma_core_rad_reduction,
             )
         )
+        pden_impurity_core_rad_total_tauE = self.pden_impurity_radiation_profile * (
+            rho
+            * create_f_rad_core_profile(
+                rho=self.plasma_profile.neprofile.profile_x,
+                radius_plasma_core_norm=self.data.impurity_radiation.radius_plasma_core_norm,
+                f_p_plasma_core_rad_reduction=f_p_plasma_core_rad_reduction_tauE,
+            )
+        )
+
+        if include_edge_radiation:
+            # Explicitly include edge radiation
+            # Edge region radiation profile
+            fradedge_profile = np.zeros_like(self.plasma_profile.neprofile.profile_x)
+            edge_mask = (
+                self.plasma_profile.neprofile.profile_x
+                >= self.data.impurity_radiation.radius_plasma_core_norm
+            )
+            fradedge_profile[edge_mask] = 1.0  # Edge region gets full value
+            pden_impurity_rad_edge_total = (
+                self.pden_impurity_radiation_profile * fradedge_profile
+            )
+
+            # Total radiation profile (core + edge)
+            pden_impurity_rad_total = (
+                pden_impurity_core_rad_total + pden_impurity_rad_edge_total
+            )
+
+            self.pden_impurity_rad_edge_profile = np.add(
+                self.pden_impurity_rad_edge_profile, pden_impurity_rad_edge_total
+            )
+        else:
+            # Old case: don't explicitly calculate edge radiation density
+            if rho_fix:
+                pden_impurity_rad_total = self.pden_impurity_radiation_profile
+
+            else:
+                pden_impurity_rad_total = (
+                    self.pden_impurity_radiation_profile
+                    * self.plasma_profile.neprofile.profile_x
+                )
 
         self.pden_impurity_rad_profile = np.add(
             self.pden_impurity_rad_profile, pden_impurity_rad_total
         )
         self.pden_impurity_core_rad_profile = np.add(
             self.pden_impurity_core_rad_profile, pden_impurity_core_rad_total
         )
+        self.pden_impurity_core_rad_profile_tauE = np.add(
+            self.pden_impurity_core_rad_profile_tauE, pden_impurity_core_rad_total_tauE
+        )
+        global pden_impurity_rad_profile
+        global pden_impurity_core_rad_profile
+        global pden_impurity_rad_edge_profile
+        global pden_impurity_core_rad_profile_tauE
+
+        pden_impurity_rad_profile = self.pden_impurity_rad_profile
+        pden_impurity_core_rad_profile = self.pden_impurity_core_rad_profile
+        pden_impurity_rad_edge_profile = self.pden_impurity_rad_edge_profile
+        pden_impurity_core_rad_profile_tauE = self.pden_impurity_core_rad_profile_tauE
 
     def integrate_radiation_loss_profiles(self):
         """Integrate the radiation loss profiles using the Simpson rule.
@@ -734,16 +817,37 @@ def integrate_radiation_loss_profiles(self):
         # but are correct:
         # see github.com/ukaea/PROCESS/issues/3968#issuecomment-3491154712
         # and github.com/ukaea/PROCESS/issues/3968#issuecomment-4935567006
+
+        # Old case: rho multiplication already performed incorrectly in power
+        # density calculation: don't multiply again here
+        # New case (rho_fix): multiply correct power density here by rho
+        # for integration
+        rho = 1.0
+        if rho_fix:
+            rho = self.plasma_profile.neprofile.profile_x
         self.pden_impurity_rad_total_mw = 2.0e-6 * integrate.simpson(
-            self.pden_impurity_rad_profile,
+            self.pden_impurity_rad_profile * rho,
             x=self.plasma_profile.neprofile.profile_x,
             dx=self.plasma_profile.neprofile.profile_dx,
         )
         self.pden_impurity_core_rad_total_mw = 2.0e-6 * integrate.simpson(
-            self.pden_impurity_core_rad_profile,
+            self.pden_impurity_core_rad_profile * rho,
             x=self.plasma_profile.neprofile.profile_x,
             dx=self.plasma_profile.neprofile.profile_dx,
         )
+        self.pden_impurity_core_rad_total_tauE_mw = 2.0e-6 * integrate.simpson(
+            self.pden_impurity_core_rad_profile_tauE * rho,
+            x=self.plasma_profile.neprofile.profile_x,
+            dx=self.plasma_profile.neprofile.profile_dx,
+        )
+
+        if include_edge_radiation:
+            # Integrate edge explicitly
+            self.pden_impurity_rad_edge_total_mw = 2.0e-6 * integrate.simpson(
+                self.pden_impurity_rad_edge_profile * rho,
+                x=self.plasma_profile.neprofile.profile_x,
+                dx=self.plasma_profile.neprofile.profile_dx,
+            )
 
     def calculate_imprad(self):
         """Call the map function to calculate impurity radiation parameters for each

process/models/physics/physics.py

@@ -751,6 +751,9 @@ def run(self):
         self.data.physics.pden_plasma_core_rad_mw = radpwrdata.pden_plasma_core_rad_mw
         self.data.physics.pden_plasma_outer_rad_mw = radpwrdata.pden_plasma_outer_rad_mw
         self.data.physics.pden_plasma_rad_mw = radpwrdata.pden_plasma_rad_mw
+        self.data.physics.pden_plasma_core_rad_tauE_mw = (
+            radpwrdata.pden_plasma_core_rad_tauE_mw
+        )
 
         self.data.physics.p_plasma_sync_mw = (
             self.data.physics.pden_plasma_sync_mw * self.data.physics.vol_plasma
@@ -871,6 +874,10 @@ def run(self):
 
         # Calculate transport losses and energy confinement time using the
         # chosen scaling law
+        # Reduce pden to effectively modify the highly radiative regime confinement
+        # time scaling
+        reduced_pden = self.data.physics.pden_plasma_core_rad_tauE_mw
+
         confinement_time_data = self.confinement.calculate_confinement_time(
             m_fuel_amu=self.data.physics.m_fuel_amu,
             p_alpha_total_mw=self.data.physics.p_alpha_total_mw,
@@ -887,7 +894,7 @@ def run(self):
             p_non_alpha_charged_mw=self.data.physics.p_non_alpha_charged_mw,
             p_hcd_injected_total_mw=self.data.current_drive.p_hcd_injected_total_mw,
             plasma_current=self.data.physics.plasma_current,
-            pden_plasma_core_rad_mw=self.data.physics.pden_plasma_core_rad_mw,
+            pden_plasma_core_rad_mw=reduced_pden,
             rmajor=self.data.physics.rmajor,
             rminor=self.data.physics.rminor,
             temp_plasma_electron_density_weighted_kev=self.data.physics.temp_plasma_electron_density_weighted_kev,