diff --git a/documentation/source/development/add-vars.md b/documentation/source/development/add-vars.md index f6cbd1a625..e925ea093d 100644 --- a/documentation/source/development/add-vars.md +++ b/documentation/source/development/add-vars.md @@ -1,147 +1,170 @@ -# Guide for adding Variables & Constraints +# Guide for adding variables and constraints -Specific instructions must be followed to add an input, iteration variable, -optimisation figure of merit and constraints to the `PROCESS` code. +A guide for how to add new inputs, constraints, figures of merit, and scan variables. - **At all times the [`PROCESS` style guide](../development/standards.md) must be used.** +**At all times ensure new variables and functions adhere to the [`PROCESS` style guide](../development/standards.md).** -!!! note - As the code is quickly converging towards a wholly Python codebase the respective files may change in type from `.f90` to `.py`. +All of these features rely on 'variables' which belong to a 'data structure'. All of the data structures can be found in `process/data_structure`. ------------------ +In general, a variable within a data structure could act as an: -## Add an input +- Input variable: is specified by the user in an `IN.DAT` and its value is not changed once PROCESS is running. +- Iteration variable: is modified by the solver to try and optimise for some figure of merit. +- Scan variable: is sequentially modified by the `Scan` class to some `IN.DAT`-defined values. +- Intermediate variable: is calculated within a model and then used within other models. +- Output variable: is calculated within a model and then written out to the `MFILE.DAT` and `OUT.DAT`. -To add a `PROCESS` input, please follow below: +It is advised that a variable is either used to define a particular PROCESS run (input variable, iteration variable, or scan variable) or mutated within a PROCESS run (output variable or intermediate variable). Mixing the two classes of variable (e.g. having a variable that can be input but is also mutated within a model) will lead to confusing, dangerous, and incorrect results. -1. Choose the most relevant module `XX` and add the variable in the `XX_variables` defined in `XX_variables.f90`. - -2. Add a description of the input variable below the declaration, using the FORD formatting described in the standards section specifying the units. - -3. Specify a sensible default value in the `init_XX_variables()` function within the corresponding model `.py` main file - -4. Add the parameter to the `INPUT_VARIABLES` dictionary in `input.py`. -Here is an example of the code to add: - +----------------- + +## Add a new variable +You may need to add a variable to PROCESS when changing or creating models. In most cases, you will want to add your variable to an existing data structure. Creating an entirely new data structure is beyond the scope of this guide, so please seek support from the PROCESS maintainers. + +For example, if you are adding a new variable that relates to the blanket model, you would add the variable to `process/data_structure/blanket_variables.py` as part of the `BlanketData` dataclass. + +```python +@dataclass(slots=True) +class BlanketData: + <... existing variables ...> -Variable definition example in `tfcoil_variables.f90`: -```fortran - real(dp) :: rho_tf_joints - !! TF joints surfacic resistivity [ohm.m] - !! Feldmetal joints assumed. + my_new_blanket_variable: float = 0.0 + """my variable description [m]""" ``` -Variable initialization example in `tf_coil.py`: +Here, `[m]` is the units and should be replaced with the appropriate units for the variable being added. + +This variable could then be used within a model + ```python - def init_tfcoil_variables(): - ... - tfv.rho_tf_joints = 2.5e-10 +self.data.blanket.my_new_blanket_variable = 1.0 +... +another_variable = self.data.blanket.my_new_blanket_variable / 2.0 ``` -Code example in the `input.py` file: +----------------- + +## Add a new input +Adding an input in PROCESS means that some variable in a data structure can be set from the `IN.DAT`. Inputs are defined in the `process/core/input.py` file in the `INPUT_VARIABLES` dictionary. Adding a new entry to this dictionary will create a new input. +Continuing with the example from the previous section: ```python - INPUT_VARIABLES = { +INPUT_VARIABLES = { ... - "rho_tf_joints": InputVariable("tfcoil", float, range=(0.0, 0.01)), + "my_new_blanket_variable": InputVariable("blanket", float), +} ``` ------------------ +`InputVariable` has several additional fields to support validation and the parsing of arrays; please consult the dataclass for these additional arguments. -## Add an iteration variable +You would replace `"blanket"` with the name of the data structure your specific variable belongs to (found by looking at `DataStructure` in `process/core/model.py`). + +Now, in the `IN.DAT`, you could set an initial value for `my_new_blanket_variable` by writing: -To add a `PROCESS` iteration variable please follow the steps below, in addition to the instructions for adding an input variable: +``` +my_new_blanket_variable = 1.0 +``` +Equally, you could choose not to specify the variable in the `IN.DAT` and the value would be initialised to its default value (`0.0` in this example, as we specified [earlier](#add-a-new-variable)). -1. The parameter `IPNVARS` in module `numerics` of `numerics.f90` will normally be greater than the actual number of iteration variables, and does not need to be changed. -2. Append a new iteration number key to the end of the `ITERATION_VARIABLES` dictionary in `iteration_variables.py`. The associated variable is the corresponding key value. -3. Set the variable origin file and then the associated lower and upper bounds -4. Update the `lablxc` description in `numerics.f90`. - -It should be noted that iteration variables must not be reset elsewhere in the -code. That is, they may only be assigned new values when originally -initialised (in the relevant module, or in the input file if required). -Otherwise, the numerical procedure cannot adjust the value as it requires, and -the program will fail. +----------------- + +## Add an iteration variable + +Adding an iteration variable allows the PROCESS solver to change the variable as part of the optimisation/solving loop. Iteration variables are defined in `process/core/solver/iteration_variables.py` in the `ITERATION_VARIABLES` dictionary. You would add a new entry to this dictionary to create a new iteration variable: -Here is a code snippet showing how `rmajor` is defined in `iteration_variables.py` ```python ITERATION_VARIABLES = { - ... - 3: IterationVariable("rmajor", "physics", 0.1, 50.00), + 123: IterationVariable("my_new_blanket_variable", "blanket", 0.1, 1.0), +} ``` +In this example: + +- `123` is the identifier of the iteration variable, and must be unique. +- `"blanket"` is the data structure the variable will be set on. +- `0.1` is the default lower bound of the variable. +- `1.0` is the default upper bound of the variable. + +You will often want to [add a variable to the input file](#add-a-new-input) if it is an iteration variable. That way, you can specify the initial value of the iteration variable in the `IN.DAT`. + +!!! note + Iteration variables are the exception to the best-practice of not allowing inputs to be mutated. The solver will obviously need to modify iteration variables away from their default value, however it remains true that _models_ should not mutate inputs. Furthermore, models should not change the values of iteration variables set by the solver. + +The iteration variable can be enabled in the `IN.DAT` by: +``` +ixc = 123 + +my_new_blanket_variable = 0.5 +``` + +Again, note that omitting the line specifying the initial value for an iteration variable means the initial value would be set to that variable's default value that we specified [earlier](#add-a-new-variable). + ----------------- ## Add a figure of merit -New figures of merit are added to `PROCESS` in the following way: +A figure of merit is the scalar that the optimiser (e.g. VMCON) will try and minimise or maximise. The figures of merit are specified in `process/core/solver/objectives.py` in the `objective_function()` function. -1. Increment the parameter `IPNFOMS` in module `numerics` in source file `numerics.py` to accommodate the new figure of merit. - -2. Assign the new integer value and description string of the new figure of merit to the `FiguresOfMerit` enumerator in `numerics.py`. - -3. Add the new figure of merit equation to `objective_function()` in `objectives.py`, following the method used in the existing examples. The value of figure of merit case should be of order unity, so select a reasonable scaling factor if necessary. - -An example can be found below: +To add a new figure of merit, first create a new entry in the `FiguresOfMerit` enum in `process/data_structure/numerics.py`: +```python +class FiguresOfMerit(IntEnum): + ... + BLANKET_FIGURE_OF_MERIT = (20, "my FOM description") +``` +Here `20` will be the identifier of the figure of merit, and **must** be unique. +Finally, add the equation to `process/core/solver/objectives.py`: ```python -objective_function(): - ... - try: - figure_of_merit = FiguresOfMerit(abs(minmax)) - ... - if figure_of_merit == FiguresOfMerit.MAJOR_RADIUS: - objective_metric = 0.2 * data.physics.rmajor +elif figure_of_merit == FiguresOfMerit.BLANKET_FIGURE_OF_MERIT: + objective_metric = data.blanket.my_new_blanket_variable ``` ------------ +Note that you will want to scale the `objective_metric` such that it is of the order unity if the variable is not already. This is recommended (but not necessary) because most optimisers are tuned to work well on data of this scale. -## Add a scan variable +The figure of merit can be selected in the `IN.DAT`: +``` +minmax = 20 +``` +Remember, setting `minmax = -20` would minimise instead of maximise our new variable. -After following the instruction to add an input variable, you can make the variable a scan variable by following these steps: +----------------- -1. Increment the parameter `IPNSCNV` defined in `scan_variables.py` in the data_structure directory, to accommodate the new scanning variable. The incremented value will identify your scan variable. - -2. Add a short description of the new scanning variable in the `nsweep` comment in `scan_variables.py`, alongside its identification number. - -3. Update the `ScanVariables` enum in the `scan.py` file by adding a new case statement connecting the variable to the scan integer switch, the variable name and a short description. - -4. Add a comment in the corresponding variable file in the data_structure directory, eg, `data_structure/[XX]_variables.py`, to add the variable description indicating the scan switch number. - +## Add a scan variable -`nsweep` comment example: -```fortran +After following the instructions to add an input variable, you can then make a scan variable. - integer :: nsweep = 1 - !! nsweep /1/ : switch denoting quantity to scan: +First, add the variable to the `ScanVariables` enum in `process/core/scan.py`. +```python +class ScanVariables(Enum): + ... + blanket_scan_variable = ScanVariable( + "my_new_blanket_variable", "A blanket variable", 82 + ) ``` +Here, `82` is the identifier of the scan variable and must be unique. -`SCAN_VARIABLES` case example: +Next, increment the parameter `IPNSCNV` in `process/data_structure/scan_variables.py` and be sure to add a description of the scan variable in the docstring of the `nsweep` variable. +Finally, in `process/core/scan.py`, add the scan variable to the `Scan.scan_select()` method. ```python - class ScanVariables(Enum): - aspect: ScanVariable("aspect", "Aspect_ratio", 1), - pflux_div_heat_load_max_mw: ScanVariable("pflux_div_heat_load_max_mw", "Div_heat_limit_(MW/m2)", 2), - ... - Bc2_0K: ScanVariable("Bc2(0K)", "GL_NbTi Bc2(0K)", 54), - dr_shld_inboard : ScanVariable("dr_shld_inboard", "Inboard neutronic shield", 55), +match nwp: + ... + case 82: + self.data.tfcoil.my_new_blanket_variable = swp[iscn - 1] ``` ---------------- +Please see the [scan documentation](../io/input-guide.md#input-guide) for how to set up a scan `IN.DAT` + +----------------- ## Add a constraint equation -Constraint equations are added to *PROCESS* in the `process/core/solver/constraints.py` file. They are registered with the `ConstraintManager` whenever the application is run. Each equation has a unique name that is currently an integer, however upgrades to the input file format in the future will allow arbitrary hashable constraint names. +Constraint equations are added to PROCESS in the `process/core/solver/constraints.py` file. They are registered with the `ConstraintManager` whenever the application is run. Each equation has a unique name that is currently an integer. A constraint is simply added by registering the constraint to the manager using a decorator. @@ -149,24 +172,25 @@ A constraint is simply added by registering the constraint to the manager using @ConstraintManager.register_constraint(1234, "m", "=") def my_constraint_function(constraint_registration): ... ``` -The arguments to the `register_constraint` function are: -- Name (again, currently an integer) -- Unit (for output reporting purposes) -- Symbol (e.g. =, >=, <=. Again, for output reporting purposes) +The arguments to the `register_constraint()` function are: +- Name (again, currently an integer): in the example above it is `1234` +- Unit (for output reporting purposes): in the example above it is `m` +- Symbol (e.g. `=`, `>=`, `<=`. Again, for output reporting purposes): in the example above it is `=` -`my_constraint_function` should be named appropriately and return a `ConstraintResult` which contains the: -- Normalised residual error +`my_constraint_function()` should be named appropriately and return a `ConstraintResult` which contains the: + +- Constraint residual +- Normalised residual - Constraint value - Constraint bound -- Constraint residual -The recommended way to do this is using one of the functions `geq`, `leq`, or `eq` depending on whether the constraint is desired to be $v\geq b$, $v\leq b$, or $v=b$, respectively. +The recommended way to do this is using one of the functions `geq`, `leq`, or `eq` depending on whether the constraint is desired to be $v\geq b$, $v\leq b$, or $v=b$, respectively, where $v$ is the value (something that PROCESS calculates) and $b$ is the bound (the limit set by the user in the `IN.DAT`). ```python @ConstraintManager.register_constraint(1234, "m", "=") def my_constraint_function(constraint_registration): - return geq(value, bound, constraint_registration) + return eq(value, bound, constraint_registration) ``` diff --git a/process/core/input.py b/process/core/input.py index 09138f9072..ca17ab6566 100644 --- a/process/core/input.py +++ b/process/core/input.py @@ -15,7 +15,7 @@ ) from process.core.solver.constraints import ConstraintManager from process.data_structure.impurity_radiation_variables import N_IMPURITIES -from process.data_structure.numerics import IPEQNS, IPNVARS +from process.data_structure.numerics import IPNVARS from process.data_structure.pfcoil_variables import N_PF_GROUPS_MAX from process.data_structure.physics_variables import N_CONFINEMENT_SCALINGS from process.data_structure.scan_variables import IPNSCNS, IPNSCNV @@ -48,7 +48,7 @@ def _icc_additional_actions( data.numerics.n_constraints += 1 -@dataclass +@dataclass(slots=True) class InputVariable: """A variable to be parsed from the input file.""" @@ -1165,7 +1165,7 @@ def bounds(self) -> tuple[NumberType | None, NumberType | None]: "icc": InputVariable( None, int, - range=(1, IPEQNS), + choices=ConstraintManager.constraint_ids(), additional_actions=_icc_additional_actions, set_variable=False, ), diff --git a/process/core/solver/constraints.py b/process/core/solver/constraints.py index a40aa75767..481f2aca16 100644 --- a/process/core/solver/constraints.py +++ b/process/core/solver/constraints.py @@ -84,6 +84,10 @@ def num_constraints(cls): """Return the number of constraints currently in the registry""" return len(cls._constraint_registry) + @classmethod + def constraint_ids(cls): + return tuple(cls._constraint_registry.keys()) + @classmethod def register_constraint( cls, name: Hashable, units: str, symbol: ConstraintSymbolType diff --git a/process/core/solver/iteration_variables.py b/process/core/solver/iteration_variables.py index 490022a270..6324b43b62 100644 --- a/process/core/solver/iteration_variables.py +++ b/process/core/solver/iteration_variables.py @@ -1,12 +1,16 @@ +from __future__ import annotations + import logging from copy import deepcopy from dataclasses import dataclass -from typing import Any +from typing import TYPE_CHECKING, Any import numpy as np from process.core.exceptions import ProcessValueError -from process.core.model import DataStructure + +if TYPE_CHECKING: + from process.core.model import DataStructure logger = logging.getLogger(__name__) diff --git a/process/data_structure/numerics.py b/process/data_structure/numerics.py index 8daf03fb78..2bdec08e41 100644 --- a/process/data_structure/numerics.py +++ b/process/data_structure/numerics.py @@ -4,6 +4,8 @@ import numpy as np +from process.core.solver.iteration_variables import ITERATION_VARIABLES + class PROCESSRunMode(IntEnum): """Enumeration of the available PROCESS run modes, which determine the behaviour @@ -100,14 +102,12 @@ def description(self): return self._description_ -IPNVARS = 177 +IPNVARS = max(ITERATION_VARIABLES.keys()) """total number of variables available for iteration""" -IPEQNS = 92 -"""number of constraint equations available""" - -IPNFOMS = 19 -"""number of available figures of merit""" +# Set to a really large number so that it should never need to be changed +IPEQNS = 500 +"""Maximum number of constraint equations available""" @dataclass(slots=True) diff --git a/process/models/blankets/hcpb.py b/process/models/blankets/hcpb.py index 6cd5d2e862..8be6c345fb 100644 --- a/process/models/blankets/hcpb.py +++ b/process/models/blankets/hcpb.py @@ -586,12 +586,6 @@ def nuclear_heating_magnets(self, output: bool): "(p_tf_nuclear_heat_mw.)", self.data.fwbs.p_tf_nuclear_heat_mw, ) - po.ovarre( - self.outfile, - "p_fusion_total_mw", - "(p_fusion_total_mw.)", - self.data.physics.p_fusion_total_mw, - ) po.ovarre( self.outfile, "total mass of the TF coils (kg)", diff --git a/process/models/engineering/materials.py b/process/models/engineering/materials.py index d42f673b93..fef43d7e1d 100644 --- a/process/models/engineering/materials.py +++ b/process/models/engineering/materials.py @@ -2,6 +2,7 @@ import logging +import numba import numpy as np logger = logging.getLogger(__name__) @@ -48,6 +49,7 @@ def eurofer97_thermal_conductivity(temp: float, fw_th_conductivity: float) -> fl ) +@numba.njit(cache=True) def calculate_tresca_stress( stress_x: float | np.ndarray, stress_y: float | np.ndarray, @@ -71,13 +73,16 @@ def calculate_tresca_stress( Tresca stress (maximum shear stress criterion) in Pa, defined as the maximum of |stress_x - stress_y|, |stress_y - stress_z|, |stress_x - stress_z|. """ - return max( - abs(stress_x - stress_y), - abs(stress_y - stress_z), - abs(stress_x - stress_z), + return np.maximum( + np.maximum( + np.abs(stress_x - stress_y), + np.abs(stress_y - stress_z), + ), + np.abs(stress_x - stress_z), ) +@numba.njit(cache=True) def calculate_von_mises_stress( stress_x: float | np.ndarray, stress_y: float | np.ndarray, diff --git a/process/models/physics/impurity_radiation.py b/process/models/physics/impurity_radiation.py index f4b5363134..567a64c244 100644 --- a/process/models/physics/impurity_radiation.py +++ b/process/models/physics/impurity_radiation.py @@ -655,8 +655,10 @@ def integrate_radiation_loss_profiles(self): """Integrate the radiation loss profiles using the Simpson rule. Store the total values for each aspect of impurity radiation loss. """ - # 2.0e-6 converts from W/m^3 to MW/m^3 and also accounts for both sides of the - # plasma + # 1e-6 converts from W/m^3 to MW/m^3 + # The factor 2 below and and normalised radius profile_x above may be unexpected, but are correct: + # see github.com/ukaea/PROCESS/issues/3968#issuecomment-3491154712 + # and github.com/ukaea/PROCESS/issues/3968#issuecomment-4935567006 self.pden_impurity_rad_total_mw = 2.0e-6 * integrate.simpson( self.pden_impurity_rad_profile, x=self.plasma_profile.neprofile.profile_x, diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 49995a5ed7..6d0e0cfafb 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -2289,20 +2289,7 @@ def outplas(self): self.data.physics.f_alpha_ion, "OP ", ) - po.ovarre( - self.outfile, - "Ion transport (MW)", - "(p_ion_transport_loss_mw)", - self.data.physics.p_ion_transport_loss_mw, - "OP ", - ) - po.ovarre( - self.outfile, - "Electron transport (MW)", - "(p_electron_transport_loss_mw)", - self.data.physics.p_electron_transport_loss_mw, - "OP ", - ) + # Ion and electron transpor are now output belowin power accounting po.ovarre( self.outfile, "Injection power to ions (MW)", @@ -2323,6 +2310,422 @@ def outplas(self): ): po.ocmmnt(self.outfile, " (Injected power only used for start-up phase)") + # Global power imbalance output #4233 + po.oheadr(self.outfile, "Power accounting") + po.ocmmnt( + self.outfile, + "See Figure in https://ukaea.github.io/PROCESS/physics-models/plasma_power_balance", + ) + p_loss_mw = ( + self.data.current_drive.f_p_beam_orbit_loss + + self.data.current_drive.p_beam_shine_through_mw + + self.data.physics.p_fw_alpha_mw + ) + p_plasma_out = ( + self.data.physics.p_electron_transport_loss_mw + + self.data.physics.p_ion_transport_loss_mw + + p_loss_mw + + self.data.physics.p_plasma_rad_mw + ) + + p_plasma_in = ( + self.data.physics.p_alpha_total_mw + + self.data.physics.p_non_alpha_charged_mw + + self.data.current_drive.p_hcd_injected_total_mw + + self.data.physics.p_plasma_ohmic_mw + ) + p_plasma_imbalance_mw = p_plasma_in - p_plasma_out + po.oshead(self.outfile, "Plasma power balance across separatrix") + po.ocmmnt(self.outfile, "IN") + po.ovarre( + self.outfile, + "Alpha power (MW)", + "(p_alpha_total_mw)", + self.data.physics.p_alpha_total_mw, + "OP ", + ) + po.ovarre( + self.outfile, + "Power from p, 3He, T products of DD and/or D-He3 fusion (MW)", + "(p_non_alpha_charged_mw)", + self.data.physics.p_non_alpha_charged_mw, + "OP ", + ) + po.ovarre( + self.outfile, + "Injected power (MW)", + "(p_hcd_injected_total_mw)", + self.data.current_drive.p_hcd_injected_total_mw, + "OP ", + ) + po.ovarre( + self.outfile, + "Ohmic heating (MW)", + "(p_plasma_ohmic_mw)", + self.data.physics.p_plasma_ohmic_mw, + "OP ", + ) + po.ovarre( + self.outfile, + "TOTAL (MW)", + "(p_plasma_in)", + p_plasma_in, + "OP ", + ) + po.oblnkl(self.outfile) + po.ocmmnt(self.outfile, "OUT") + po.ovarre( + self.outfile, + "Net power transported by electrons (MW)", + "(p_electron_transport_loss_mw)", + self.data.physics.p_electron_transport_loss_mw, + "OP ", + ) + po.ovarre( + self.outfile, + "Net power transported by ions (MW)", + "(p_ion_transport_loss_mw)", + self.data.physics.p_ion_transport_loss_mw, + "OP ", + ) + po.ovarre( + self.outfile, + "Power lost by beam ions and unthermalised alphas (MW)", + "(p_loss_mw)", + p_loss_mw, + "OP ", + ) + po.ovarre( + self.outfile, + "Total radiation loss, including net loss by synchrotron radiation (MW)", + "(p_plasma_rad_mw)", + self.data.physics.p_plasma_rad_mw, + "OP ", + ) + po.ovarre( + self.outfile, + "TOTAL (MW)", + "(p_plasma_out)", + p_plasma_out, + "OP ", + ) + po.ovarre( + self.outfile, + "Plasma power imbalance (MW)", + "(p_plasma_imbalance_mw)", + p_plasma_imbalance_mw, + "OP ", + ) + if abs(p_plasma_imbalance_mw) > 0.1: + logger.error("Plasma power imbalance > 0.1 MW") + + po.oshead(self.outfile, "Power balance for reactor") + p_reactor_in = ( + self.data.physics.p_fusion_total_mw + + self.data.fwbs.p_blkt_multiplication_mw + + self.data.current_drive.p_hcd_injected_total_mw + + self.data.physics.p_plasma_ohmic_mw + + self.data.primary_pumping.p_fw_blkt_coolant_pump_mw + + self.data.heat_transport.p_div_coolant_pump_mw + ) + p_reactor_out = ( + self.data.heat_transport.p_plant_primary_heat_mw + + self.data.heat_transport.p_div_secondary_heat_mw + + self.data.heat_transport.p_shld_secondary_heat_mw + + self.data.fwbs.p_tf_nuclear_heat_mw + + self.data.fwbs.p_fw_hcd_nuclear_heat_mw + + self.data.fwbs.p_fw_hcd_rad_total_mw + ) + p_reactor_imbalance_mw = p_reactor_in - p_reactor_out + po.ocmmnt(self.outfile, "IN") + po.ovarre( + self.outfile, + "Fusion power (MW)", + "(p_fusion_total_mw)", + self.data.physics.p_fusion_total_mw, + "OP ", + ) + po.ovarre( + self.outfile, + "Energy multiplication in blanket and shield (MW)", + "(p_blkt_multiplication_mw)", + self.data.fwbs.p_blkt_multiplication_mw, + "OP ", + ) + po.ovarre( + self.outfile, + "Injected power (MW)", + "(p_hcd_injected_total_mw)", + self.data.current_drive.p_hcd_injected_total_mw, + "OP ", + ) + po.ovarre( + self.outfile, + "Ohmic heating (MW)", + "(p_plasma_ohmic_mw)", + self.data.physics.p_plasma_ohmic_mw, + "OP ", + ) + po.ovarre( + self.outfile, + "Power deposited by pump in coolant for FW and blanket circuit (MW)", + "(p_fw_blkt_coolant_pump_mw)", + self.data.primary_pumping.p_fw_blkt_coolant_pump_mw, + "OP ", + ) + po.ovarre( + self.outfile, + "Power deposited by pump in coolant for divertor circuit (MW)", + "(p_div_coolant_pump_mw)", + self.data.heat_transport.p_div_coolant_pump_mw, + "OP ", + ) + po.ovarre( + self.outfile, + "TOTAL (MW)", + "(p_reactor_in)", + p_reactor_in, + "OP ", + ) + po.oblnkl(self.outfile) + po.ocmmnt(self.outfile, "OUT") + po.ovarre( + self.outfile, + "Total primary thermal power used for electricity production (MW)", + "(p_plant_primary_heat_mw)", + self.data.heat_transport.p_plant_primary_heat_mw, + "OP ", + ) + po.ovarre( + self.outfile, + "Divertor thermal power not used for electricity production (MW)", + "(p_div_secondary_heat_mw)", + self.data.heat_transport.p_div_secondary_heat_mw, + "OP ", + ) + po.ovarre( + self.outfile, + "Shield thermal power not used for electricity production (MW)", + "(p_shld_secondary_heat_mw)", + self.data.heat_transport.p_shld_secondary_heat_mw, + "OP ", + ) + po.ovarre( + self.outfile, + "Nuclear heating in TF coils (MW)", + "(p_tf_nuclear_heat_mw)", + self.data.fwbs.p_tf_nuclear_heat_mw, + "OP ", + ) + po.ovarre( + self.outfile, + "Nuclear heating in H&CD systems and diagnostics (MW)", + "(p_fw_hcd_nuclear_heat_mw)", + self.data.fwbs.p_fw_hcd_nuclear_heat_mw, + "OP ", + ) + po.ovarre( + self.outfile, + "Radiation heat deposited in H&CD systems and diagnostics (MW)", + "(p_fw_hcd_rad_total_mw)", + self.data.fwbs.p_fw_hcd_rad_total_mw, + "OP ", + ) + po.ovarre( + self.outfile, + "TOTAL (MW)", + "(p_reactor_out)", + p_reactor_out, + "OP ", + ) + po.ovarre( + self.outfile, + "Reactor power imbalance (MW)", + "(p_reactor_imbalance_mw)", + p_reactor_imbalance_mw, + "OP ", + ) + if abs(p_reactor_imbalance_mw) > 0.1: + logger.error("Reactor power imbalance > 0.1 MW") + + po.oshead(self.outfile, "Electrical power balance") + p_electric_demand = ( + self.data.heat_transport.p_plant_electric_net_mw + + self.data.heat_transport.p_tf_electric_supplies_mw + + self.data.pf_coil.p_pf_electric_supplies_mw + + self.data.heat_transport.p_hcd_electric_total_mw + + self.data.heat_transport.p_coolant_pump_elec_total_mw + + self.data.heat_transport.vachtmw + + self.data.heat_transport.p_cryo_plant_electric_mw + + self.data.heat_transport.p_tritium_plant_electric_mw + + self.data.heat_transport.fachtmw + ) + p_electric_imbalance = ( + self.data.heat_transport.p_plant_electric_gross_mw - p_electric_demand + ) + po.ocmmnt(self.outfile, "GENERATION") + po.ovarre( + self.outfile, + "Gross electric output (MW)", + "(p_plant_electric_gross_mw)", + self.data.heat_transport.p_plant_electric_gross_mw, + "OP ", + ) + po.oblnkl(self.outfile) + po.ocmmnt(self.outfile, "PURPOSE") + po.ovarre( + self.outfile, + "Net electric exported (MW)", + "(p_plant_electric_net_mw)", + self.data.heat_transport.p_plant_electric_net_mw, + "OP ", + ) + po.ovarre( + self.outfile, + "TF coils electric power (MW)", + "(p_tf_electric_supplies_mw)", + self.data.heat_transport.p_tf_electric_supplies_mw, + "OP ", + ) + po.ovarre( + self.outfile, + "PF coils electric power (MW)", + "(p_pf_electric_supplies_mw)", + self.data.pf_coil.p_pf_electric_supplies_mw, + "OP ", + ) + po.ovarre( + self.outfile, + "Heating and current drive electric power (MW)", + "(p_hcd_electric_total_mw)", + self.data.heat_transport.p_hcd_electric_total_mw, + "OP ", + ) + po.ovarre( + self.outfile, + "Primary coolant pump electric power (MW)", + "(p_coolant_pump_elec_total_mw)", + self.data.heat_transport.p_coolant_pump_elec_total_mw, + "OP ", + ) + po.ovarre( + self.outfile, + "Vacuum pumps electric power (MW)", + "(vachtmw)", + self.data.heat_transport.vachtmw, + "OP ", + ) + po.ovarre( + self.outfile, + "Cryoplant electric power (MW)", + "(p_cryo_plant_electric_mw)", + self.data.heat_transport.p_cryo_plant_electric_mw, + "OP ", + ) + po.ovarre( + self.outfile, + "Tritium plant electric power (MW)", + "(p_tritium_plant_elec_mw)", + self.data.heat_transport.p_tritium_plant_electric_mw, + "OP ", + ) + po.ovarre( + self.outfile, + "All other internal electric power requirements (MW)", + "(fachtmw)", + self.data.heat_transport.fachtmw, + "OP ", + ) + po.ovarre( + self.outfile, + "TOTAL (MW)", + "(p_electric_demand)", + p_electric_demand, + "OP ", + ) + po.oblnkl(self.outfile) + po.ovarre( + self.outfile, + "Electric power imbalance (MW)", + "(p_electric_imbalance)", + p_electric_imbalance, + "OP ", + ) + if abs(p_electric_imbalance) > 0.1: + logger.error("Electric power imbalance > 0.1 MW") + + po.oshead(self.outfile, "Power balance for power plant") + po.ocmmnt(self.outfile, "IN") + p_plant_in_mw = ( + self.data.physics.p_fusion_total_mw + self.data.fwbs.p_blkt_multiplication_mw + ) + p_plant_out_mw = ( + self.data.heat_transport.p_plant_electric_net_mw + + self.data.power.p_turbine_loss_mw + + self.data.heat_transport.p_plant_secondary_heat_mw + ) + p_plant_imbalance_mw = p_plant_in_mw - p_plant_out_mw + po.ovarre( + self.outfile, + "Fusion power (MW)", + "(p_fusion_total_mw)", + self.data.physics.p_fusion_total_mw, + "OP ", + ) + po.ovarre( + self.outfile, + "Energy multiplication in blanket and shield (MW)", + "(p_blkt_multiplication_mw)", + self.data.fwbs.p_blkt_multiplication_mw, + "OP ", + ) + po.ovarre( + self.outfile, + "TOTAL (MW)", + "(p_nuclear_total_mw)", + p_plant_in_mw, + "OP ", + ) + po.oblnkl(self.outfile) + po.ocmmnt(self.outfile, "OUT") + po.ovarre( + self.outfile, + "Net electric (MW)", + "(p_plant_electric_net_mw)", + self.data.heat_transport.p_plant_electric_net_mw, + "OP ", + ) + po.ovarre( + self.outfile, + "Heat rejected by main power conversion circuit (MW)", + "(p_turbine_loss_mw)", + self.data.power.p_turbine_loss_mw, + "OP ", + ) + po.ovarre( + self.outfile, + "Heat rejected by other circuits (secondary heat) (MW)", + "(p_plant_secondary_heat_mw)", + self.data.heat_transport.p_plant_secondary_heat_mw, + "OP ", + ) + po.ovarre( + self.outfile, + "TOTAL (MW)", + "(p_plant_out_mw)", + p_plant_out_mw, + "OP ", + ) + po.ovarre( + self.outfile, + "Power plant overall imbalance (MW)", + "(p_plant_imbalance_mw)", + p_plant_imbalance_mw, + "OP ", + ) + if abs(p_plant_imbalance_mw) > 0.1: + logger.error("Power plant overall imbalance > 0.1 MW") + self.exhaust.output() if self.data.stellarator.istell == 0: diff --git a/process/models/tfcoil/base.py b/process/models/tfcoil/base.py index 3fd5ed6deb..0163ee3695 100644 --- a/process/models/tfcoil/base.py +++ b/process/models/tfcoil/base.py @@ -19,6 +19,10 @@ from process.data_structure.build_variables import TFCSRadialConfiguration from process.data_structure.pfcoil_variables import PFConductorModel from process.data_structure.physics_variables import DivertorNumberModels +from process.models.engineering.materials import ( + calculate_tresca_stress, + calculate_von_mises_stress, +) logger = logging.getLogger(__name__) @@ -3135,19 +3139,20 @@ def stresscl( # Tresca / Von Mises yield criteria calculations # ----------------------------- # Array equation - s_shear_tf = np.maximum( - np.absolute(sig_tf_r - sig_tf_z), np.absolute(sig_tf_z - sig_tf_t) + + s_shear_tf = calculate_tresca_stress( + stress_x=sig_tf_r, stress_y=sig_tf_t, stress_z=sig_tf_z ) # Array equation - sig_tf_vmises = np.sqrt( - 0.5e0 - * ( - (sig_tf_r - sig_tf_t) ** 2 - + (sig_tf_r - sig_tf_z) ** 2 - + (sig_tf_z - sig_tf_t) ** 2 - ) + sig_tf_vmises = calculate_von_mises_stress( + stress_x=sig_tf_r, + stress_y=sig_tf_t, + stress_z=sig_tf_z, + stress_shear_xy=0.0, + stress_shear_yz=0.0, + stress_shear_zx=0.0, ) # Array equation @@ -3163,9 +3168,23 @@ def stresscl( ): # Addaped Von-mises stress calculation to WP strucure [Pa] - svmxz = sigvm(0.0e0, sig_tf_t[ii], sig_tf_z[ii], 0.0e0, 0.0e0, 0.0e0) + svmxz = calculate_von_mises_stress( + stress_x=0.0e0, + stress_y=sig_tf_t[ii], + stress_z=sig_tf_z[ii], + stress_shear_xy=0.0e0, + stress_shear_yz=0.0e0, + stress_shear_zx=0.0e0, + ) - svmyz = sigvm(sig_tf_r[ii], 0.0e0, sig_tf_z[ii], 0.0e0, 0.0e0, 0.0e0) + svmyz = calculate_von_mises_stress( + stress_x=sig_tf_r[ii], + stress_y=0.0e0, + stress_z=sig_tf_z[ii], + stress_shear_xy=0.0e0, + stress_shear_yz=0.0e0, + stress_shear_zx=0.0e0, + ) sig_tf_vmises[ii] = max(svmxz, svmyz) # Maximum shear stress for the Tresca yield criterion using CEA @@ -3702,44 +3721,6 @@ def eyoung_parallel( return eyoung_j_3, a_3, poisson_j_perp_3 -@numba.njit(cache=True) -def sigvm(sx: float, sy: float, sz: float, txy: float, txz: float, tyz: float) -> float: - """Calculates Von Mises stress in a TF coil - - This routine calculates the Von Mises combination of - stresses (Pa) in a TF coil. - - Parameters - ---------- - sx : - In-plane stress in X direction [Pa] - sy : - In-plane stress in Y direction [Pa] - sz : - In-plane stress in Z direction [Pa] - txy : - Out-of-plane stress in X-Y plane [Pa] - txz : - Out-of-plane stress in X-Z plane [Pa] - tyz : - Out-of-plane stress in Y-Z plane [Pa] - - Returns - ------- - : - Von Mises combination of stresses (Pa) in a TF coil. - """ - return np.sqrt( - 0.5 - * ( - (sx - sy) ** 2 - + (sx - sz) ** 2 - + (sz - sy) ** 2 - + 6 * (txy**2 + txz**2 + tyz**2) - ) - ) - - @numba.njit(cache=True, error_model="numpy") def extended_plane_strain( nu_t, diff --git a/tests/unit/models/engineering/test_materials.py b/tests/unit/models/engineering/test_materials.py index 6b05f78652..ec3f576900 100644 --- a/tests/unit/models/engineering/test_materials.py +++ b/tests/unit/models/engineering/test_materials.py @@ -1,6 +1,11 @@ +import numpy as np import pytest -from process.models.engineering.materials import eurofer97_thermal_conductivity +from process.models.engineering.materials import ( + calculate_tresca_stress, + calculate_von_mises_stress, + eurofer97_thermal_conductivity, +) def test_eurofer97_thermal_conductivity(monkeypatch, process_models): @@ -9,3 +14,147 @@ def test_eurofer97_thermal_conductivity(monkeypatch, process_models): assert eurofer97_thermal_conductivity( 1900.0, process_models.data.fwbs.fw_th_conductivity ) == pytest.approx(326.70406785462256) + + +def test_calculate_tresca_stress_floats(): + """Test Tresca stress calculation with float inputs.""" + result = calculate_tresca_stress(100.0, 50.0, 25.0) + assert result == pytest.approx(75.0) + + +def test_calculate_tresca_stress_ndarrays(): + """Test Tresca stress calculation with ndarray inputs.""" + stress_x = np.array([100.0, 200.0]) + stress_y = np.array([50.0, 100.0]) + stress_z = np.array([25.0, 50.0]) + + result = calculate_tresca_stress(stress_x, stress_y, stress_z) + expected = np.array([75.0, 150.0]) + + np.testing.assert_array_almost_equal(result, expected) + + +def test_calculate_tresca_stress_mixed_float_ndarray(): + """Test Tresca stress calculation with mixed float and ndarray inputs.""" + stress_x = 100.0 + stress_y = np.array([50.0, 75.0]) + stress_z = 25.0 + + result = calculate_tresca_stress(stress_x, stress_y, stress_z) + expected = np.array([75.0, 75.0]) + + np.testing.assert_array_almost_equal(result, expected) + + +def test_calculate_tresca_stress_zeros(): + """Test Tresca stress calculation with zero values.""" + result = calculate_tresca_stress(0.0, 0.0, 0.0) + assert result == pytest.approx(0.0) + + +def test_calculate_tresca_stress_zeros_ndarray(): + """Test Tresca stress calculation with zero ndarray values.""" + stress_x = np.array([0.0, 100.0]) + stress_y = np.array([0.0, 50.0]) + stress_z = np.array([0.0, 25.0]) + + result = calculate_tresca_stress(stress_x, stress_y, stress_z) + expected = np.array([0.0, 75.0]) + + np.testing.assert_array_almost_equal(result, expected) + + +def test_calculate_von_mises_stress_floats(): + """Test von Mises stress calculation with float inputs.""" + result = calculate_von_mises_stress(100.0, 50.0, 25.0, 10.0, 5.0, 2.0) + expected = np.sqrt( + 0.5 + * ( + (100.0 - 50.0) ** 2 + + (50.0 - 25.0) ** 2 + + (25.0 - 100.0) ** 2 + + 6 * (10.0**2 + 5.0**2 + 2.0**2) + ) + ) + assert result == pytest.approx(expected) + + +def test_calculate_von_mises_stress_ndarrays(): + """Test von Mises stress calculation with ndarray inputs.""" + stress_x = np.array([100.0, 200.0]) + stress_y = np.array([50.0, 100.0]) + stress_z = np.array([25.0, 50.0]) + stress_shear_xy = np.array([10.0, 20.0]) + stress_shear_yz = np.array([5.0, 10.0]) + stress_shear_zx = np.array([2.0, 4.0]) + + result = calculate_von_mises_stress( + stress_x, stress_y, stress_z, stress_shear_xy, stress_shear_yz, stress_shear_zx + ) + expected = np.sqrt( + 0.5 + * ( + (stress_x - stress_y) ** 2 + + (stress_y - stress_z) ** 2 + + (stress_z - stress_x) ** 2 + + 6 * (stress_shear_xy**2 + stress_shear_yz**2 + stress_shear_zx**2) + ) + ) + + np.testing.assert_array_almost_equal(result, expected) + + +def test_calculate_von_mises_stress_mixed_float_ndarray(): + """Test von Mises stress calculation with mixed float and ndarray inputs.""" + stress_x = 100.0 + stress_y = np.array([50.0, 75.0]) + stress_z = 25.0 + stress_shear_xy = 10.0 + stress_shear_yz = np.array([5.0, 7.5]) + stress_shear_zx = 2.0 + + result = calculate_von_mises_stress( + stress_x, stress_y, stress_z, stress_shear_xy, stress_shear_yz, stress_shear_zx + ) + expected = np.sqrt( + 0.5 + * ( + (stress_x - stress_y) ** 2 + + (stress_y - stress_z) ** 2 + + (stress_z - stress_x) ** 2 + + 6 * (stress_shear_xy**2 + stress_shear_yz**2 + stress_shear_zx**2) + ) + ) + + np.testing.assert_array_almost_equal(result, expected) + + +def test_calculate_von_mises_stress_zeros(): + """Test von Mises stress calculation with zero values.""" + result = calculate_von_mises_stress(0.0, 0.0, 0.0, 0.0, 0.0, 0.0) + assert result == pytest.approx(0.0) + + +def test_calculate_von_mises_stress_zeros_ndarray(): + """Test von Mises stress calculation with zero ndarray values.""" + stress_x = np.array([0.0, 100.0]) + stress_y = np.array([0.0, 50.0]) + stress_z = np.array([0.0, 25.0]) + stress_shear_xy = np.array([0.0, 10.0]) + stress_shear_yz = np.array([0.0, 5.0]) + stress_shear_zx = np.array([0.0, 2.0]) + + result = calculate_von_mises_stress( + stress_x, stress_y, stress_z, stress_shear_xy, stress_shear_yz, stress_shear_zx + ) + expected = np.sqrt( + 0.5 + * ( + (stress_x - stress_y) ** 2 + + (stress_y - stress_z) ** 2 + + (stress_z - stress_x) ** 2 + + 6 * (stress_shear_xy**2 + stress_shear_yz**2 + stress_shear_zx**2) + ) + ) + + np.testing.assert_array_almost_equal(result, expected) diff --git a/tests/unit/models/tfcoil/test_tfcoil.py b/tests/unit/models/tfcoil/test_tfcoil.py index da98774ee6..7e314c5e55 100644 --- a/tests/unit/models/tfcoil/test_tfcoil.py +++ b/tests/unit/models/tfcoil/test_tfcoil.py @@ -1950,21 +1950,6 @@ def test_eyoung_parallel_array(eyoungparallelarrayparam, monkeypatch): ) -@pytest.mark.parametrize( - ("sx", "sy", "sz", "expected"), - [ - (0, -3.2e8, 2.4e8, 486621002.42385757), - (-2.8e8, 0, 2.4e8, 450777106.7833858), - ], -) -def test_sigvm(sx, sy, sz, expected): - # could not find an example of a use in PROCESS where - # tx, ty, or tz were anything other than 0 - ret = tfcoil_module.sigvm(sx, sy, sz, 0, 0, 0) - - assert ret == pytest.approx(expected) - - @pytest.mark.parametrize( ( "ind_tf_coil",