diff --git a/process/core/caller.py b/process/core/caller.py index 1634777e53..353d592103 100644 --- a/process/core/caller.py +++ b/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 diff --git a/process/data_structure/numerics.py b/process/data_structure/numerics.py index 57505f1c39..23d13c4841 100644 --- a/process/data_structure/numerics.py +++ b/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""" diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index 22f742e3ad..075a1c1842 100644 --- a/process/data_structure/physics_variables.py +++ b/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)""" diff --git a/process/main.py b/process/main.py index 2501505685..a8df683be6 100644 --- a/process/main.py +++ b/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 " diff --git a/process/models/costs/costs.py b/process/models/costs/costs.py index 208f98cbe0..afc89b580f 100644 --- a/process/models/costs/costs.py +++ b/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 diff --git a/process/models/physics/confinement_time.py b/process/models/physics/confinement_time.py index 29829b65d4..253a06eb02 100644 --- a/process/models/physics/confinement_time.py +++ b/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 diff --git a/process/models/physics/impurity_radiation.py b/process/models/physics/impurity_radiation.py index dd8d1fcbac..3a4c07811b 100644 --- a/process/models/physics/impurity_radiation.py +++ b/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,6 +670,9 @@ 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 ) @@ -667,6 +680,7 @@ def __init__(self, plasma_profile: PlasmaProfile, data_structure: DataStructure) 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,18 +719,75 @@ 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 @@ -724,6 +795,18 @@ def calculate_radiation_loss_profiles(self): 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 diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 4b3bef40be..7e125ca687 100644 --- a/process/models/physics/physics.py +++ b/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, diff --git a/process/models/physics/radiation_power.py b/process/models/physics/radiation_power.py index 0b19e1451b..bc38022188 100644 --- a/process/models/physics/radiation_power.py +++ b/process/models/physics/radiation_power.py @@ -14,6 +14,7 @@ from process.models.physics.plasma_profiles import PlasmaProfile logger = logging.getLogger(__name__) +include_edge_radiation = False @dataclass @@ -24,6 +25,7 @@ class RadpwrData: pden_plasma_core_rad_mw: float pden_plasma_outer_rad_mw: float pden_plasma_rad_mw: float + pden_plasma_core_rad_tauE_mw: float def calculate_radiation_powers( @@ -103,9 +105,14 @@ def calculate_radiation_powers( imp_rad = impurity.ImpurityRadiation(plasma_profile, data_structure) imp_rad.calculate_imprad() - pden_plasma_outer_rad_mw = ( - imp_rad.pden_impurity_rad_total_mw - imp_rad.pden_impurity_core_rad_total_mw - ) + if impurity.include_edge_radiation: + # Use calculated edge radiation + pden_plasma_outer_rad_mw = imp_rad.pden_impurity_rad_edge_total_mw + else: + # Use difference between total and core radiation + pden_plasma_outer_rad_mw = ( + imp_rad.pden_impurity_rad_total_mw - imp_rad.pden_impurity_core_rad_total_mw + ) # Synchrotron radiation power/volume; assumed to be from core only. pden_plasma_sync_mw = psync_albajar_fidone( @@ -127,6 +134,9 @@ def calculate_radiation_powers( pden_plasma_core_rad_mw = ( imp_rad.pden_impurity_core_rad_total_mw + pden_plasma_sync_mw ) + pden_plasma_core_rad_tauE_mw = ( + imp_rad.pden_impurity_core_rad_total_tauE_mw + pden_plasma_sync_mw + ) # Total radiation power/volume. pden_plasma_rad_mw = imp_rad.pden_impurity_rad_total_mw + pden_plasma_sync_mw @@ -136,6 +146,7 @@ def calculate_radiation_powers( pden_plasma_core_rad_mw, pden_plasma_outer_rad_mw, pden_plasma_rad_mw, + pden_plasma_core_rad_tauE_mw, )