diff --git a/CMakeLists.txt b/CMakeLists.txt index c44d8e64ff..8830394df4 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -88,7 +88,6 @@ LIST(APPEND PROCESS_SRCS heat_transport_variables.f90 buildings_variables.f90 water_usage_variables.f90 - constraint_equations.f90 constants.f90 build_variables.f90 current_drive_variables.f90 diff --git a/documentation/proc-pages/development/add-vars.md b/documentation/proc-pages/development/add-vars.md index 12852133c6..cc7457d864 100644 --- a/documentation/proc-pages/development/add-vars.md +++ b/documentation/proc-pages/development/add-vars.md @@ -187,78 +187,32 @@ After following the instruction to add an input variable, you can make the varia ## Add a constraint equation -Constraint equations are added to *PROCESS* in the following way: +Constraint equations are added to *PROCESS* in the `process/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. -1.

- Increment the parameter `ipeqns` in module `numerics` in the source file - `numerics.f90` in order to accommodate the new constraint. -

-2.

- Add a line to `lablcc` in the source file `numerics.f90` decribing the - constraint equation. -

-3.

- Add a line to the FORD description of `lablcc` the source file `numerics.f90`. -

-4.

- Add a new Fortran `case` statement to routine `CONSTRAINT_EQNS` in source - file `constraint_equations.f90`. -

-5.

- Then add a new subrountine including the `constraints` module ensuring that - all the variables used in the formula are contained in the modules specified - via `use, XX only: XX` statements. Use a similar formulation to that used - for the existing constraint equations, remembering that the code will try - to force `cc(i)` to be zero. -

-6.

- If the constraint is using a f-value, notify the constraint equation - number on the f-value description. -

- -

- Remember that if an inequality is being added, a new f-value iteration - variable may also need to be added to the code. -

+A constraint is simply added by registering the constraint to the manager using a decorator. -```fortran - do i = i1,i2 - ! The constraint value in physical units is - ! a) for consistency equations, the quantity to be equated, or - ! b) for limit equations, the limiting value. - ! The symbol is = for a consistency equation, < for an upper limit - ! or > for a lower limit. - select case (icc(i)) - ... - ! Equation for fusion power upper limit - case (9); call constraint_eqn_009(args) +```python +@ConstraintManager.register_constraint(1234, "m", "=") +def my_constraint_function(): ... ``` +The arguments to the `register_constraint` function are: -```fortran - subroutine constraint_eqn_009(args) - !! Equation for fusion power upper limit - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for fusion power upper limit - !! #=# physics - !! #=#=# fp_fusion_total_max_mw, p_fusion_total_max_mw - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! fp_fusion_total_max_mw : input real : f-value for maximum fusion power - !! p_fusion_total_max_mw : input real : maximum fusion power (MW) - !! p_fusion_total_mw : input real : fusion power (MW) - use constraint_variables, only: fp_fusion_total_max_mw, p_fusion_total_max_mw - use physics_variables, only: p_fusion_total_mw - implicit none - type (constraint_args_type), intent(out) :: args - - args%cc = 1.0D0 - fp_fusion_total_max_mw * p_fusion_total_max_mw/p_fusion_total_mw - args%con = p_fusion_total_max_mw * (1.0D0 - args%cc) - args%err = p_fusion_total_mw * args%cc - args%symbol = '<' - args%units = 'MW' - - end subroutine constraint_eqn_009 -``` +- Name (again, currently an integer) +- Unit (for output reporting purposes) +- Symbol (e.g. =, >=, <=. Again, for output reporting purposes) + + +`my_constraint_function` should be named appropriately and return a `ConstraintResult` which contains the: + +- Normalised residual error +- Constraint value +- Constraint error +```python +@ConstraintManager.register_constraint(1234, "m", "=") +def my_constraint_function(): + normalised_residual = ... + value = ... + error = ... + return ConstraintResult(normalised_residual, value, error) +``` \ No newline at end of file diff --git a/examples/single_model_evaluation.ipynb b/examples/single_model_evaluation.ipynb index abbee3cab6..7a46661367 100644 --- a/examples/single_model_evaluation.ipynb +++ b/examples/single_model_evaluation.ipynb @@ -238,6 +238,9 @@ "metadata": {}, "outputs": [], "source": [ + "from process.constraints import ConstraintManager\n", + "\n", + "\n", "def run_impurities(w_imp_fracs):\n", " \"\"\"Calculate responses to W impurities.\"\"\"\n", " n = w_imp_fracs.shape[0]\n", @@ -253,7 +256,7 @@ " single_run.models.physics.physics()\n", "\n", " # Evaluate constraint equation 15 (L-H threshold constraint)\n", - " con15_value, _, _, _, _ = process.fortran.constraints.constraint_eqn_015()\n", + " con15_value = ConstraintManager.evaluate_constraint(15).normalised_residual\n", "\n", " # Need to copy values\n", " p_plasma_rad_mw[i] = process.fortran.physics_variables.p_plasma_rad_mw.item()\n", @@ -372,7 +375,7 @@ ], "metadata": { "kernelspec": { - "display_name": "process", + "display_name": ".venv", "language": "python", "name": "python3" }, @@ -386,7 +389,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.9" + "version": "3.10.15" } }, "nbformat": 4, diff --git a/process/caller.py b/process/caller.py index 9721fa7367..3df4190087 100644 --- a/process/caller.py +++ b/process/caller.py @@ -7,6 +7,7 @@ import numpy as np from tabulate import tabulate +import process.constraints as constraints from process import fortran as ft from process.final import finalise from process.io.mfile import MFile @@ -75,7 +76,7 @@ def call_models(self, xc: np.ndarray, m: int) -> tuple[float, np.ndarray]: self._call_models_once(xc) # Evaluate objective function and constraints objf = objective_function(ft.numerics.minmax) - conf, _, _, _, _ = ft.constraints.constraint_eqns(m, -1) + conf, _, _, _, _ = constraints.constraint_eqns(m, -1) if objf_prev is None and conf_prev is None: # First run: run again to check idempotence diff --git a/process/constraints.py b/process/constraints.py new file mode 100644 index 0000000000..0fd0dd3789 --- /dev/null +++ b/process/constraints.py @@ -0,0 +1,2430 @@ +from collections.abc import Callable, Hashable +from dataclasses import dataclass +from typing import ClassVar, Literal + +import numpy as np + +import process.fortran as fortran +from process.exceptions import ProcessError, ProcessValueError + +ConstraintSymbolType = Literal["=", ">=", "<="] + + +@dataclass +class ConstraintResult: + """The constraint quantities given the current state of the code + (aka given an evaluation at the point x). + """ + + normalised_residual: float + """The normalised residual of the constraint.""" + constraint_value: float + """The value of the constraint (in the physical units).""" + constraint_error: float + """The residual error of the constraint (in the physical units).""" + + +@dataclass +class ConstraintRegistration: + """Contains the constraint equation and metadata about the constraint.""" + + name: Hashable + """The name (often a number) of the constraint. It can be any hashable e.g. a string.""" + constraint_equation: Callable[[], ConstraintResult] + """The constraint equation that, when called, returns the normalised resiudal, + constraint value, and constraint error. + """ + units: str + """The units of the constraint value and error.""" + symbol: ConstraintSymbolType + """The type of constraint (<=, >=, ==). Only used for writing output diagnostics, + this does not impact the calculations. + """ + + +class ConstraintManager: + """A singleton class that manages the registration of constraint equations + and metadata. + + This class maintains an internal registry of constraints indexed by their names. + Classmethods provide access to this registry or to directly evaluate constraints. + """ + + _constraint_registry: ClassVar[dict[Hashable, ConstraintRegistration]] = {} + """An internal registry of the PROCESS constraint equations""" + + def __init__(self): + raise NotImplementedError(f"{self.__class__.__name__} cannot be instantiated.") + + @classmethod + def num_constraints(cls): + """Return the number of constraints currently in the registry""" + return len(cls._constraint_registry) + + @classmethod + def register_constraint( + cls, name: Hashable, units: str, symbol: ConstraintSymbolType + ) -> Callable[[], Callable[[], ConstraintResult]]: + """A decorator to add a constraint equation with metadata to the registry. + + The decorator should wrap a function with no argument which returns a + ConstraintResult. + + :param name: the name of the constraint and how it can be indexed from the registry + :type name: Hashable + :param units: the units of the constraint written to the output files + :type units: str + :param symbol: the symbol of the constraint written to the output files + :type symbol: str + """ + + def wrapper(func: Callable[[], ConstraintResult]): + if name in cls._constraint_registry: + raise ValueError(f"Constraint {name} already exists.") + cls._constraint_registry[name] = ConstraintRegistration( + name, func, units, symbol + ) + + return func + + return wrapper + + @classmethod + def get_constraint(cls, name: Hashable): + """Retrieves a constraint registration from the registry given its name. + Returns None if no constraint with the name exists. + + :param name: the name of the constraint + :type name: Hashable + :returns: the constraint registration object + :rtype: ConstraintRegistration | None + """ + return cls._constraint_registry.get(name) + + @classmethod + def evaluate_constraint(cls, name: Hashable): + """Evalutes a constraint with a given name. + :param name: the name of the constraint + :type name: Hashable + :returns: the result of evaluating the constraint + :rtype: ConstraintResult | None + """ + registration = cls.get_constraint(name) + + if registration is not None: + return registration.constraint_equation() + + return None + + +@ConstraintManager.register_constraint(1, "", "=") +def constraint_equation_1(): + """Relationship between beta, temperature (keV) and density + + author: J Morris + + beta: total plasma beta + beta_{ft}: fast alpha beta component + beta_{NBI}: neutral beam beta component + n_e: electron density [/m3] + n_i: total ion density [/m3] + T_e: density weighted average electron temperature [keV] + T_i: density weighted average ion temperature [keV] + B_{tot}: total toroidal + poloidal field [T] + """ + cc = ( + 1.0 + - ( + fortran.physics_variables.beta_fast_alpha + + fortran.physics_variables.beta_beam + + 2.0e3 + * fortran.constants.rmu0 + * fortran.constants.electron_charge + * ( + fortran.physics_variables.dene * fortran.physics_variables.ten + + fortran.physics_variables.nd_ions_total + * fortran.physics_variables.tin + ) + / fortran.physics_variables.btot**2 + ) + / fortran.physics_variables.beta + ) + return ConstraintResult( + normalised_residual=cc, + constraint_value=(fortran.physics_variables.beta * (1.0 - cc)), + constraint_error=(fortran.physics_variables.beta * cc), + ) + + +@ConstraintManager.register_constraint(2, "MW/m3", "=") +def constraint_equation_2(): + """author: J. Morris + + i_rad_loss: switch for radiation loss term usage in power balance (see User Guide): + - 0 total power lost is scaling power plus radiation (needed for ipedestal=2,3) + - 1 total power lost is scaling power plus core radiation only + - 2 total power lost is scaling power only, with no additional + allowance for radiation. This is not recommended for power plant models. + + i_plasma_ignited: switch for ignition assumption: + - 0 do not assume plasma ignition; + - 1 assume ignited (but include auxiliary power in costs) + + pden_electron_transport_loss_mw: electron transport power per volume (MW/m3) + pden_ion_transport_loss_mw: ion transport power per volume (MW/m3) + pden_plasma_rad_mw: total radiation power per volume (MW/m3) + pden_plasma_core_rad_mw: total core radiation power per volume (MW/m3) + f_alpha_plasma: fraction of alpha power deposited in plasma + pden_alpha_total_mw: alpha power per volume (MW/m3) + pden_non_alpha_charged_mw: non-alpha charged particle fusion power per volume (MW/m3) + pden_plasma_ohmic_mw: ohmic heating power per volume (MW/m3) + p_hcd_injected_total_mw: total auxiliary injected power (MW) + vol_plasma: plasma volume (m3) + """ + # pscaling: total transport power per volume (MW/m3) + + pscaling = ( + fortran.physics_variables.pden_electron_transport_loss_mw + + fortran.physics_variables.pden_ion_transport_loss_mw + ) + # Total power lost is scaling power plus radiation: + if fortran.physics_variables.i_rad_loss == 0: + pnumerator = pscaling + fortran.physics_variables.pden_plasma_rad_mw + elif fortran.physics_variables.i_rad_loss == 1: + pnumerator = pscaling + fortran.physics_variables.pden_plasma_core_rad_mw + else: + pnumerator = pscaling + + # if plasma not ignited include injected power + if fortran.physics_variables.i_plasma_ignited == 0: + pdenom = ( + fortran.physics_variables.f_alpha_plasma + * fortran.physics_variables.pden_alpha_total_mw + + fortran.physics_variables.pden_non_alpha_charged_mw + + fortran.physics_variables.pden_plasma_ohmic_mw + + fortran.current_drive_variables.p_hcd_injected_total_mw + / fortran.physics_variables.vol_plasma + ) + else: + # if plasma ignited + pdenom = ( + fortran.physics_variables.f_alpha_plasma + * fortran.physics_variables.pden_alpha_total_mw + + fortran.physics_variables.pden_non_alpha_charged_mw + + fortran.physics_variables.pden_plasma_ohmic_mw + ) + + cc = 1.0 - pnumerator / pdenom + + return ConstraintResult(cc, pdenom * (1.0 - cc), pdenom * cc) + + +@ConstraintManager.register_constraint(3, "MW/m3", "=") +def constraint_equation_3(): + """Global power balance equation for ions + i_plasma_ignited: switch for ignition assumption + - 0 do not assume plasma ignition; + - 1 assume ignited (but include auxiliary power in costs) + + pden_ion_transport_loss_mw: ion transport power per volume (MW/m3) + piepv: ion/electron equilibration power per volume (MW/m3) + f_alpha_plasma: fraction of alpha power deposited in plasma + f_pden_alpha_ions_mw: alpha power per volume to ions (MW/m3) + p_hcd_injected_ions_mw: auxiliary injected power to ions (MW) + vol_plasma: plasma volume (m3) + """ + # No assume plasma ignition: + if fortran.physics_variables.i_plasma_ignited == 0: + cc = 1.0 - ( + fortran.physics_variables.pden_ion_transport_loss_mw + + fortran.physics_variables.piepv + ) / ( + fortran.physics_variables.f_alpha_plasma + * fortran.physics_variables.f_pden_alpha_ions_mw + + fortran.current_drive_variables.p_hcd_injected_ions_mw + / fortran.physics_variables.vol_plasma + ) + return ConstraintResult( + cc, + ( + fortran.physics_variables.f_alpha_plasma + * fortran.physics_variables.f_pden_alpha_ions_mw + + fortran.current_drive_variables.p_hcd_injected_ions_mw + / fortran.physics_variables.vol_plasma + ) + * (1.0 - cc), + ( + fortran.physics_variables.f_alpha_plasma + * fortran.physics_variables.f_pden_alpha_ions_mw + + fortran.current_drive_variables.p_hcd_injected_ions_mw + / fortran.physics_variables.vol_plasma + ) + * cc, + ) + + # Plasma ignited: + cc = 1.0 - ( + fortran.physics_variables.pden_ion_transport_loss_mw + + fortran.physics_variables.piepv + ) / ( + fortran.physics_variables.f_alpha_plasma + * fortran.physics_variables.f_pden_alpha_ions_mw + ) + return ConstraintResult( + cc, + ( + fortran.physics_variables.f_alpha_plasma + * fortran.physics_variables.f_pden_alpha_ions_mw + ) + * (1.0 - cc), + ( + fortran.physics_variables.f_alpha_plasma + * fortran.physics_variables.f_pden_alpha_ions_mw + ) + * cc, + ) + + +@ConstraintManager.register_constraint(4, "MW/m3", "=") +def constraint_equation_4(): + """Global power balance equation for electrons + author: P B Lloyd, CCFE, Culham Science Centre + + i_rad_loss: switch for radiation loss term usage in power balance + - 0 total power lost is scaling power plus radiation (needed for ipedestal=2,3) + - 1 total power lost is scaling power plus core radiation only + - 2 total power lost is scaling power only, with no additional + allowance for radiation. This is not recommended for power plant models. + + i_plasma_ignited: switch for ignition assumption: + - 0 do not assume plasma ignition; + - 1 assume ignited (but include auxiliary power in costs) + + pden_electron_transport_loss_mw: electron transport power per volume (MW/m3) + pden_plasma_rad_mw: total radiation power per volume (MW/m3) + pden_plasma_core_rad_mw: total core radiation power per volume (MW/m3) + f_alpha_plasma: fraction of alpha power deposited in plasma + f_pden_alpha_electron_mw: alpha power per volume to electrons (MW/m3) + piepv: ion/electron equilibration power per volume (MW/m3) + p_hcd_injected_electrons_mw: auxiliary injected power to electrons (MW) + vol_plasma: plasma volume (m3) + """ + # pscaling: total transport power per volume (MW/m3) + + pscaling = fortran.physics_variables.pden_electron_transport_loss_mw + # Total power lost is scaling power plus radiation: + if fortran.physics_variables.i_rad_loss == 0: + pnumerator = pscaling + fortran.physics_variables.pden_plasma_rad_mw + elif fortran.physics_variables.i_rad_loss == 1: + pnumerator = pscaling + fortran.physics_variables.pden_plasma_core_rad_mw + else: + pnumerator = pscaling + + # if plasma not ignited include injected power + if fortran.physics_variables.i_plasma_ignited == 0: + pdenom = ( + fortran.physics_variables.f_alpha_plasma + * fortran.physics_variables.f_pden_alpha_electron_mw + + fortran.physics_variables.piepv + + fortran.current_drive_variables.p_hcd_injected_electrons_mw + / fortran.physics_variables.vol_plasma + ) + else: + # if plasma ignited + pdenom = ( + fortran.physics_variables.f_alpha_plasma + * fortran.physics_variables.f_pden_alpha_electron_mw + + fortran.physics_variables.piepv + ) + + cc = 1.0 - pnumerator / pdenom + return ConstraintResult(cc, pdenom * (1.0 - cc), pdenom * cc) + + +@ConstraintManager.register_constraint(5, "/m3", "<=") +def constraint_equation_5(): + """Equation for density upper limit + author: P B Lloyd, CCFE, Culham Science Centre + + fdene: f-value for density limit + dene: electron density (/m3) + dnelimt: density limit (/m3) + dnla: line averaged electron density (m-3) + + i_density_limit: + - 1 old ASDEX; + - 2 Borrass model for ITER (I); + - 3 Borrass model for ITER (II); + - 4 JET edge radiation; + - 5 JET simplified; + - 6 Hugill-Murakami Mq limit; + - 7 Greenwald limit + """ + # Apply Greenwald limit to line-averaged density + if fortran.physics_variables.i_density_limit == 7: + return ConstraintResult( + fortran.physics_variables.dnla / fortran.physics_variables.dnelimt + - 1.0 * fortran.constraint_variables.fdene, + fortran.constraint_variables.fdene * fortran.physics_variables.dnelimt, + fortran.constraint_variables.fdene * fortran.physics_variables.dnelimt + - fortran.physics_variables.dnla, + ) + + cc = ( + fortran.physics_variables.dene / fortran.physics_variables.dnelimt + - 1.0 * fortran.constraint_variables.fdene + ) + return ConstraintResult( + cc, + fortran.physics_variables.dnelimt * (1.0 - cc), + fortran.physics_variables.dene * cc, + ) + + +@ConstraintManager.register_constraint(6, "", "<=") +def constraint_equation_6(): + """Equation for epsilon beta-poloidal upper limit + author: P B Lloyd, CCFE, Culham Science Centre + + fbeta_poloidal_eps: f-value for epsilon beta-poloidal + beta_poloidal_eps_max: maximum (eps*beta_poloidal) + eps: inverse aspect ratio + beta_poloidal: poloidal beta + """ + cc = ( + (fortran.physics_variables.eps * fortran.physics_variables.beta_poloidal) + / fortran.physics_variables.beta_poloidal_eps_max + - 1.0 * fortran.constraint_variables.fbeta_poloidal_eps + ) + return ConstraintResult( + cc, + fortran.physics_variables.beta_poloidal_eps_max * (1.0 - cc), + (fortran.physics_variables.eps * fortran.physics_variables.beta_poloidal) * cc, + ) + + +@ConstraintManager.register_constraint(7, "/m3", "=") +def constraint_equation_7(): + """Equation for hot beam ion density + + i_plasma_ignited: switch for ignition assumption: + - 0 do not assume plasma ignition + - 1 assume ignited (but include auxiliary power in costs) + O + bviously, i_plasma_ignited must be zero if current drive is required. + If i_plasma_ignited=1, any auxiliary power is assumed to be used only + during plasma start-up, and is excluded from all steady-state + power balance calculations. + beam_density_out: hot beam ion density from calculation (/m3) + nd_beam_ions: hot beam ion density, variable (/m3) + """ + if fortran.physics_variables.i_plasma_ignited == 1: + raise ProcessValueError( + "Do not use constraint equation 7 if i_plasma_ignited=1" + ) + + cc = ( + 1.0 + - fortran.physics_variables.beam_density_out + / fortran.physics_variables.nd_beam_ions + ) + return ConstraintResult( + cc, + fortran.physics_variables.nd_beam_ions * (1.0 - cc), + fortran.physics_variables.nd_beam_ions * cc, + ) + + +@ConstraintManager.register_constraint(8, "MW/m2", "<=") +def constraint_equation_8(): + """Equation for neutron wall load upper limit + + fwalld: f-value for maximum wall load + walalw: allowable wall-load (MW/m2) + pflux_fw_neutron_mw: average neutron wall load (MW/m2) + """ + return ConstraintResult( + ( + fortran.physics_variables.pflux_fw_neutron_mw + / fortran.constraint_variables.walalw + - 1.0 * fortran.constraint_variables.fwalld + ), + fortran.constraint_variables.fwalld * fortran.constraint_variables.walalw, + fortran.constraint_variables.fwalld * fortran.constraint_variables.walalw + - fortran.physics_variables.pflux_fw_neutron_mw, + ) + + +@ConstraintManager.register_constraint(9, "MW", "<=") +def constraint_equation_9(): + """Equation for fusion power upper limit + + fp_fusion_total_max_mw: f-value for maximum fusion power + p_fusion_total_max_mw: maximum fusion power (MW) + p_fusion_total_mw: fusion power (MW) + """ + cc = ( + fortran.physics_variables.p_fusion_total_mw + / fortran.constraint_variables.p_fusion_total_max_mw + - 1.0 * fortran.constraint_variables.fp_fusion_total_max_mw + ) + return ConstraintResult( + cc, + fortran.constraint_variables.p_fusion_total_max_mw * (1.0 - cc), + fortran.physics_variables.p_fusion_total_mw * cc, + ) + + +@ConstraintManager.register_constraint(11, "m", "=") +def constraint_equation_11(): + """Equation for radial build + author: P B Lloyd, CCFE, Culham Science Centre + + rbld: sum of thicknesses to the major radius (m) + rmajor: plasma major radius (m) + """ + cc = 1.0 - fortran.build_variables.rbld / fortran.physics_variables.rmajor + return ConstraintResult( + cc, + fortran.physics_variables.rmajor * (1.0 - cc), + fortran.physics_variables.rmajor * cc, + ) + + +@ConstraintManager.register_constraint(12, "V.sec", ">=") +def constraint_equation_12(): + """Equation for volt-second capability lower limit + author: P B Lloyd, CCFE, Culham Science Centre + + vs_plasma_total_required: total V-s needed (Wb) + vs_plasma_total_required (lower limit) is positive; vs_cs_pf_total_pulse (available) is negative + fvs: f-value for flux-swing (V-s) requirement (STEADY STATE) + vs_cs_pf_total_pulse: total flux swing for pulse (Wb) + """ + # vs_cs_pf_total_pulse is negative, requires sign change + cc = ( + 1.0 + - fortran.constraint_variables.fvs + * (-fortran.pfcoil_variables.vs_cs_pf_total_pulse) + / fortran.physics_variables.vs_plasma_total_required + ) + + return ConstraintResult( + cc, + fortran.pfcoil_variables.vs_plasma_total_required * (1.0 - cc), + fortran.pfcoil_variables.vs_plasma_total_required * cc, + ) + + +@ConstraintManager.register_constraint(13, "sec", ">=") +def constraint_equation_13(): + """Equation for burn time lower limit + + author: P B Lloyd, CCFE, Culham Science Centre + + ft_burn: f-value for minimum burn time + t_burn: burn time (s) (calculated if i_pulsed_plant=1) + t_burn_min: minimum burn time (s) + """ + return ConstraintResult( + 1.0 + - fortran.constraint_variables.ft_burn + * fortran.times_variables.t_burn + / fortran.constraint_variables.t_burn_min, + fortran.constraint_variables.t_burn_min / fortran.constraint_variables.ft_burn, + fortran.constraint_variables.t_burn_min / fortran.constraint_variables.ft_burn + - fortran.times_variables.t_burn, + ) + + +@ConstraintManager.register_constraint(15, "MW", ">=") +def constraint_equation_15(): + """Equation for L-H power threshold limit + author: P B Lloyd, CCFE, Culham Science Centre + + fl_h_threshold: f-value for L-H power threshold + p_l_h_threshold_mw: L-H mode power threshold (MW) + p_plasma_separatrix_mw: power to conducted to the divertor region (MW) + """ + return ConstraintResult( + 1.0 + - fortran.constraint_variables.fl_h_threshold + * fortran.physics_variables.p_plasma_separatrix_mw + / fortran.physics_variables.p_l_h_threshold_mw, + fortran.physics_variables.p_l_h_threshold_mw, + fortran.physics_variables.p_l_h_threshold_mw + - fortran.physics_variables.p_plasma_separatrix_mw + / fortran.constraint_variables.fl_h_threshold, + ) + + +@ConstraintManager.register_constraint(16, "MW", ">=") +def constraint_equation_16(): + """Equation for net electric power lower limit + author: P B Lloyd, CCFE, Culham Science Centre + + fpnetel: f-value for net electric power + p_plant_electric_net_mw: net electric power (MW) + pnetelin: required net electric power (MW) + """ + return ConstraintResult( + 1.0 + - fortran.constraint_variables.fpnetel + * fortran.heat_transport_variables.p_plant_electric_net_mw + / fortran.constraint_variables.pnetelin, + fortran.constraint_variables.pnetelin, + fortran.heat_transport_variables.p_plant_electric_net_mw + - fortran.constraint_variables.pnetelin / fortran.constraint_variables.fpnetel, + ) + + +@ConstraintManager.register_constraint(14, "", "=") +def constraint_equation_14(): + """Equation to fix number of NBI decay lengths to plasma centre + author: P B Lloyd, CCFE, Culham Science Centre + + n_beam_decay_lengths_core: neutral beam e-decay lengths to plasma centre + tbeamin: permitted neutral beam e-decay lengths to plasma centre + """ + cc = ( + 1.0 + - fortran.current_drive_variables.n_beam_decay_lengths_core + / fortran.current_drive_variables.tbeamin + ) + return ConstraintResult( + cc, + fortran.current_drive_variables.tbeamin * (1.0 - cc), + fortran.current_drive_variables.tbeamin * cc, + ) + + +@ConstraintManager.register_constraint(17, "MW/m3", "<=") +def constraint_equation_17(): + """Equation for radiation power upper limit + author: P B Lloyd, CCFE, Culham Science Centre + + f_alpha_plasma: fraction of alpha power deposited in plasma + p_hcd_injected_total_mw: total auxiliary injected power (MW) + vol_plasma: plasma volume (m3) + pden_alpha_total_mw: alpha power per volume (MW/m3) + pden_non_alpha_charged_mw: non-alpha charged particle fusion power per volume (MW/m3) + pden_plasma_ohmic_mw: ohmic heating power per volume (MW/m3) + fradpwr: f-value for core radiation power limit + pden_plasma_rad_mw: total radiation power per volume (MW/m3) + """ + # Maximum possible power/vol_plasma that can be radiated (local) + pradmaxpv = ( + fortran.current_drive_variables.p_hcd_injected_total_mw + / fortran.physics_variables.vol_plasma + + fortran.physics_variables.pden_alpha_total_mw + * fortran.physics_variables.f_alpha_plasma + + fortran.physics_variables.pden_non_alpha_charged_mw + + fortran.physics_variables.pden_plasma_ohmic_mw + ) + + cc = ( + fortran.physics_variables.pden_plasma_rad_mw / pradmaxpv + - 1.0 * fortran.constraint_variables.fradpwr + ) + return ConstraintResult( + cc, pradmaxpv * (1.0 - cc), fortran.physics_variables.pden_plasma_rad_mw * cc + ) + + +@ConstraintManager.register_constraint(18, "MW/m2", "<=") +def constraint_equation_18(): + """Equation for divertor heat load upper limit + author: P B Lloyd, CCFE, Culham Science Centre + + fpflux_div_heat_load_mw: f-value for divertor heat load + pflux_div_heat_load_max_mw: heat load limit (MW/m2) + pflux_div_heat_load_mw: divertor heat load (MW/m2) + """ + cc = ( + fortran.divertor_variables.pflux_div_heat_load_mw + / fortran.divertor_variables.pflux_div_heat_load_max_mw + - 1.0 * fortran.constraint_variables.fpflux_div_heat_load_mw + ) + return ConstraintResult( + cc, + fortran.divertor_variables.pflux_div_heat_load_max_mw * (1.0 - cc), + fortran.divertor_variables.pflux_div_heat_load_mw * cc, + ) + + +@ConstraintManager.register_constraint(19, "MVA", "<=") +def constraint_equation_19(): + """Equation for MVA (power) upper limit: resistive TF coil set + author: P B Lloyd, CCFE, Culham Science Centre + + p_cp_resistive_mw: peak resistive TF coil inboard leg power (total) (MW) + p_tf_leg_resistive_mw: TF coil outboard leg resistive power (total) (MW) + fmva: f-value for maximum MVA + mvalim: MVA limit for resistive TF coil set (total) (MW) + """ + totmva = ( + fortran.tfcoil_variables.p_cp_resistive_mw + + fortran.tfcoil_variables.p_tf_leg_resistive_mw + ) + + cc = ( + totmva / fortran.constraint_variables.mvalim + - 1.0 * fortran.constraint_variables.fmva + ) + return ConstraintManager( + cc, fortran.constraint_variables.mvalim * (1.0 - cc), totmva * cc + ) + + +@ConstraintManager.register_constraint(20, "m", "<=") +def constraint_equation_20(): + """Equation for neutral beam tangency radius upper limit + author: P B Lloyd, CCFE, Culham Science Centre + + fportsz: f-value for neutral beam tangency radius limit + rtanmax: maximum tangency radius for centreline of beam (m) + rtanbeam: neutral beam centreline tangency radius (m) + """ + cc = ( + fortran.current_drive_variables.rtanbeam + / fortran.current_drive_variables.rtanmax + - 1.0 * fortran.constraint_variables.fportsz + ) + return ConstraintResult( + cc, + fortran.current_drive_variables.rtanmax * (1.0 - cc), + fortran.current_drive_variables.rtanbeam * cc, + ) + + +@ConstraintManager.register_constraint(21, "", ">=") +def constraint_equation_21(): + """Equation for minor radius lower limit + author: P B Lloyd, CCFE, Culham Science Centre + + frminor: f-value for minor radius limit + rminor: plasma minor radius (m) + aplasmin: minimum minor radius (m) + """ + cc = ( + 1.0 + - fortran.constraint_variables.frminor + * fortran.physics_variables.rminor + / fortran.build_variables.aplasmin + ) + return ConstraintResult( + cc, + fortran.build_variables.aplasmin * (1.0 - cc), + fortran.build_variables.aplasmin * cc, + ) + + +@ConstraintManager.register_constraint(23, "m", "<=") +def constraint_equation_23(): + """Equation for conducting shell radius / rminor upper limit + author: P B Lloyd, CCFE, Culham Science Centre + + rminor: plasma minor radius (m) + dr_fw_plasma_gap_outboard: gap between plasma and first wall, outboard side (m) + dr_fw_outboard: outboard first wall thickness, initial estimate (m) + dr_blkt_outboard: outboard blanket thickness (m) + fr_conducting_wall: f-value for conducting wall radius / rminor limit + f_r_conducting_wall: maximum ratio of conducting wall distance to plasma minor radius for vertical stability + """ + # conducting shell radius (m) + rcw = ( + fortran.physics_variables.rminor + + fortran.build_variables.dr_fw_plasma_gap_outboard + + fortran.build_variables.dr_fw_outboard + + fortran.build_variables.dr_blkt_outboard + ) + + cc = ( + rcw + / ( + fortran.physics_variables.f_r_conducting_wall + * fortran.physics_variables.rminor + ) + - 1.0 * fortran.constraint_variables.fr_conducting_wall + ) + return ConstraintManager( + cc, + fortran.physics_variables.f_r_conducting_wall + * fortran.physics_variables.rminor + * (1.0 - cc), + rcw * cc, + ) + + +@ConstraintManager.register_constraint(24, "", "<=") +def constraint_equation_24(): + """Equation for beta upper limit + author: P B Lloyd, CCFE, Culham Science Centre + + i_beta_component: switch for beta limit scaling (constraint equation 24): + - 0 apply limit to total beta; + - 1 apply limit to thermal beta; + - 2 apply limit to thermal + neutral beam beta + - 3 apply limit to toroidal beta + istell: switch for stellarator option (set via device.dat): + - 0 use tokamak model; + - 1 use stellarator model + fbeta_max: f-value for beta limit + beta_max: allowable beta + beta: total plasma beta (calculated if ipedestal =3) + beta_fast_alpha: fast alpha beta component + beta_beam: neutral beam beta component + bt: toroidal field + btot: total field + """ + # Include all beta components: relevant for both tokamaks and stellarators + if ( + fortran.physics_variables.i_beta_component == 0 + or fortran.stellarator_variables.istell != 0 + ): + cc = ( + fortran.physics_variables.beta / fortran.physics_variables.beta_max + - 1.0 * fortran.constraint_variables.fbeta_max + ) + con = fortran.physics_variables.beta_max + err = ( + fortran.physics_variables.beta_max + - fortran.physics_variables.beta / fortran.constraint_variables.fbeta_max + ) + # Here, the beta limit applies to only the thermal component, not the fast alpha or neutral beam parts + elif fortran.physics_variables.i_beta_component == 1: + cc = ( + ( + fortran.physics_variables.beta + - fortran.physics_variables.beta_fast_alpha + - fortran.physics_variables.beta_beam + ) + / fortran.physics_variables.beta_max + - 1.0 * fortran.constraint_variables.fbeta_max + ) + con = fortran.physics_variables.beta_max + err = ( + fortran.physics_variables.beta_max + - ( + fortran.physics_variables.beta + - fortran.physics_variables.beta_fast_alpha + - fortran.physics_variables.beta_beam + ) + / fortran.constraint_variables.fbeta_max + ) + # Beta limit applies to thermal + neutral beam: components of the total beta, i.e. excludes alphas + elif fortran.physics_variables.i_beta_component == 2: + cc = ( + (fortran.physics_variables.beta - fortran.physics_variables.beta_fast_alpha) + / fortran.physics_variables.beta_max + - 1.0 * fortran.constraint_variables.fbeta_max + ) + con = fortran.physics_variables.beta_max * (1.0 - cc) + err = ( + fortran.physics_variables.beta - fortran.physics_variables.beta_fast_alpha + ) * cc + # Beta limit applies to toroidal beta + elif fortran.physics_variables.i_beta_component == 3: + cc = ( + ( + fortran.physics_variables.beta + * (fortran.physics_variables.btot / fortran.physics_variables.bt) ** 2 + ) + / fortran.physics_variables.beta_max + - 1.0 * fortran.constraint_variables.fbeta_max + ) + con = fortran.physics_variables.beta_max + err = ( + fortran.physics_variables.beta_max + - ( + fortran.physics_variables.beta + * (fortran.physics_variables.btot / fortran.physics_variables.bt) ** 2 + ) + / fortran.constraint_variables.fbeta_max + ) + + return ConstraintResult(cc, con, err) + + +@ConstraintManager.register_constraint(25, "T", "<=") +def constraint_equation_25(): + """Equation for peak toroidal field upper limit + author: P B Lloyd, CCFE, Culham Science Centre + + fpeakb: f-value for maximum toroidal field + bmxlim: maximum peak toroidal field (T) + b_tf_inboard_peak: mean peak field at TF coil (T) + """ + cc = ( + fortran.tfcoil_variables.b_tf_inboard_peak / fortran.constraint_variables.bmxlim + - 1.0 * fortran.constraint_variables.fpeakb + ) + return ConstraintResult( + cc, + fortran.constraint_variables.bmxlim * (1.0 - cc), + fortran.tfcoil_variables.b_tf_inboard_peak * cc, + ) + + +@ConstraintManager.register_constraint(26, "A/m2", "<=") +def constraint_equation_26(): + """Equation for Central Solenoid current density upper limit at EOF + author: P B Lloyd, CCFE, Culham Science Centre + + fjohcreal: f-value for central solenoid current at end-of-flattop + j_cs_critical_flat_top_end: allowable central solenoid current density at end of flat-top (A/m2) + j_cs_flat_top_end: central solenoid overall current density at end of flat-top (A/m2) + """ + return ConstraintResult( + fortran.pfcoil_variables.j_cs_flat_top_end + / fortran.pfcoil_variables.j_cs_critical_flat_top_end + - 1.0 * fortran.constraint_variables.fjohc, + fortran.pfcoil_variables.j_cs_critical_flat_top_end, + fortran.pfcoil_variables.j_cs_critical_flat_top_end + - fortran.pfcoil_variables.j_cs_flat_top_end + / fortran.constraint_variables.fjohc, + ) + + +@ConstraintManager.register_constraint(27, "A/m2", "<=") +def constraint_equation_27(): + """Equation for Central Solenoid current density upper limit at BOP + author: P B Lloyd, CCFE, Culham Science Centre + + fjohc0: f-value for central solenoid current at beginning of pulse + j_cs_critical_pulse_start: allowable central solenoid current density at beginning of pulse (A/m2) + j_cs_pulse_start: central solenoid overall current density at beginning of pulse (A/m2) + """ + return ConstraintResult( + fortran.pfcoil_variables.j_cs_pulse_start + / fortran.pfcoil_variables.j_cs_critical_pulse_start + - 1.0 * fortran.constraint_variables.fjohc0, + fortran.pfcoil_variables.j_cs_critical_pulse_start, + fortran.pfcoil_variables.j_cs_critical_pulse_start + - fortran.pfcoil_variables.j_cs_pulse_start + / fortran.constraint_variables.fjohc0, + ) + + +@ConstraintManager.register_constraint(28, "", ">=") +def constraint_equation_28(): + """Equation for fusion gain (big Q) lower limit + author: P B Lloyd, CCFE, Culham Science Centre + + fqval: pf-value for Q + bigq: Fusion gain; P_fusion / (P_injection + P_ohmic) + bigqmin: minimum fusion gain Q + i_plasma_ignited : input integer : switch for ignition assumption: + - 0 do not assume plasma ignition; + - 1 assume ignited (but include auxiliary power in costs) + Obviously, ignite must be zero if current drive is required. + If i_plasma_ignited=1, any auxiliary power is assumed to be used only + during plasma start-up, and is excluded from all steady-state + power balance calculations. + """ + if fortran.physics_variables.i_plasma_ignited != 0: + raise ProcessValueError("Do not use constraint 28 if i_plasma_ignited=1") + + cc = ( + 1.0 + - fortran.constraint_variables.fqval + * fortran.current_drive_variables.bigq + / fortran.constraint_variables.bigqmin + ) + return ConstraintResult( + cc, + fortran.constraint_variables.bigqmin * (1.0 - cc), + fortran.constraint_variables.bigqmin * cc, + ) + + +@ConstraintManager.register_constraint(29, "m", "=") +def constraint_equation_29(): + """Equation for inboard major radius: This is a consistency equation + author: P B Lloyd, CCFE, Culham Science Centre + + rmajor: plasma major radius (m) (iteration variable 3) + rminor: plasma minor radius (m) + rinboard: plasma inboard radius (m) + """ + cc = ( + 1.0 + - (fortran.physics_variables.rmajor - fortran.physics_variables.rminor) + / fortran.build_variables.rinboard + ) + return ConstraintResult( + cc, + fortran.build_variables.rinboard * (1.0 - cc), + fortran.build_variables.rinboard * cc, + ) + + +@ConstraintManager.register_constraint(30, "MW", "<=") +def constraint_equation_30(): + """Equation for injection power upper limit + author: P B Lloyd, CCFE, Culham Science Centre + + p_hcd_injected_total_mw: total auxiliary injected power (MW) + fpinj: f-value for injection power + p_hcd_injected_max: Maximum allowable value for injected power (MW) + """ + return ConstraintResult( + fortran.current_drive_variables.p_hcd_injected_total_mw + / fortran.current_drive_variables.p_hcd_injected_max + - 1.0 * fortran.constraint_variables.fpinj, + fortran.current_drive_variables.p_hcd_injected_max, + fortran.current_drive_variables.p_hcd_injected_max + - fortran.current_drive_variables.p_hcd_injected_total_mw + / fortran.constraint_variables.fpinj, + ) + + +@ConstraintManager.register_constraint(31, "Pa", "<=") +def constraint_equation_31(): + """Equation for TF coil case stress upper limit (SCTF) + author: P B Lloyd, CCFE, Culham Science Centre + + fstrcase: f-value for TF coil case stress + sig_tf_case_max: Allowable maximum shear stress in TF coil case (Tresca criterion) (Pa) + sig_tf_case: Constrained stress in TF coil case (Pa) + """ + return ConstraintResult( + fortran.tfcoil_variables.sig_tf_case / fortran.tfcoil_variables.sig_tf_case_max + - 1.0 * fortran.constraint_variables.fstrcase, + fortran.tfcoil_variables.sig_tf_case_max, + fortran.tfcoil_variables.sig_tf_case_max + - fortran.tfcoil_variables.sig_tf_case / fortran.constraint_variables.fstrcase, + ) + + +@ConstraintManager.register_constraint(32, "Pa", "<=") +def constraint_equation_32(): + """Equation for TF coil conduit stress upper limit (SCTF) + author: P B Lloyd, CCFE, Culham Science Centre + + fstrcond: f-value for TF coil conduit stress + sig_tf_wp_max: Allowable maximum shear stress in TF coil conduit (Tresca criterion) (Pa) + sig_tf_wp: Constrained stress in TF conductor conduit (Pa) + """ + return ConstraintResult( + fortran.tfcoil_variables.sig_tf_wp / fortran.tfcoil_variables.sig_tf_wp_max + - 1.0 * fortran.constraint_variables.fstrcond, + fortran.tfcoil_variables.sig_tf_wp_max, + fortran.tfcoil_variables.sig_tf_wp_max + - fortran.tfcoil_variables.sig_tf_wp / fortran.constraint_variables.fstrcond, + ) + + +@ConstraintManager.register_constraint(33, "A/m2", "<=") +def constraint_equation_33(): + """Equation for TF coil operating/critical J upper limit (SCTF) + author: P B Lloyd, CCFE, Culham Science Centre + args : output structure : residual error; constraint value; + + fiooic: f-value for TF coil operating current / critical + jwdgcrt: critical current density for winding pack (A/m2) + j_tf_wp: winding pack current density (A/m2) + """ + if fortran.constraint_variables.fiooic > 0.7: + fortran.error_handling.report_error(285) + + cc = ( + fortran.tfcoil_variables.j_tf_wp / fortran.tfcoil_variables.jwdgcrt + - 1.0 * fortran.constraint_variables.fiooic + ) + return ConstraintResult( + cc, + fortran.tfcoil_variables.jwdgcrt * (1.0 - cc), + fortran.tfcoil_variables.j_tf_wp * cc, + ) + + +@ConstraintManager.register_constraint(34, "V", "<=") +def constraint_equation_34(): + """Equation for TF coil dump voltage upper limit (SCTF) + author: P B Lloyd, CCFE, Culham Science Centre + + fvdump: f-value for dump voltage + vdalw: max voltage across TF coil during quench (kV) + vtfskv: voltage across a TF coil during quench (kV) + """ + return ConstraintResult( + fortran.tfcoil_variables.vtfskv / fortran.tfcoil_variables.vdalw + - 1.0 * fortran.constraint_variables.fvdump, + fortran.tfcoil_variables.vdalw, + fortran.tfcoil_variables.vdalw - fortran.tfcoil_variables.vtfskv, + ) + + +@ConstraintManager.register_constraint(35, "A/m2", "<=") +def constraint_equation_35(): + """Equation for TF coil J_wp/J_prot upper limit (SCTF) + author: P B Lloyd, CCFE, Culham Science Centre + + fjprot: f-value for TF coil winding pack current density + jwdgpro: allowable TF coil winding pack current density, for dump temperature + rise protection (A/m2) + j_tf_wp: winding pack current density (A/m2) + """ + return ConstraintResult( + fortran.tfcoil_variables.j_tf_wp / fortran.tfcoil_variables.jwdgpro + - 1.0 * fortran.constraint_variables.fjprot, + fortran.tfcoil_variables.jwdgpro, + fortran.tfcoil_variables.j_tf_wp - fortran.tfcoil_variables.jwdgpro, + ) + + +@ConstraintManager.register_constraint(36, "K", ">=") +def constraint_equation_36(): + """Equation for TF coil s/c temperature margin lower limit (SCTF) + author: P B Lloyd, CCFE, Culham Science Centre + + ftmargtf: f-value for TF coil temperature margin + tmargtf: TF coil temperature margin (K) + tmargmin_tf: minimum allowable temperature margin : TF coils (K) + """ + return ConstraintResult( + 1.0 + - fortran.constraint_variables.ftmargtf + * fortran.tfcoil_variables.tmargtf + / fortran.tfcoil_variables.tmargmin_tf, + fortran.tfcoil_variables.tmargmin_tf, + fortran.tfcoil_variables.tmargmin_tf - fortran.tfcoil_variables.tmargtf, + ) + + +@ConstraintManager.register_constraint(37, "1E20 A/Wm2", "<=") +def constraint_equation_37(): + """Equation for current drive gamma upper limit + author: P B Lloyd, CCFE, Culham Science Centre + + fgamcd: f-value for current drive gamma + gammax: maximum current drive gamma + eta_cd_norm_hcd_primary: normalised current drive efficiency (1.0e20 A/W-m2) + """ + cc = ( + fortran.current_drive_variables.eta_cd_norm_hcd_primary + / fortran.constraint_variables.gammax + - 1.0 * fortran.constraint_variables.fgamcd + ) + return ConstraintResult( + cc, + fortran.constraint_variables.gammax * (1.0 - cc), + fortran.current_drive_variables.eta_cd_norm_hcd_primary * cc, + ) + + +@ConstraintManager.register_constraint(39, "K", "<=") +def constraint_equation_39(): + """Equation for first wall temperature upper limit + author: P B Lloyd, CCFE, Culham Science Centre + + ftpeak: f-value for first wall peak temperature + temp_fw_max: maximum temperature of first wall material (K) (i_thermal_electric_conversion>1) + temp_fw_peak: peak first wall temperature (K) + """ + if fortran.fwbs_variables.temp_fw_peak < 1.0: + raise ProcessValueError( + "temp_fw_peak = 0 implies i_pulsed_plant=0; do not use constraint 39 if i_pulsed_plant=0" + ) + cc = ( + fortran.fwbs_variables.temp_fw_peak / fortran.fwbs_variables.temp_fw_max + - 1.0 * fortran.constraint_variables.ftpeak + ) + return ConstraintResult( + cc, + fortran.fwbs_variables.temp_fw_max * (1.0 - cc), + fortran.fwbs_variables.temp_fw_peak * cc, + ) + + +@ConstraintManager.register_constraint(40, "MW", ">=") +def constraint_equation_40(): + """Equation for auxiliary power lower limit + author: P B Lloyd, CCFE, Culham Science Centre + + fauxmn: f-value for minimum auxiliary power + p_hcd_injected_total_mw: total auxiliary injected power (MW) + auxmin: minimum auxiliary power (MW) + """ + cc = ( + 1.0 + - fortran.constraint_variables.fauxmn + * fortran.current_drive_variables.p_hcd_injected_total_mw + / fortran.constraint_variables.auxmin + ) + return ConstraintResult( + cc, + fortran.constraint_variables.auxmin * (1.0 - cc), + fortran.constraint_variables.auxmin * cc, + ) + + +@ConstraintManager.register_constraint(41, "sec", ">=") +def constraint_equation_41(): + """Equation for plasma current ramp-up time lower limit + author: P B Lloyd, CCFE, Culham Science Centre + + ft_current_ramp_up: f-value for plasma current ramp-up time + t_current_ramp_up: plasma current ramp-up time for current initiation (s) + t_current_ramp_up_min: minimum plasma current ramp-up time (s) + """ + cc = ( + 1.0 + - fortran.constraint_variables.ft_current_ramp_up + * fortran.times_variables.t_current_ramp_up + / fortran.constraint_variables.t_current_ramp_up_min + ) + return ConstraintResult( + cc, + fortran.constraint_variables.t_current_ramp_up_min * (1.0 - cc), + fortran.constraint_variables.t_current_ramp_up_min * cc, + ) + + +@ConstraintManager.register_constraint(42, "sec", ">=") +def constraint_equation_42(): + """Equation for cycle time lower limit + author: P B Lloyd, CCFE, Culham Science Centre + + ftcycl: f-value for cycle time + t_cycle: full cycle time (s) + tcycmn: minimum cycle time (s) + """ + if fortran.constraint_variables.tcycmn < 1.0: + raise ProcessValueError( + "tcycmn = 0 implies that i_pulsed_plant=0; do not use constraint 42 if i_pulsed_plant=0" + ) + + cc = ( + 1.0 + - fortran.constraint_variables.ftcycl + * fortran.times_variables.t_cycle + / fortran.constraint_variables.tcycmn + ) + return ConstraintResult( + cc, + fortran.constraint_variables.tcycmn * (1.0 - cc), + fortran.constraint_variables.tcycmn * cc, + ) + + +@ConstraintManager.register_constraint(43, "deg C", "=") +def constraint_equation_43(): + """Equation for average centrepost temperature: This is a consistency equation (TART) + author: P B Lloyd, CCFE, Culham Science Centre + + temp_cp_average: average temp of TF coil inboard leg conductor (C)e + tcpav2: centrepost average temperature (C) (for consistency) + itart: switch for spherical tokamak (ST) models: + - 0 use conventional aspect ratio models; + - 1 use spherical tokamak models + """ + if fortran.physics_variables.itart == 0: + raise ProcessValueError("Do not use constraint 43 if itart=0") + + if fortran.tfcoil_variables.i_tf_sup == 0: + temp_cp_average = fortran.tfcoil_variables.temp_cp_average - 273.15 + tcpav2 = fortran.tfcoil_variables.tcpav2 - 273.15 + else: + temp_cp_average = fortran.tfcoil_variables.temp_cp_average + tcpav2 = fortran.tfcoil_variables.tcpav2 + + cc = 1.0 - temp_cp_average / tcpav2 + + return ConstraintResult(cc, tcpav2 * (1.0 - cc), tcpav2 * cc) + + +@ConstraintManager.register_constraint(44, "deg C", "<=") +def constraint_equation_44(): + """Equation for centrepost temperature upper limit (TART) + author: P B Lloyd, CCFE, Culham Science Centre + + fptemp: f-value for peak centrepost temperature + ptempalw: maximum peak centrepost temperature (K) + tcpmax: peak centrepost temperature (K) + itart: switch for spherical tokamak (ST) models: + - 0: use conventional aspect ratio models; + - 1: use spherical tokamak models + """ + if fortran.physics_variables.itart == 0: + raise ProcessValueError("Do not use constraint 44 if itart=0") + + if fortran.tfcoil_variables.i_tf_sup == 0: # ! Copper case + ptempalw = fortran.tfcoil_variables.ptempalw - 273.15 + tcpmax = fortran.tfcoil_variables.tcpmax - 273.15 + else: + ptempalw = fortran.tfcoil_variables.ptempalw + tcpmax = fortran.tfcoil_variables.tcpmax + + cc = tcpmax / ptempalw - 1.0 * fortran.constraint_variables.fptemp + return ConstraintResult(cc, ptempalw * (1.0 - cc), tcpmax * cc) + + +@ConstraintManager.register_constraint(45, "", ">=") +def constraint_manager_45(): + """Equation for edge safety factor lower limit (TART) + author: P B Lloyd, CCFE, Culham Science Centre + + fq: f-value for edge safety factor + q95 : safety factor 'near' plasma edge + (unless i_plasma_current = 2 (ST current scaling), in which case q = mean edge safety factor qbar) + q95_min: lower limit for edge safety factor + itart : input integer : switch for spherical tokamak (ST) models: + - 0 use conventional aspect ratio models; + - 1 use spherical tokamak models""" + if fortran.physics_variables.itart == 0: + raise ProcessValueError("Do not use constraint 45 if itart=0") + + cc = ( + 1.0 + - fortran.constraint_variables.fq + * fortran.physics_variables.q95 + / fortran.physics_variables.q95_min + ) + return ConstraintResult( + cc, + fortran.physics_variables.q95_min * (1.0 - cc), + fortran.physics_variables.q95_min * cc, + ) + + +@ConstraintManager.register_constraint(46, "", "<=") +def constraint_equation_46(): + """Equation for Ip/Irod upper limit (TART) + author: P B Lloyd, CCFE, Culham Science Centre + + eps: inverse aspect ratio + fipir: f-value for Ip/Irod upper limit + c_tf_total: total (summed) current in TF coils (A) + plasma_current: plasma current (A) + itart: switch for spherical tokamak (ST) models: + - 0: use conventional aspect ratio models; + - 1: use spherical tokamak models + """ + if fortran.physics_variables.itart == 0: + raise ProcessValueError("Do not use constraint 46 if itart=0") + + # maximum ratio of plasma current to centrepost current + cratmx = 1.0 + 4.91 * (fortran.physics_variables.eps - 0.62) + cc = ( + fortran.physics_variables.plasma_current / fortran.tfcoil_variables.c_tf_total + ) / cratmx - 1.0 * fortran.constraint_variables.fipir + + return ConstraintResult( + cc, + cratmx * (1.0 - cc), + fortran.physics_variables.plasma_current + / fortran.tfcoil_variables.c_tf_total + * cc, + ) + + +@ConstraintManager.register_constraint(48, "", "<=") +def constraint_equation_48(): + """Equation for poloidal beta upper limit + author: P B Lloyd, CCFE, Culham Science Centre + + fbeta_poloidal: rf-value for poloidal beta + beta_poloidal_max: maximum poloidal beta + beta_poloidal: poloidal beta + """ + cc = ( + fortran.physics_variables.beta_poloidal + / fortran.constraint_variables.beta_poloidal_max + - 1.0 * fortran.constraint_variables.fbeta_poloidal + ) + return ConstraintResult( + cc, + fortran.constraint_variables.beta_poloidal_max * (1.0 - cc), + fortran.physics_variables.beta_poloidal * cc, + ) + + +@ConstraintManager.register_constraint(50, "Hz", "<=") +def constraint_equation_50(): + """IFE option: Equation for repetition rate upper limit + author: P B Lloyd, CCFE, Culham Science Centre + author: S I Muldrew, CCFE, Culham Science Centre + """ + cc = ( + fortran.ife_variables.reprat / fortran.ife_variables.rrmax + - 1.0 * fortran.ife_variables.frrmax + ) + return ConstraintResult( + cc, fortran.ife_variables.rrmax * (1.0 - cc), fortran.ife_variables.reprat * cc + ) + + +@ConstraintManager.register_constraint(51, "V.s", "=") +def constraint_equation_51(): + """Equation to enforce startup flux = available startup flux + author: P B Lloyd, CCFE, Culham Science Centre + + vs_plasma_res_ramp: resistive losses in startup V-s (Wb) + vs_plasma_ind_ramp: internal and external plasma inductance V-s (Wb)) + vs_cs_pf_total_ramp: total flux swing for startup (Wb) + """ + cc = 1.0 - fortran.pfcoil_variables.fvs_cs_pf_total_ramp * abs( + ( + fortran.physics_variables.vs_plasma_res_ramp + + fortran.physics_variables.vs_plasma_ind_ramp + ) + / fortran.pfcoil_variables.vs_cs_pf_total_ramp + ) + return ConstraintResult( + cc, + fortran.pfcoil_variables.vs_cs_pf_total_ramp * (1.0 - cc), + fortran.pfcoil_variables.vs_cs_pf_total_ramp * cc, + ) + + +@ConstraintManager.register_constraint(52, "", ">=") +def constraint_equation_52(): + """Equation for tritium breeding ratio lower limit + author: P B Lloyd, CCFE, Culham Science Centre + + ftbr: f-value for minimum tritium breeding ratio + tbr: tritium breeding ratio (i_blanket_type=2,3 (KIT HCPB/HCLL)) + tbrmin: minimum tritium breeding ratio (If i_blanket_type=1, tbrmin=minimum 5-year time-averaged tritium breeding ratio) + """ + cc = ( + 1.0 + - fortran.constraint_variables.ftbr + * fortran.fwbs_variables.tbr + / fortran.constraint_variables.tbrmin + ) + return ConstraintResult( + cc, + fortran.constraint_variables.tbrmin * (1.0 - cc), + fortran.constraint_variables.tbrmin * cc, + ) + + +@ConstraintManager.register_constraint(53, "neutron/m2", "<=") +def constraint_equation_53(): + """Equation for fast neutron fluence on TF coil upper limit + author: P B Lloyd, CCFE, Culham Science Centre + + fflutf: f-value for maximum TF coil nuclear heating + nflutfmax: max fast neutron fluence on TF coil (n/m2) + nflutf: peak fast neutron fluence on TF coil superconductor (n/m2) + """ + cc = ( + fortran.fwbs_variables.nflutf / fortran.constraint_variables.nflutfmax + - 1.0 * fortran.constraint_variables.fflutf + ) + return ConstraintResult( + cc, + fortran.constraint_variables.nflutfmax * (1.0 - cc), + fortran.fwbs_variables.nflutf * cc, + ) + + +@ConstraintManager.register_constraint(54, "MW/m3", "<=") +def constraint_equation_54(): + """Equation for peak TF coil nuclear heating upper limit + author: P B Lloyd, CCFE, Culham Science Centre + + fptfnuc: f-value for maximum TF coil nuclear heating + ptfnucmax: maximum nuclear heating in TF coil (MW/m3) + ptfnucpm3: nuclear heating in the TF coil (MW/m3) (blktmodel>0) + """ + cc = ( + fortran.fwbs_variables.ptfnucpm3 / fortran.constraint_variables.ptfnucmax + - 1.0 * fortran.constraint_variables.fptfnuc + ) + return ConstraintResult( + cc, + fortran.constraint_variables.ptfnucmax * (1.0 - cc), + fortran.fwbs_variables.ptfnucpm3 * cc, + ) + + +@ConstraintManager.register_constraint(56, "MW/m", "<=") +def constraint_equation_56(): + """Equation for power through separatrix / major radius upper limit + author: P B Lloyd, CCFE, Culham Science Centre + + fpsepr: f-value for maximum Psep/R limit + pseprmax: maximum ratio of power crossing the separatrix to plasma major radius (Psep/R) (MW/m) + p_plasma_separatrix_mw: power to be conducted to the divertor region (MW) + rmajor: plasma major radius (m) + """ + cc = ( + ( + fortran.physics_variables.p_plasma_separatrix_mw + / fortran.physics_variables.rmajor + ) + / fortran.constraint_variables.pseprmax + - 1.0 * fortran.constraint_variables.fpsepr + ) + return ConstraintResult( + cc, + fortran.constraint_variables.pseprmax * (1.0 - cc), + ( + fortran.physics_variables.p_plasma_separatrix_mw + / fortran.physics_variables.rmajor + ) + * cc, + ) + + +@ConstraintManager.register_constraint(59, "", "<=") +def constraint_equation_59(): + """Equation for neutral beam shine-through fraction upper limit + author: P B Lloyd, CCFE, Culham Science Centre + + fnbshinef: f-value for maximum neutral beam shine-through fraction + f_p_beam_shine_through_max: maximum neutral beam shine-through fraction + f_p_beam_shine_through: neutral beam shine-through fraction + """ + cc = ( + fortran.current_drive_variables.f_p_beam_shine_through + / fortran.constraint_variables.f_p_beam_shine_through_max + - 1.0 * fortran.constraint_variables.fnbshinef + ) + return ConstraintResult( + cc, + fortran.constraint_variables.f_p_beam_shine_through_max * (1.0 - cc), + fortran.current_drive_variables.f_p_beam_shine_through * cc, + ) + + +@ConstraintManager.register_constraint(60, "K", ">=") +def constraint_equation_60(): + """Equation for Central Solenoid s/c temperature margin lower limit + author: P B Lloyd, CCFE, Culham Science Centre + + ftmargoh: f-value for central solenoid temperature margin + temp_cs_margin: Central solenoid temperature margin (K) + tmargmin_cs: Minimum allowable temperature margin : CS (K) + """ + return ConstraintResult( + 1.0 + - fortran.constraint_variables.ftmargoh + * fortran.pfcoil_variables.temp_cs_margin + / fortran.tfcoil_variables.tmargmin_cs, + fortran.tfcoil_variables.tmargmin_cs, + fortran.tfcoil_variables.tmargmin_cs - fortran.pfcoil_variables.temp_cs_margin, + ) + + +@ConstraintManager.register_constraint(61, "", ">=") +def constraint_equation_61(): + """Equation for availability lower limit + author: P B Lloyd, CCFE, Culham Science Centre + + favail: F-value for minimum availability + cfactr: Total plant availability fraction + avail_min: Minimum availability + """ + cc = ( + 1.0 + - fortran.cost_variables.favail + * fortran.cost_variables.cfactr + / fortran.cost_variables.avail_min + ) + return ConstraintResult( + cc, + fortran.cost_variables.avail_min * (1.0 - cc), + fortran.cost_variables.cfactr * cc, + ) + + +@ConstraintManager.register_constraint(62, "", ">=") +def constraint_equation_62(): + """Lower limit on f_alpha_energy_confinement the ratio of alpha particle to energy confinement times + author: P B Lloyd, CCFE, Culham Science Centre + + falpha_energy_confinement: f-value for lower limit on f_alpha_energy_confinement the ratio of alpha particle to energy confinement + t_alpha_confinement: alpha particle confinement time (s) + t_energy_confinement: global thermal energy confinement time (sec) + f_alpha_energy_confinement_min: Lower limit on f_alpha_energy_confinement the ratio of alpha particle to energy confinement times + """ + cc = ( + 1.0 + - fortran.constraint_variables.falpha_energy_confinement + * ( + fortran.physics_variables.t_alpha_confinement + / fortran.physics_variables.t_energy_confinement + ) + / fortran.constraint_variables.f_alpha_energy_confinement_min + ) + return ConstraintResult( + cc, + fortran.constraint_variables.f_alpha_energy_confinement_min, + ( + fortran.physics_variables.t_alpha_confinement + / fortran.physics_variables.t_energy_confinement + ) + * cc, + ) + + +@ConstraintManager.register_constraint(63, "", "<=") +def constraint_equation_63(): + """Upper limit on niterpump (vacuum_model = simple) + author: P B Lloyd, CCFE, Culham Science Centre + + fniterpump: f-value for constraint that number of pumps < tfno + tfno: number of TF coils (default = 50 for stellarators) + niterpump: number of high vacuum pumps (real number), each with the throughput + """ + cc = ( + fortran.vacuum_variables.niterpump / fortran.tfcoil_variables.n_tf_coils + - 1.0 * fortran.constraint_variables.fniterpump + ) + return ConstraintResult( + cc, + fortran.tfcoil_variables.n_tf_coils, + fortran.tfcoil_variables.n_tf_coils * cc, + ) + + +@ConstraintManager.register_constraint(64, "", "<=") +def constraint_equation_64(): + """Upper limit on Zeff + author: P B Lloyd, CCFE, Culham Science Centre + + fzeffmax: f-value for maximum zeff + zeffmax: maximum value for Zeff + zeff: plasma effective charge + """ + cc = ( + fortran.physics_variables.zeff / fortran.constraint_variables.fzeffmax + - 1.0 * fortran.constraint_variables.ffzeffmax + ) + return ConstraintResult( + cc, + fortran.constraint_variables.fzeffmax, + fortran.constraint_variables.fzeffmax * cc, + ) + + +@ConstraintManager.register_constraint(65, "Pa", "<=") +def constraint_equation_65(): + """Upper limit on stress of the vacuum vessel that occurs when the TF coil quenches. + author: Timothy Nunn, UKAEA + + fmaxvvstress: f-value for constraint on maximum VV stress + max_vv_stress: Maximum permitted stress of the VV (Pa) + vv_stress_quench: Stress of the VV (Pa) + """ + cc = ( + fortran.sctfcoil_module.vv_stress_quench + / fortran.tfcoil_variables.max_vv_stress + - 1.0 * fortran.constraint_variables.fmaxvvstress + ) + return ConstraintResult( + cc, + fortran.tfcoil_variables.max_vv_stress, + fortran.tfcoil_variables.max_vv_stress * cc, + ) + + +@ConstraintManager.register_constraint(66, "MW", "<=") +def constrain_equation_66(): + """Upper limit on rate of change of energy in poloidal field + author: P B Lloyd, CCFE, Culham Science Centre + + fpoloidalpower: f-value for constraint on rate of change of energy in poloidal field + maxpoloidalpower: Maximum permitted absolute rate of change of stored energy in poloidal field (MW) + peakpoloidalpower: Peak absolute rate of change of stored energy in poloidal field (MW) (11/01/16) + """ + cc = ( + fortran.pf_power_variables.peakpoloidalpower + / fortran.pf_power_variables.maxpoloidalpower + - 1.0 * fortran.constraint_variables.fpoloidalpower + ) + return ConstraintResult( + cc, + fortran.pf_power_variables.maxpoloidalpower, + fortran.pf_power_variables.maxpoloidalpower * cc, + ) + + +@ConstraintManager.register_constraint(67, "MW/m2", "<=") +def constraint_equation_67(): + """Simple upper limit on radiation wall load + author: P B Lloyd, CCFE, Culham Science Centre + + fradwall: f-value for upper limit on radiation wall load + pflux_fw_rad_max: Maximum permitted radiation wall load (MW/m^2) + pflux_fw_rad_max_mw: Peak radiation wall load (MW/m^2) + """ + cc = ( + fortran.constraint_variables.pflux_fw_rad_max_mw + / fortran.constraint_variables.pflux_fw_rad_max + - 1.0 * fortran.constraint_variables.fradwall + ) + return ConstraintResult( + cc, + fortran.constraint_variables.pflux_fw_rad_max, + fortran.constraint_variables.pflux_fw_rad_max * cc, + ) + + +@ConstraintManager.register_constraint(68, "MWT/m", "<=") +def constraint_equation_68(): + """Upper limit on Psep scaling (PsepB/qAR) + author: P B Lloyd, CCFE, Culham Science Centre + + fpsepbqar: f-value for upper limit on psepbqar, maximum Psep*Bt/qAR limit + psepbqarmax: maximum permitted value of ratio of Psep*Bt/qAR (MWT/m) + p_plasma_separatrix_mw: Power to conducted to the divertor region (MW) + bt: toroidal field on axis (T) (iteration variable 2) + q95: safety factor q at 95% flux surface + aspect: aspect ratio (iteration variable 1) + rmajor: plasma major radius (m) (iteration variable 3) + i_q95_fixed: Switch that allows for fixing q95 only in this constraint. + q95_fixed: fixed safety factor q at 95% flux surface + """ + if fortran.constraint_variables.i_q95_fixed == 1: + cc = ( + ( + ( + fortran.physics_variables.p_plasma_separatrix_mw + * fortran.physics_variables.bt + ) + / ( + fortran.constraint_variables.q95_fixed + * fortran.physics_variables.aspect + * fortran.physics_variables.rmajor + ) + ) + / fortran.constraint_variables.psepbqarmax + - 1.0 * fortran.constraint_variables.fpsepbqar + ) + err = ( + fortran.physics_variables.p_plasma_separatrix_mw + * fortran.physics_variables.bt + ) / ( + fortran.constraint_variables.q95_fixed + * fortran.physics_variables.aspect + * fortran.physics_variables.rmajor + ) - fortran.constraint_variables.psepbqarmax + else: + cc = ( + ( + ( + fortran.physics_variables.p_plasma_separatrix_mw + * fortran.physics_variables.bt + ) + / ( + fortran.physics_variables.q95 + * fortran.physics_variables.aspect + * fortran.physics_variables.rmajor + ) + ) + / fortran.constraint_variables.psepbqarmax + - 1.0 * fortran.constraint_variables.fpsepbqar + ) + err = ( + fortran.physics_variables.p_plasma_separatrix_mw + * fortran.physics_variables.bt + ) / ( + fortran.physics_variables.q95 + * fortran.physics_variables.aspect + * fortran.physics_variables.rmajor + ) - fortran.constraint_variables.psepbqarmax + + return ConstraintResult(cc, fortran.constraint_variables.psepbqarmax, err) + + +@ConstraintManager.register_constraint(72, "Pa", "<=") +def constraint_equation_72(): + """Upper limit on central Solenoid Tresca yield stress + author: P B Lloyd, CCFE, Culham Science Centre + + In the case if the bucked and wedged option ( i_tf_bucking >= 2 ) the constrained + stress is the largest the largest stress of the + - CS stress at maximum current (conservative as the TF inward pressure is not taken + into account) + - CS stress at flux swing (no current in CS) from the TF inward pressure + This allow to cover the 2 worst stress scenario in the bucked and wedged design + Otherwise (free standing TF), the stress limits are only set by the CS stress at max current + Reverse the sign so it works as an inequality constraint (tmp_cc > 0) + This will have no effect if it is used as an equality constraint because it will be squared. + + foh_stress: f-value for Tresca yield criterion in Central Solenoid + alstroh: allowable hoop stress in Central Solenoid structural material (Pa) + s_shear_cs_peak: Maximum shear stress coils/central solenoid (Pa) + sig_tf_cs_bucked: Maximum shear stress in CS case at flux swing (no current in CS) + can be significant for the bucked and weged design + i_tf_bucking: switch for TF structure design + """ + # bucked and wedged desing + if ( + fortran.tfcoil_variables.i_tf_bucking >= 2 + and fortran.build_variables.i_tf_inside_cs == 0 + ): + cc = ( + max( + fortran.pfcoil_variables.s_shear_cs_peak, + fortran.tfcoil_variables.sig_tf_cs_bucked, + ) + / fortran.pfcoil_variables.alstroh + - 1.0 * fortran.constraint_variables.foh_stress + ) + err = fortran.pfcoil_variables.alstroh - max( + fortran.pfcoil_variables.s_shear_cs_peak, + fortran.tfcoil_variables.sig_tf_cs_bucked, + ) + # Free standing CS + else: + cc = ( + fortran.pfcoil_variables.s_shear_cs_peak / fortran.pfcoil_variables.alstroh + - 1.0 * fortran.constraint_variables.foh_stress + ) + err = ( + fortran.pfcoil_variables.alstroh - fortran.pfcoil_variables.s_shear_cs_peak + ) + + return ConstraintResult(cc, fortran.pfcoil_variables.alstroh, err) + + +@ConstraintManager.register_constraint(73, "MW", ">=") +def constraint_equation_73(): + """Lower limit to ensure separatrix power is greater than the L-H power + auxiliary power + Related to constraint 15 + author: P B Lloyd, CCFE, Culham Science Centre + + fplhsep: F-value for Psep >= Plh + Paux : for consistency of two values of separatrix power + p_l_h_threshold_mw: L-H mode power threshold (MW) + p_plasma_separatrix_mw: power to be conducted to the divertor region (MW) + p_hcd_injected_total_mw : inout real : total auxiliary injected power (MW) + """ + cc = ( + 1.0 + - fortran.physics_variables.fplhsep + * fortran.physics_variables.p_plasma_separatrix_mw + / ( + fortran.physics_variables.p_l_h_threshold_mw + + fortran.current_drive_variables.p_hcd_injected_total_mw + ) + ) + return ConstraintResult( + cc, + fortran.physics_variables.p_plasma_separatrix_mw, + fortran.physics_variables.p_plasma_separatrix_mw * cc, + ) + + +@ConstraintManager.register_constraint(74, "K", "<=") +def constraint_equation_74(): + """Upper limit to ensure TF coil quench temperature < tmax_croco + ONLY used for croco HTS coil + author: P B Lloyd, CCFE, Culham Science Centre + + fcqt: f-value: TF coil quench temparature remains below tmax_croco + croco_quench_temperature: CroCo strand: Actual temp reached during a quench (K) + tmax_croco: CroCo strand: maximum permitted temp during a quench (K) + """ + cc = ( + fortran.tfcoil_variables.croco_quench_temperature + / fortran.tfcoil_variables.tmax_croco + - 1.0 * fortran.constraint_variables.fcqt + ) + return ConstraintResult( + cc, + fortran.tfcoil_variables.croco_quench_temperature, + fortran.tfcoil_variables.croco_quench_temperature * cc, + ) + + +@ConstraintManager.register_constraint(75, "A/m2", "<=") +def constraint_equation_75(): + """Upper limit to ensure that TF coil current / copper area < Maximum value + ONLY used for croco HTS coil + author: P B Lloyd, CCFE, Culham Science Centre + + copperA_m2: TF coil current / copper area (A/m2) + copperA_m2_max: Maximum TF coil current / copper area (A/m2) + f_coppera_m2: f-value for TF coil current / copper area < copperA_m2_max + """ + cc = ( + fortran.rebco_variables.coppera_m2 / fortran.rebco_variables.coppera_m2_max + - 1.0 * fortran.rebco_variables.f_coppera_m2 + ) + return ConstraintResult( + cc, fortran.rebco_variables.coppera_m2, fortran.rebco_variables.coppera_m2 * cc + ) + + +@ConstraintManager.register_constraint(76, "m-3", "<=") +def constraint_equation_76(): + """Upper limit for Eich critical separatrix density model: Added for issue 558 + author: P B Lloyd, CCFE, Culham Science Centre + + Eich critical separatrix density model + Added for issue 558 with ref to http://iopscience.iop.org/article/10.1088/1741-4326/aaa340/pdf + + alpha_crit: critical ballooning parameter value + nesep_crit: critical electron density at separatrix [m-3] + kappa: plasma separatrix elongation (calculated if i_plasma_geometry = 1-5, 7 or 9) + triang: plasma separatrix triangularity (calculated if i_plasma_geometry = 1, 3-5 or 7) + aspect: aspect ratio (iteration variable 1) + p_plasma_separatrix_mw: power to conducted to the divertor region (MW) + dlimit(7)array : density limit (/m3) as calculated using various models + fnesep: f-value for Eich critical separatrix density + """ + # TODO: why on earth are these variables being set here!? Should they be local? + fortran.physics_variables.alpha_crit = (fortran.physics_variables.kappa**1.2) * ( + 1.0 + 1.5 * fortran.physics_variables.triang + ) + fortran.physics_variables.nesep_crit = ( + 5.9 + * fortran.physics_variables.alpha_crit + * (fortran.physics_variables.aspect ** (-2.0 / 7.0)) + * (((1.0 + (fortran.physics_variables.kappa**2.0)) / 2.0) ** (-6.0 / 7.0)) + * ((fortran.physics_variables.p_plasma_separatrix_mw * 1.0e6) ** (-11.0 / 70.0)) + * fortran.physics_variables.dlimit[6] + ) + + cc = ( + fortran.physics_variables.nesep / fortran.physics_variables.nesep_crit + - 1.0 * fortran.constraint_variables.fnesep + ) + return ConstraintResult( + cc, fortran.physics_variables.nesep, fortran.physics_variables.nesep * cc + ) + + +@ConstraintManager.register_constraint(77, "A/turn", "<=") +def constraint_equation_77(): + """Equation for maximum TF current per turn upper limit + author: P B Lloyd, CCFE, Culham Science Centre + + fcpttf: f-value for TF coil current per turn + cpttf_max : allowable TF coil current per turn [A/turn] + c_tf_turn : TF coil current per turn [A/turn] + """ + cc = ( + fortran.tfcoil_variables.c_tf_turn / fortran.tfcoil_variables.cpttf_max + - 1.0 * fortran.constraint_variables.fcpttf + ) + return ConstraintResult( + cc, fortran.tfcoil_variables.cpttf_max, fortran.tfcoil_variables.cpttf_max * cc + ) + + +@ConstraintManager.register_constraint(78, "", ">=") +def constraint_equation_78(): + """Equation for Reinke criterion, divertor impurity fraction lower limit + author: P B Lloyd, CCFE, Culham Science Centre + + freinke : input : f-value for Reinke criterion (itv 147) + fzmin : input : minimum impurity fraction from Reinke model + fzactual : input : actual impurity fraction + """ + cc = ( + 1.0 + - fortran.constraint_variables.freinke + * fortran.reinke_variables.fzactual + / fortran.reinke_variables.fzmin + ) + return ConstraintResult( + cc, + fortran.reinke_variables.fzmin * (1.0 - cc), + fortran.reinke_variables.fzmin * cc, + ) + + +@ConstraintManager.register_constraint(79, "A/turn", "<=") +def constraint_equation_79(): + """Equation for maximum CS field + author: P B Lloyd, CCFE, Culham Science Centre + + fb_cs_limit_max: F-value for CS mmax field (cons. 79, itvar 149) + b_cs_limit_max: Central solenoid max field limit [T] + b_cs_peak_pulse_start: maximum field in central solenoid at beginning of pulse (T) + b_cs_peak_flat_top_end: maximum field in central solenoid at end of flat-top (EoF) (T) + (Note: original code has "b_cs_peak_flat_top_end/b_cs_peak_pulse_start | peak CS field [T]".) + """ + cc = ( + max( + fortran.pfcoil_variables.b_cs_peak_flat_top_end, + fortran.pfcoil_variables.b_cs_peak_pulse_start, + ) + / fortran.pfcoil_variables.b_cs_limit_max + - 1.0 * fortran.pfcoil_variables.fb_cs_limit_max + ) + return ConstraintResult( + cc, + fortran.pfcoil_variables.b_cs_limit_max, + max( + fortran.pfcoil_variables.b_cs_peak_flat_top_end, + fortran.pfcoil_variables.b_cs_peak_pulse_start, + ) + * cc, + ) + + +@ConstraintManager.register_constraint(80, "MW", ">=") +def constraint_equation_80(): + """Equation for p_plasma_separatrix_mw lower limit + author: J Morris, Culham Science Centre + args : output structure : residual error; constraint value; residual error in physical units; + output string; units string + Lower limit p_plasma_separatrix_mw + #=# physics + #=#=# fp_plasma_separatrix_min_mw, p_plasma_separatrix_mw + Logic change during pre-factoring: err, symbol, units will be assigned only if present. + fp_plasma_separatrix_min_mw : input : F-value for lower limit on p_plasma_separatrix_mw (cons. 80, itvar 153) + p_plasma_separatrix_min_mw : input : Minimum power crossing separatrix p_plasma_separatrix_mw [MW] + p_plasma_separatrix_mw : input : Power crossing separatrix [MW] + """ + cc = ( + 1.0 + - fortran.physics_variables.fp_plasma_separatrix_min_mw + * fortran.physics_variables.p_plasma_separatrix_mw + / fortran.constraint_variables.p_plasma_separatrix_min_mw + ) + return ConstraintResult( + cc, + fortran.constraint_variables.p_plasma_separatrix_min_mw, + fortran.constraint_variables.p_plasma_separatrix_min_mw * cc, + ) + + +@ConstraintManager.register_constraint(81, "m-3", ">=") +def constraint_equation_81(): + """Lower limit to ensure central density is larger that the pedestal one + author: S Kahn, Culham Science Centre + args : output structure : residual error; constraint value; + residual error in physical units; output string; units string + Lower limit ne0 > neped + !#=# physics + !#=#=# ne0, neped + Logic change during pre-factoring: err, symbol, units will be + assigned only if present. + fne0 : input : F-value for constraint on ne0 > neped + ne0 : input : Central electron density [m-3] + neped : input : Electron density at pedestal [m-3] + """ + cc = ( + 1.0 + - fortran.physics_variables.fne0 + * fortran.physics_variables.ne0 + / fortran.physics_variables.neped + ) + return ConstraintResult( + cc, fortran.physics_variables.fne0, fortran.physics_variables.fne0 * cc + ) + + +@ConstraintManager.register_constraint(82, "m", ">=") +def constraint_equation_82(): + """Equation for toroidal consistency of stellarator build + author: J Lion, IPP Greifswald + + ftoroidalgap: f-value for constraint toroidalgap > dx_tf_inboard_out_toroidal + toroidalgap: minimal gap between two stellarator coils + dx_tf_inboard_out_toroidal: total toroidal width of a tf coil + """ + return ConstraintResult( + 1.0 + - fortran.tfcoil_variables.ftoroidalgap + * fortran.tfcoil_variables.toroidalgap + / fortran.tfcoil_variables.dx_tf_inboard_out_toroidal, + fortran.tfcoil_variables.toroidalgap, + fortran.tfcoil_variables.toroidalgap + - fortran.tfcoil_variables.dx_tf_inboard_out_toroidal + / fortran.tfcoil_variables.ftoroidalgap, + ) + + +@ConstraintManager.register_constraint(83, "m", ">=") +def constraint_equation_83(): + """Equation for radial consistency of stellarator build + author: J Lion, IPP Greifswald + + f_avspace: f-value for constraint available_radial_space > required_radial_space + available_radial_space: avaible space in radial direction as given by each s.-configuration + required_radial_space: required space in radial direction + """ + cc = ( + 1.0 + - fortran.build_variables.f_avspace + * fortran.build_variables.available_radial_space + / fortran.build_variables.required_radial_space + ) + return ConstraintResult( + cc, + fortran.build_variables.available_radial_space * (1.0 - cc), + fortran.build_variables.required_radial_space * cc, + ) + + +@ConstraintManager.register_constraint(84, "", ">=") +def constraint_equation_84(): + """Equation for the lower limit of beta + author: J Lion, IPP Greifswald + + fbeta_min: f-value for constraint beta-beta_fast_alpha > beta_min + beta_min: Lower limit for beta + beta: plasma beta + """ + cc = ( + 1.0 + - fortran.constraint_variables.fbeta_min + * fortran.physics_variables.beta + / fortran.physics_variables.beta_min + ) + return ConstraintResult( + cc, + fortran.physics_variables.beta_min * (1.0 - cc), + fortran.physics_variables.beta * cc, + ) + + +@ConstraintManager.register_constraint(85, "years", "=") +def constraint_equation_85(): + """Equality constraint for the centerpost (CP) lifetime + Author : S Kahn + + Depending on the chosen option i_cp_lifetime: + - 0 : The CP full power year lifelime is set by the user (cplife_input) + - 1 : The CP lifelime is equal to the divertor one + - 2 : The CP lifetime is equal to the breeding blankets one + - 3 : The CP lifetime is equal to the plant one + + cplife: calculated CP full power year lifetime (years) + life_blkt_fpy: calculated first wall/blanket power year lifetime (years) + divlife: calculated divertor power year lifetime (years) + i_cp_lifetime: switch chosing which plant element the CP + the CP lifetime must equate + """ + # The CP lifetime is equal to the the divertor one + if fortran.cost_variables.i_cp_lifetime == 0: + cc = 1.0 - fortran.cost_variables.cplife / fortran.cost_variables.cplife_input + + elif fortran.cost_variables.i_cp_lifetime == 1: + cc = 1.0 - fortran.cost_variables.cplife / fortran.cost_variables.divlife + + # The CP lifetime is equal to the tritium breeding blankets / FW one + elif fortran.cost_variables.i_cp_lifetime == 2: + cc = 1.0 - fortran.cost_variables.cplife / fortran.fwbs_variables.life_blkt_fpy + + elif fortran.cost_variables.i_cp_lifetime == 3: + cc = 1.0 - fortran.cost_variables.cplife / fortran.cost_variables.tlife + + return ConstraintResult( + cc, + fortran.cost_variables.divlife * (1.0 - cc), + fortran.cost_variables.divlife * cc, + ) + + +@ConstraintManager.register_constraint(86, "m", "<=") +def constraint_equation_86(): + """Upper limit on the turn edge length in the TF winding pack + Author : S Kahn + + t_turn_tf: TF coil turn edge length including turn insulation [m] + f_t_turn_tf: f-value for TF turn edge length constraint + t_turn_tf_max: TF turn edge length including turn insulation upper limit [m] + """ + cc = ( + fortran.tfcoil_variables.t_turn_tf / fortran.tfcoil_variables.t_turn_tf_max + - 1.0 * fortran.tfcoil_variables.f_t_turn_tf + ) + return ConstraintResult( + cc, + fortran.tfcoil_variables.t_turn_tf_max * (1.0 - cc), + fortran.tfcoil_variables.t_turn_tf_max * cc, + ) + + +@ConstraintManager.register_constraint(87, "MW", "<=") +def constraint_equation_87(): + """Equation for TF coil cryogenic power upper limit + author: S. Kahn, CCFE, Culham Science Centre + + p_cryo_plant_electric_mw: cryogenic plant power (MW) + f_crypmw: f-value for maximum cryogenic plant power + p_cryo_plant_electric_max_mw: Maximum cryogenic plant power (MW) + """ + cc = ( + fortran.heat_transport_variables.p_cryo_plant_electric_mw + / fortran.heat_transport_variables.p_cryo_plant_electric_max_mw + - 1.0 * fortran.heat_transport_variables.f_crypmw + ) + return ConstraintResult( + cc, + fortran.heat_transport_variables.p_cryo_plant_electric_max_mw * (1.0 - cc), + fortran.heat_transport_variables.p_cryo_plant_electric_mw * cc, + ) + + +@ConstraintManager.register_constraint(88, "", "<=") +def constraint_equation_88(): + """Equation for TF coil vertical strain upper limit (absolute value) + author: CPS Swanson, PPPL, USA + + fstr_wp: f-value for TF coil strain + str_wp_max: Allowable maximum TF coil vertical strain + str_wp: Constrained TF coil vertical strain + """ + return ConstraintResult( + abs(fortran.tfcoil_variables.str_wp) / fortran.tfcoil_variables.str_wp_max + - 1.0 * fortran.constraint_variables.fstr_wp, + fortran.tfcoil_variables.str_wp_max, + fortran.tfcoil_variables.str_wp_max + - abs(fortran.tfcoil_variables.str_wp) / fortran.constraint_variables.fstr_wp, + ) + + +@ConstraintManager.register_constraint(89, "A/m2", "<=") +def constraint_equation_89(): + """Upper limit to ensure that the Central Solenoid [OH] coil current / copper area < Maximum value + author: G Turkington, CCFE, Culham Science Centre + + copperaoh_m2: CS coil current at EOF / copper area [A/m2] + copperaoh_m2_max: maximum coil current / copper area [A/m2] + f_copperaoh_m2: f-value for CS coil current / copper area + """ + cc = ( + fortran.rebco_variables.copperaoh_m2 / fortran.rebco_variables.copperaoh_m2_max + - 1.0 * fortran.rebco_variables.f_copperaoh_m2 + ) + return ConstraintResult( + cc, + fortran.rebco_variables.copperaoh_m2, + fortran.rebco_variables.copperaoh_m2 * cc, + ) + + +@ConstraintManager.register_constraint(90, "", ">=") +def constraint_equation_90(): + """Lower limit for CS coil stress load cycles + author: A. Pearce, G Turkington CCFE, Culham Science Centre + + fncycle: f-value for constraint n_cycle > n_cycle_min + n_cycle: Allowable number of cycles for CS + n_cycle_min: Minimum required cycles for CS + """ + if ( + fortran.cost_variables.ibkt_life == 1 + and fortran.cs_fatigue_variables.bkt_life_csf == 1 + ): + fortran.cs_fatigue_variables.n_cycle_min = fortran.cost_variables.bktcycles + + cc = ( + 1.0 + - fortran.constraint_variables.fncycle + * fortran.cs_fatigue_variables.n_cycle + / fortran.cs_fatigue_variables.n_cycle_min + ) + return ConstraintResult( + cc, + fortran.cs_fatigue_variables.n_cycle_min * (1.0 - cc), + fortran.cs_fatigue_variables.n_cycle * cc, + ) + + +@ConstraintManager.register_constraint(91, "MW", ">=") +def constraint_equation_91(): + """Lower limit to ensure ECRH te is greater than required te for ignition + at lower values for n and B. Or if the design point is ECRH heatable (if i_plasma_ignited==0) + stellarators only (but in principle usable also for tokamaks). + author: J Lion, IPP Greifswald + + fecrh_ignition: f-value for constraint powerht_local > powerscaling + max_gyrotron_frequency: Max. av. gyrotron frequency + te0_ecrh_achievable: Max. achievable electron temperature at ignition point + """ + # Achievable ECRH te needs to be larger than needed te for igntion + if fortran.physics_variables.i_plasma_ignited == 0: + cc = ( + 1.0 + - fortran.constraint_variables.fecrh_ignition + * ( + fortran.stellarator_variables.powerht_constraint + + fortran.current_drive_variables.p_hcd_primary_extra_heat_mw + ) + / fortran.stellarator_variables.powerscaling_constraint + ) + else: + cc = ( + 1.0 + - fortran.constraint_variables.fecrh_ignition + * fortran.stellarator_variables.powerht_constraint + / fortran.stellarator_variables.powerscaling_constraint + ) + + return ConstraintResult( + cc, + fortran.stellarator_variables.powerscaling_constraint * (1.0 - cc), + fortran.stellarator_variables.powerht_constraint * cc, + ) + + +@ConstraintManager.register_constraint(92, "", "=") +def constraint_equation_92(): + """Equation for checking is D/T ratio is consistent, and sums to 1. + author: G Turkington, UKAEA + + f_deuterium: fraction of deuterium ions + f_tritium: fraction of tritium ions + f_helium3: fraction of helium-3 ions + """ + f_deuterium = 1.0 - ( + fortran.physics_variables.f_tritium + fortran.physics_variables.f_helium3 + ) + cc = 1.0 - ( + f_deuterium + + fortran.physics_variables.f_tritium + + fortran.physics_variables.f_helium3 + ) + + return ConstraintResult(cc, 1.0, cc) + + +def constraint_eqns(m: int, ieqn: int): + """Evaluates the constraints given the current state of PROCESS. + + :param m: The number of constraints to evaluate + :param ieqn: Evaluates the 'ieqn'th constraint equation (index starts at 1) + or all equations if <= 0 + """ + + if ieqn > 0: + i1 = ieqn - 1 + i2 = ieqn + else: + i1 = 0 + i2 = m + + cc, con, err, symbol, units = [], [], [], [], [] + + for i in range(i1, i2): + constraint_id = fortran.numerics.icc[i].item() + try: + constraint = ConstraintManager.get_constraint(constraint_id) + + if constraint is None: + tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units = getattr( + fortran.constraints, f"constraint_eqn_{constraint_id:03d}" + )() + else: + result = constraint.constraint_equation() + tmp_cc, tmp_con, tmp_err = ( + result.normalised_residual, + result.constraint_value, + result.constraint_error, + ) + tmp_symbol, tmp_units = constraint.symbol, constraint.units + + except AttributeError as e: + error_msg = f"Constraint equation {i + 1} cannot be found" + raise ProcessError(error_msg) from e + + if np.isnan(tmp_cc) or np.isinf(tmp_cc) or abs(tmp_cc) > 9.99e99: + error_msg = ( + f"Constraint equation {constraint_id} returned an invalid residual" + ) + + raise ProcessValueError(error_msg, cc=tmp_cc) + + # Reverse the sign so it works as an inequality constraint (cc(i) > 0) + # This will have no effect if it is used as an equality constraint because it will be squared. + cc.append(-tmp_cc) + con.append(tmp_con) + err.append(tmp_err) + symbol.append(tmp_symbol) + units.append(tmp_units) + + return np.array(cc), np.array(con), np.array(err), symbol, units + + +def init_constraint_variables(): + """Initialise the constraint variables""" + fortran.constraint_variables.auxmin = 0.1 + fortran.constraint_variables.beta_poloidal_max = 0.19 + fortran.constraint_variables.bigqmin = 10.0 + fortran.constraint_variables.bmxlim = 12.0 + fortran.constraint_variables.fauxmn = 1.0 + fortran.constraint_variables.fbeta_poloidal_eps = 1.0 + fortran.constraint_variables.fbeta_poloidal = 1.0 + fortran.constraint_variables.fbeta_max = 1.0 + fortran.constraint_variables.fbeta_min = 1.0 + fortran.constraint_variables.fcpttf = 1.0 + fortran.constraint_variables.fr_conducting_wall = 1.0 + fortran.constraint_variables.fdene = 1.0 + fortran.constraint_variables.fdtmp = 1.0 + fortran.constraint_variables.fflutf = 1.0 + fortran.constraint_variables.fp_fusion_total_max_mw = 1.0 + fortran.constraint_variables.fgamcd = 1.0 + fortran.constraint_variables.fpflux_div_heat_load_mw = 1.0 + fortran.constraint_variables.fiooic = 0.5 + fortran.constraint_variables.fipir = 1.0 + fortran.constraint_variables.q95_fixed = 3.0 + fortran.constraint_variables.fjohc = 1.0 + fortran.constraint_variables.fjohc0 = 1.0 + fortran.constraint_variables.fjprot = 1.0 + fortran.constraint_variables.fl_h_threshold = 1.0 + fortran.constraint_variables.fmva = 1.0 + fortran.constraint_variables.fnbshinef = 1.0 + fortran.constraint_variables.fncycle = 1.0 + fortran.constraint_variables.fnesep = 1.0 + fortran.constraint_variables.foh_stress = 1.0 + fortran.constraint_variables.fpeakb = 1.0 + fortran.constraint_variables.fpinj = 1.0 + fortran.constraint_variables.fpnetel = 1.0 + fortran.constraint_variables.fportsz = 1.0 + fortran.constraint_variables.fpsepbqar = 1.0 + fortran.constraint_variables.fpsepr = 1.0 + fortran.constraint_variables.fptemp = 1.0 + fortran.constraint_variables.fptfnuc = 1.0 + fortran.constraint_variables.fq = 1.0 + fortran.constraint_variables.fqval = 1.0 + fortran.constraint_variables.fradpwr = 0.99 + fortran.constraint_variables.fradwall = 1.0 + fortran.constraint_variables.freinke = 1.0 + fortran.constraint_variables.frminor = 1.0 + fortran.constraint_variables.fstrcase = 1.0 + fortran.constraint_variables.fstrcond = 1.0 + fortran.constraint_variables.fstr_wp = 1.0 + fortran.constraint_variables.fmaxvvstress = 1.0 + fortran.constraint_variables.ftbr = 1.0 + fortran.constraint_variables.ft_burn = 1.0 + fortran.constraint_variables.ftcycl = 1.0 + fortran.constraint_variables.ftmargoh = 1.0 + fortran.constraint_variables.ftmargtf = 1.0 + fortran.constraint_variables.ft_current_ramp_up = 1.0 + fortran.constraint_variables.ftpeak = 1.0 + fortran.constraint_variables.fvdump = 1.0 + fortran.constraint_variables.fvs = 1.0 + fortran.constraint_variables.fvvhe = 1.0 + fortran.constraint_variables.fwalld = 1.0 + fortran.constraint_variables.fzeffmax = 1.0 + fortran.constraint_variables.gammax = 2.0 + fortran.constraint_variables.i_q95_fixed = 0 + fortran.constraint_variables.pflux_fw_rad_max = 1.0 + fortran.constraint_variables.mvalim = 40.0 + fortran.constraint_variables.f_p_beam_shine_through_max = 1e-3 + fortran.constraint_variables.nflutfmax = 1.0e23 + fortran.constraint_variables.p_plasma_separatrix_min_mw = 150.0 + fortran.constraint_variables.f_fw_rad_max = 3.33 + fortran.constraint_variables.pflux_fw_rad_max_mw = 0.0 + fortran.constraint_variables.pnetelin = 1.0e3 + fortran.constraint_variables.p_fusion_total_max_mw = 1.5e3 + fortran.constraint_variables.psepbqarmax = 9.5 + fortran.constraint_variables.pseprmax = 25.0 + fortran.constraint_variables.ptfnucmax = 1e-3 + fortran.constraint_variables.tbrmin = 1.1 + fortran.constraint_variables.t_burn_min = 1.0 + fortran.constraint_variables.tcycmn = 0.0 + fortran.constraint_variables.t_current_ramp_up_min = 1.0 + fortran.constraint_variables.vvhealw = 1.0 + fortran.constraint_variables.walalw = 1.0 + fortran.constraint_variables.f_alpha_energy_confinement_min = 5.0 + fortran.constraint_variables.falpha_energy_confinement = 1.0 + fortran.constraint_variables.fniterpump = 1.0 + fortran.constraint_variables.zeffmax = 3.6 + fortran.constraint_variables.fpoloidalpower = 1.0 + fortran.constraint_variables.fpsep = 1.0 + fortran.constraint_variables.fcqt = 1.0 + fortran.constraint_variables.fecrh_ignition = 1.0 diff --git a/process/final.py b/process/final.py index 87d3c46819..d0d4e495a2 100644 --- a/process/final.py +++ b/process/final.py @@ -2,13 +2,13 @@ from tabulate import tabulate +import process.constraints as constraints from process import output as op from process import ( process_output as po, ) from process.fortran import ( constants, - constraints, numerics, ) from process.objectives import objective_function @@ -67,20 +67,15 @@ def output_evaluation(): f2py_compatible_to_string(i) for i in numerics.lablcc[numerics.icc[: numerics.neqns + numerics.nineqns] - 1] ] - units = [f2py_compatible_to_string(i) for i in units] - physical_constraint = [ - f"{c} {u}" for c, u in zip(value.tolist(), units, strict=False) - ] - physical_residual = [ - f"{c} {u}" for c, u in zip(residual.tolist(), units, strict=False) - ] + physical_constraint = [f"{c} {u}" for c, u in zip(value, units, strict=False)] + physical_residual = [f"{c} {u}" for c, u in zip(residual, units, strict=False)] table_data = { "Constraint Name": labels, - "Constraint Type": symbols.tolist(), + "Constraint Type": symbols, "Physical constraint": physical_constraint, "Constraint residual": physical_residual, - "Normalised residual": residual_error.tolist(), + "Normalised residual": residual_error, } po.write(constants.nout, tabulate(table_data, headers="keys")) diff --git a/process/init.py b/process/init.py index a51c835a59..8a9cad0072 100644 --- a/process/init.py +++ b/process/init.py @@ -12,6 +12,7 @@ from process.blanket_library import init_blanket_library, init_primary_pumping_variables from process.build import init_build_variables from process.buildings import init_buildings_variables +from process.constraints import ConstraintManager, init_constraint_variables from process.costs import init_cost_variables from process.cs_fatigue import init_cs_fatigue_variables from process.current_drive import init_current_drive_variables @@ -271,7 +272,7 @@ def init_all_module_vars(): init_vacuum_variables() init_pf_power_variables() init_build_variables() - fortran.constraint_variables.init_constraint_variables() + init_constraint_variables() init_pulse_variables() init_rebco_variables() init_reinke_variables() @@ -1134,7 +1135,7 @@ def check_process(inputs): # noqa: ARG001 def set_active_constraints(): """Set constraints provided in the input file as 'active'""" num_constraints = 0 - for i in range(fortran.numerics.ipeqns): + for i in range(ConstraintManager.num_constraints()): if fortran.numerics.icc[i] != 0: fortran.numerics.active_constraints[fortran.numerics.icc[i] - 1] = True num_constraints += 1 diff --git a/process/input.py b/process/input.py index 64a77d215e..7a13934fe2 100644 --- a/process/input.py +++ b/process/input.py @@ -9,6 +9,7 @@ from warnings import warn from process import fortran +from process.constraints import ConstraintManager from process.exceptions import ProcessValidationError, ProcessValueError from process.utilities.f2py_string_patch import ( f2py_compatible_to_string, @@ -104,10 +105,10 @@ def __post_init__(self): "maxcal": InputVariable(fortran.global_variables, int, range=(0, 10000)), "minmax": InputVariable(fortran.numerics, int), "neqns": InputVariable( - fortran.numerics, int, range=(1, fortran.numerics.ipeqns.item()) + fortran.numerics, int, range=(1, ConstraintManager.num_constraints()) ), "nineqns": InputVariable( - fortran.numerics, int, range=(1, fortran.numerics.ipeqns.item()) + fortran.numerics, int, range=(1, ConstraintManager.num_constraints()) ), "alphaj": InputVariable(fortran.physics_variables, float, range=(0.0, 10.0)), "alphan": InputVariable(fortran.physics_variables, float, range=(0.0, 10.0)), diff --git a/process/scan.py b/process/scan.py index 1e634bd44a..5e334ca88c 100644 --- a/process/scan.py +++ b/process/scan.py @@ -3,6 +3,7 @@ import numpy as np from tabulate import tabulate +import process.constraints as constraints import process.process_output as process_output from process.caller import write_output_files from process.exceptions import ProcessValueError @@ -10,7 +11,6 @@ build_variables, constants, constraint_variables, - constraints, cost_variables, cs_fatigue_variables, current_drive_variables, @@ -470,8 +470,8 @@ def post_optimise(self, ifail: int): equality_constraint_table.append([ name, sym[i], - f"{con2[i]} {f2py_compatible_to_string(lab[i])}", - f"{err[i]} {f2py_compatible_to_string(lab[i])}", + f"{con2[i]} {lab[i]}", + f"{err[i]} {lab[i]}", con1[i], ]) process_output.ovarre( @@ -513,8 +513,8 @@ def post_optimise(self, ifail: int): inequality_constraint_table.append([ name, sym[i], - f"{con2[i]} {f2py_compatible_to_string(lab[i])}", - f"{err[i]} {f2py_compatible_to_string(lab[i])}", + f"{con2[i]} {lab[i]}", + f"{err[i]} {lab[i]}", ]) process_output.ovarre( constants.mfile, diff --git a/process/utilities/errorlist.json b/process/utilities/errorlist.json index 0b7290dacd..34f27887fb 100644 --- a/process/utilities/errorlist.json +++ b/process/utilities/errorlist.json @@ -12,8 +12,8 @@ "errors": [ { "no": 1, - "level": 3, - "message": "CONSTRAINTS: Do not use constraint 7 if i_plasma_ignited=1" + "level": 0, + "message": "REMOVED" }, { "no": 2, @@ -27,38 +27,38 @@ }, { "no": 4, - "level": 3, - "message": "CONSTRAINTS: Do not use constraint 28 if i_plasma_ignited=1" + "level": 0, + "message": "REMOVED" }, { "no": 5, - "level": 3, - "message": "CONSTRAINTS: temp_fw_peak = 0 ==> i_pulsed_plant=0; do not use constraint 39 if i_pulsed_plant=0" + "level": 0, + "message": "REMOVED" }, { "no": 6, - "level": 3, - "message": "CONSTRAINTS: tcycmn = 0 ==> i_pulsed_plant=0; do not use constraint 42 if i_pulsed_plant=0" + "level": 0, + "message": "REMOVED" }, { "no": 7, - "level": 3, - "message": "CONSTRAINTS: Do not use constraint 43 if itart=0" + "level": 0, + "message": "REMOVED" }, { "no": 8, - "level": 3, - "message": "CONSTRAINTS: Do not use constraint 44 if itart=0" + "level": 0, + "message": "REMOVED" }, { "no": 9, - "level": 3, - "message": "CONSTRAINTS: Do not use constraint 45 if itart=0" + "level": 0, + "message": "REMOVED" }, { "no": 10, - "level": 3, - "message": "CONSTRAINTS: Do not use constraint 46 if itart=0" + "level": 0, + "message": "REMOVED" }, { "no": 11, @@ -67,18 +67,18 @@ }, { "no": 12, - "level": 3, - "message": "CONSTRAINTS: Do not use constraint 50 if ife=0" + "level": 0, + "message": "REMOVED" }, { "no": 13, - "level": 3, - "message": "CONSTRAINTS: No such constraint equation number..." + "level": 0, + "message": "REMOVED" }, { "no": 14, - "level": 3, - "message": "CONSTRAINTS: NaN/infty error for constraint equation..." + "level": 0, + "message": "REMOVED" }, { "no": 15, @@ -202,12 +202,12 @@ }, { "no": 39, - "level": 1, + "level": 0, "message": "REMOVED" }, { "no": 40, - "level": 2, + "level": 0, "message": "REMOVED" }, { @@ -257,7 +257,7 @@ }, { "no": 50, - "level": 3, + "level": 0, "message": "REMOVED" }, { @@ -292,23 +292,23 @@ }, { "no": 57, - "level": 3, - "message": "CONVXC: Illegal iteration variable number" + "level": 0, + "message": "REMOVED" }, { "no": 58, "level": 3, - "message": "CONVXC: Iteration variable is zero; change its initial value or lower bound" + "message": "REMOVED" }, { "no": 59, "level": 3, - "message": "CONVXC: NaN error for iteration variable" + "message": "REMOVED" }, { "no": 60, - "level": 3, - "message": "CONVXC: scale(i) = 0 for iteration variable" + "level": 0, + "message": "REMOVED" }, { "no": 61, @@ -1187,8 +1187,8 @@ }, { "no": 236, - "level": 3, - "message": "CHECK: The 1/R toroidal B field dependency constraint is being depreciated." + "level": 0, + "message": "REMOVED" }, { "no": 237, @@ -1452,8 +1452,8 @@ }, { "no": 289, - "level": 3, - "message": "CHECK: Constraint equation 22 has been deprecated." + "level": 0, + "message": "REMOVED" } ] } diff --git a/source/fortran/constraint_equations.f90 b/source/fortran/constraint_equations.f90 deleted file mode 100755 index 1df6e6967d..0000000000 --- a/source/fortran/constraint_equations.f90 +++ /dev/null @@ -1,3442 +0,0 @@ -! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - -module constraints - !! author: J Morris - !! - !! Module defining the constraint equations and the routine that evaluates the - !! constraint equations. - - ! Import modules -#ifndef dp - use, intrinsic :: iso_fortran_env, only: dp=>real64 -#endif - use error_handling, only: report_error, idiags, fdiags - - implicit none - - public :: constraint_eqns - -! type constraint_args_type -! real(dp) :: cc -! !! Residual error in constraint equation -! real(dp) :: con -! !! constraint value for constraint equation in physical units -! real(dp) :: err -! !! residual error in constraint equation in physical units -! character(len=1) :: symbol -! !! `=<`, `>`, `<` symbol for constraint equation denoting its type -! character(len=10) :: units -! !! constraint units in constraint equation -! end type - -contains - - subroutine constraint_eqns(m,ieqn,cc,con,err,symbol,units) - !! Routine that formulates the constraint equations - !! - !! **author: P J Knight** (UKAEA) - !! - !! **author: J Morris** (UKAEA) - !! - !! if `ieqn` is zero or negative, evaluate all the constraint equations, otherwise - !! evaluate only the `ieqn`th equation. The code attempts to make \( cc(i) = 0 \) for all - !! \( i \in \{1,\cdots,m\} \) equations. All relevant consistency equations should be active in - !! order to make a self-consistent machine. - !! - !! **References** - !! - !! 1. - use numerics, only: icc - - implicit none - - integer, intent(in) :: m - !! Number of constraint equations to solve - - integer, intent(in) :: ieqn - !! Switch for constraint equations to evaluate; - - real(dp), dimension(m), intent(out) :: cc - !! Residual error in equation i - - real(dp), optional, dimension(m), intent(out) :: con - !! constraint value for equation i in physical units - - real(dp), optional, dimension(m), intent(out) :: err - !! residual error in equation i in physical units - - character(len=1), optional, dimension(m), intent(out) :: symbol - !! `=<`, `>`, `<` symbol for equation i denoting its type - - character*10, optional, dimension(m), intent(out) :: units - !! constraint units in equation i - - ! Local variables - integer :: i, i1, i2 - - real(dp) :: tmp_cc = 0 - !! Residual error in constraint equation - real(dp) :: tmp_con = 0 - !! constraint value for constraint equation in physical units - real(dp) :: tmp_err = 0 - !! residual error in constraint equation in physical units - character(len=1) :: tmp_symbol = ' ' - !! `=<`, `>`, `<` symbol for constraint equation denoting its type - character(len=10) :: tmp_units = ' ' - !! constraint units in constraint equation - - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - ! If ieqn is positive, only evaluate the 'ieqn'th constraint residue, - ! otherwise evaluate all m constraint residues - if (ieqn > 0) then - i1 = ieqn ; i2 = ieqn - else - i1 = 1 ; i2 = m - end if - - ! Consistency (equality) constraints should converge to zero. - do i = i1,i2 - - ! The constraint value in physical units is - ! a) for consistency equations, the quantity to be equated, or - ! b) for limit equations, the limiting value. - ! The symbol is = for a consistency equation, < for an upper limit - ! or > for a lower limit. - select case (icc(i)) - - ! Relationship between beta, temperature (keV) and density - case (1); call constraint_eqn_001(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Global plasma power balance equation - case (2); call constraint_eqn_002(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Global power balance equation for ions - case (3); call constraint_eqn_003(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Global power balance equation for electrons - case (4); call constraint_eqn_004(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for density upper limit - case (5); call constraint_eqn_005(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for epsilon beta-poloidal upper limit - case (6); call constraint_eqn_006(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for hot beam ion density - case (7); call constraint_eqn_007(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for neutron wall load upper limit - case (8); call constraint_eqn_008(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for fusion power upper limit - case (9); call constraint_eqn_009(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Obsolete - case (10); call constraint_eqn_010(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for radial build - case (11); call constraint_eqn_011(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for volt-second capability lower limit - case (12); call constraint_eqn_012(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for burn time lower limit - case (13); call constraint_eqn_013(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation to fix number of NBI decay lengths to plasma centre - case (14); call constraint_eqn_014(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for L-H power threshold limit - case (15); call constraint_eqn_015(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for net electric power lower limit - case (16); call constraint_eqn_016(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for radiation power upper limit - case (17); call constraint_eqn_017(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for divertor heat load upper limit - case (18); call constraint_eqn_018(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for MVA upper limit - case (19); call constraint_eqn_019(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for neutral beam tangency radius upper limit - case (20); call constraint_eqn_020(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for minor radius lower limit - case (21); call constraint_eqn_021(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Obsolete - case (22); call constraint_eqn_022(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for conducting shell radius / rminor upper limit - case (23); call constraint_eqn_023(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for beta upper limit - case (24); call constraint_eqn_024(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for peak toroidal field upper limit - case (25); call constraint_eqn_025(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for Central Solenoid current density upper limit at EOF - case (26); call constraint_eqn_026(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for Central Solenoid current density upper limit at BOP - case (27); call constraint_eqn_027(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for fusion gain (big Q) lower limit - case (28); call constraint_eqn_028(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for inboard major radius - case (29); call constraint_eqn_029(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for injection power upper limit - case (30); call constraint_eqn_030(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for TF coil case stress upper limit (SCTF) - case (31); call constraint_eqn_031(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for TF coil conduit stress upper limit (SCTF) - case (32); call constraint_eqn_032(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for TF coil operating/critical J upper limit (SCTF) - case (33); call constraint_eqn_033(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for TF coil dump voltage upper limit (SCTF) - case (34); call constraint_eqn_034(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for TF coil J_wp/J_prot upper limit (SCTF) - case (35); call constraint_eqn_035(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for TF coil s/c temperature margin lower limit (SCTF) - case (36); call constraint_eqn_036(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for current drive gamma upper limit - case (37); call constraint_eqn_037(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for first wall temperature upper limit - case (39); call constraint_eqn_039(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for auxiliary power lower limit - case (40); call constraint_eqn_040(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for plasma current ramp-up time lower limit - case (41); call constraint_eqn_041(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for cycle time lower limit - case (42); call constraint_eqn_042(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for average centrepost temperature - case (43); call constraint_eqn_043(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for centrepost temperature upper limit (TART) - case (44); call constraint_eqn_044(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for edge safety factor lower limit (TART) - case (45); call constraint_eqn_045(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for Ip/Irod upper limit (TART) - case (46); call constraint_eqn_046(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for TF coil toroidal thickness upper limit - case (47); call constraint_eqn_047(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for poloidal beta upper limit - case (48); call constraint_eqn_048(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Issue #508 Remove RFP option Equation to ensure reversal parameter F is negative - case (49); call constraint_eqn_049(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! IFE option: Equation for repetition rate upper limit - case (50); call constraint_eqn_050(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation to enforce startup flux = available startup flux - case (51); call constraint_eqn_051(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for tritium breeding ratio lower limit - case (52); call constraint_eqn_052(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for fast neutron fluence on TF coil upper limit - case (53); call constraint_eqn_053(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for peak TF coil nuclear heating upper limit - case (54); call constraint_eqn_054(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for helium concentration in vacuum vessel upper limit - case (55); call constraint_eqn_055(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for power through separatrix / major radius upper limit - case (56); call constraint_eqn_056(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Obsolete - case (57); call constraint_eqn_057(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Obsolete - case (58); call constraint_eqn_058(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for neutral beam shine-through fraction upper limit - case (59); call constraint_eqn_059(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for Central Solenoid s/c temperature margin lower limit - case (60); call constraint_eqn_060(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for availability limit - case (61); call constraint_eqn_061(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Lower limit on f_alpha_energy_confinement the ratio of alpha particle to energy confinement times - case (62); call constraint_eqn_062(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Upper limit on niterpump (vacuum_model = simple) - case (63); call constraint_eqn_063(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Upper limit on Zeff - case (64); call constraint_eqn_064(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Limit TF dump time to calculated quench time - case (65); call constraint_eqn_065(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Limit on rate of change of energy in poloidal field - case (66); call constraint_eqn_066(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Simple upper limit on radiation wall load - case (67); call constraint_eqn_067(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! New Psep scaling (PsepB/qAR) - case (68); call constraint_eqn_068(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Ensure separatrix power is less than value from Kallenbach divertor - case (69); call constraint_eqn_069(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Separatrix temperature consistency - case (70); call constraint_eqn_070(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Separatrix density consistency - case (71); call constraint_eqn_071(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Central Solenoid Tresca yield criterion - case (72); call constraint_eqn_072(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! ensure separatrix power is greater than the L-H power + auxiliary power - case (73); call constraint_eqn_073(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! ensure TF coil quench temperature < tmax_croco - case (74); call constraint_eqn_074(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! ensure that TF coil current / copper area < Maximum value ONLY used for croco HTS coil - case (75); call constraint_eqn_075(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Eich critical separatrix density model - case (76); call constraint_eqn_076(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for maximum TF current per turn upper limit - case (77); call constraint_eqn_077(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for Reinke criterion, divertor impurity fraction lower limit - case (78); call constraint_eqn_078(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Equation for maximum CS field - case (79); call constraint_eqn_079(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Lower limit p_plasma_separatrix_mw - case (80); call constraint_eqn_080(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Constraint equation making sure that ne(0) > ne(ped) - case (81); call constraint_eqn_081(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Constraint equation making sure that stellarator coils dont touch in toroidal direction - case (82); call constraint_eqn_082(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Constraint ensuring radial build consistency for stellarators - case (83); call constraint_eqn_083(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Constraint for lower limit of beta - case (84); call constraint_eqn_084(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Constraint for CP lifetime - case (85); call constraint_eqn_085(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Constraint for turn dimension - case (86); call constraint_eqn_086(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Constraint for cryogenic power - case (87); call constraint_eqn_087(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Constraint for TF coil strain - case (88); call constraint_eqn_088(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! ensure that OH coil current / copper area < Maximum value ONLY used for croco HTS coil - case (89); call constraint_eqn_089(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Constraint for minimum CS stress load cycles - case (90); call constraint_eqn_090(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Constraint for indication of ECRH ignitability - case (91); call constraint_eqn_091(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - ! Constraint for D/T ratio - case (92); call constraint_eqn_092(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - case default - - idiags(1) = icc(i) - call report_error(13) - tmp_cc = 0 - tmp_con = 0 - tmp_err = 0 - tmp_symbol = ' ' - tmp_units = ' ' - - end select - - cc(i) = tmp_cc - if (present(con)) then - con(i) = tmp_con - if (present(err)) err(i) = tmp_err - if (present(symbol)) symbol(i) = tmp_symbol - if (present(units)) units(i) = tmp_units - end if - - if (isnan(cc(i)).or.(abs(cc(i))>9.99D99)) then - - ! Add debugging lines as appropriate... - select case (icc(i)) - - ! Relationship between beta, temperature (keV) and density (consistency equation) - case (1); call constraint_err_001() - ! Equation for net electric power lower limit - case (16); call constraint_err_016() - ! Equation for injection power upper limit - case (30); call constraint_err_030() - ! Limit on rate of change of energy in poloidal field - case (66); call constraint_err_066() - - end select - - idiags(1) = icc(i) ; fdiags(1) = cc(i) - call report_error(14) - - end if - - end do - ! Issue 505 Reverse the sign so it works as an inequality constraint (cc(i) > 0) - ! This will have no effect if it is used as an equality constraint because it will be squared. - cc = -cc - - end subroutine constraint_eqns - - !--- Error-handling routines - - subroutine constraint_err_001() - !! Error in: Relationship between beta, temperature (keV) and density (consistency equation) - !! author: P B Lloyd, CCFE, Culham Science Centre - use physics_variables, only: beta_fast_alpha, beta_beam, dene, ten, nd_ions_total, tin, btot, beta - write(*,*) 'beta_fast_alpha = ', beta_fast_alpha - write(*,*) 'beta_beam = ', beta_beam - write(*,*) 'dene = ', dene - write(*,*) 'ten = ', ten - write(*,*) 'nd_ions_total = ', nd_ions_total - write(*,*) 'tin = ', tin - write(*,*) 'btot = ',btot - write(*,*) 'beta = ', beta - end subroutine - - subroutine constraint_err_016() - !! Error in: Equation for net electric power lower limit - !! author: P B Lloyd, CCFE, Culham Science Centre - use constraint_variables, only: fpnetel, pnetelin - use heat_transport_variables, only: p_plant_electric_net_mw - implicit none - write(*,*) 'fpnetel = ', fpnetel - write(*,*) 'p_plant_electric_net_mw = ', p_plant_electric_net_mw - write(*,*) 'pnetelin = ', pnetelin - end subroutine - - subroutine constraint_err_030() - !! Error in: Equation for injection power upper limit - !! author: P B Lloyd, CCFE, Culham Science Centre - use current_drive_variables, only: p_hcd_injected_total_mw, p_hcd_injected_max - use constraint_variables, only: fpinj - implicit none - write(*,*) 'fpinj = ', fpinj - write(*,*) 'p_hcd_injected_max = ', p_hcd_injected_max - write(*,*) 'p_hcd_injected_total_mw = ', p_hcd_injected_total_mw - end subroutine - - subroutine constraint_err_066() - !! Error in: Limit on rate of change of energy in poloidal field - !! author: P B Lloyd, CCFE, Culham Science Centre - use constraint_variables, only: fpoloidalpower - use pf_power_variables, only: maxpoloidalpower, peakpoloidalpower - implicit none - write(*,*) 'fpoloidalpower = ', fpoloidalpower - write(*,*) 'maxpoloidalpower = ', maxpoloidalpower - write(*,*) 'peakpoloidalpower = ', peakpoloidalpower - end subroutine constraint_err_066 - - !--- - - subroutine constraint_eqn_001(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! author: J Morris - !! category: equality constraint - !! - !! Relationship between beta, temperature (keV) and density - !! - !! \begin{equation} - !! c_i = 1 - \frac{1}{\beta}\left( \beta_{ft} + \beta_{NBI} + 2 \times 10^3 \mu_0 e - !! \left( \frac{n_e T_e + n_i T_i}{B_{tot}^2} \right) \right) - !! \end{equation} - !! - !! - \( \beta \) -- total plasma beta - !! - \( \beta_{ft} \) -- fast alpha beta component - !! - \( \beta_{NBI} \) -- neutral beam beta component - !! - \( n_e \) -- electron density [m\(^3\)] - !! - \( n_i \) -- total ion density [m\(^3\)] - !! - \( T_e \) -- density weighted average electron temperature [keV] - !! - \( T_i \) -- density weighted average ion temperature [keV] - !! - \( B_{tot} \) -- total toroidal + poloidal field [T] - - use physics_variables, only: beta_fast_alpha, beta_beam, dene, ten, nd_ions_total, tin, btot, beta - use constants, only: electron_charge,rmu0 - - implicit none - - ! type(constraint_args_type), intent(out) :: args - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - !! constraint derived type - - tmp_cc = 1.0D0 - (beta_fast_alpha + beta_beam + & - 2.0D3*rmu0*electron_charge * (dene*ten + nd_ions_total*tin)/btot**2 )/beta - tmp_con = beta * (1.0D0 - tmp_cc) - tmp_err = beta * tmp_cc - tmp_symbol = '=' - tmp_units = '' - - end subroutine constraint_eqn_001 - - subroutine constraint_eqn_002(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! author: J. Morris - !! category: equality constraint - !! - !! Global plasma power balance equation - !! - !! \begin{equation} c_i = - !! \end{equation} - !! - !! i_rad_loss : input integer : switch for radiation loss term usage in power balance (see User Guide): - !! i_plasma_ignited : input integer : switch for ignition assumption: - !! pden_electron_transport_loss_mw : input real : electron transport power per volume (MW/m3) - !! pden_ion_transport_loss_mw : input real : ion transport power per volume (MW/m3) - !! pden_plasma_rad_mw : input real : total radiation power per volume (MW/m3) - !! pden_plasma_core_rad_mw : input real : total core radiation power per volume (MW/m3) - !! f_alpha_plasma : input real : fraction of alpha power deposited in plasma - !! pden_alpha_total_mw : input real : alpha power per volume (MW/m3) - !! pden_non_alpha_charged_mw : input real : non-alpha charged particle fusion power per volume (MW/m3) - !! pden_plasma_ohmic_mw : input real : ohmic heating power per volume (MW/m3) - !! p_hcd_injected_total_mw : input real : total auxiliary injected power (MW) - !! vol_plasma : input real : plasma volume (m3) - - use physics_variables, only: i_rad_loss, i_plasma_ignited, pden_electron_transport_loss_mw, pden_ion_transport_loss_mw, pden_plasma_rad_mw, & - pden_plasma_core_rad_mw, f_alpha_plasma, pden_alpha_total_mw, pden_non_alpha_charged_mw, & - pden_plasma_ohmic_mw, vol_plasma - use current_drive_variables, only: p_hcd_injected_total_mw - - implicit none - - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - !! constraint derived type - - ! pscaling : Local real : total transport power per volume (MW/m3) - real(dp) :: pscaling - real(dp) :: pnumerator, pdenom - pscaling = pden_electron_transport_loss_mw + pden_ion_transport_loss_mw - ! Total power lost is scaling power plus radiation: - if (i_rad_loss == 0) then - pnumerator = pscaling + pden_plasma_rad_mw - else if (i_rad_loss == 1) then - pnumerator = pscaling + pden_plasma_core_rad_mw - else - pnumerator = pscaling - end if - - ! if plasma not ignited include injected power - if (i_plasma_ignited == 0) then - pdenom = f_alpha_plasma*pden_alpha_total_mw + pden_non_alpha_charged_mw + pden_plasma_ohmic_mw + p_hcd_injected_total_mw/vol_plasma - else - ! if plasma ignited - pdenom = f_alpha_plasma*pden_alpha_total_mw + pden_non_alpha_charged_mw + pden_plasma_ohmic_mw - end if - - tmp_cc = 1.0D0 - pnumerator / pdenom - tmp_con = pdenom * (1.0D0 - tmp_cc) - tmp_err = pdenom * tmp_cc - tmp_symbol = '=' - tmp_units = 'MW/m3' - - end subroutine constraint_eqn_002 - - subroutine constraint_eqn_003(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Global power balance equation for ions - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Global power balance equation for ions - !! This is a consistency equation (NBI) - !! #=# physics - !! #=#=# consistency - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! i_plasma_ignited : input integer : switch for ignition assumption: - !! pden_ion_transport_loss_mw : input real : ion transport power per volume (MW/m3) - !! piepv : input real : ion/electron equilibration power per volume (MW/m3) - !! f_alpha_plasma : input real : fraction of alpha power deposited in plasma - !! f_pden_alpha_ions_mw : input real : alpha power per volume to ions (MW/m3) - !! p_hcd_injected_ions_mw : input real : auxiliary injected power to ions (MW) - !! vol_plasma : input real : plasma volume (m3) - use physics_variables, only: i_plasma_ignited, pden_ion_transport_loss_mw, piepv, f_alpha_plasma, f_pden_alpha_ions_mw, vol_plasma - use current_drive_variables, only: p_hcd_injected_ions_mw - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - ! No assume plasma ignition: - if (i_plasma_ignited == 0) then - tmp_cc = 1.0D0 - (pden_ion_transport_loss_mw + piepv) / (f_alpha_plasma*f_pden_alpha_ions_mw + p_hcd_injected_ions_mw/vol_plasma) - tmp_con = (f_alpha_plasma*f_pden_alpha_ions_mw + p_hcd_injected_ions_mw/vol_plasma) * (1.0D0 - tmp_cc) - tmp_err = (f_alpha_plasma*f_pden_alpha_ions_mw + p_hcd_injected_ions_mw/vol_plasma) * tmp_cc - tmp_symbol = '=' - tmp_units = 'MW/m3' - ! Plasma ignited: - else - tmp_cc = 1.0D0 - (pden_ion_transport_loss_mw+piepv) / (f_alpha_plasma*f_pden_alpha_ions_mw) - tmp_con = (f_alpha_plasma*f_pden_alpha_ions_mw) * (1.0D0 - tmp_cc) - tmp_err = (f_alpha_plasma*f_pden_alpha_ions_mw) * tmp_cc - tmp_symbol = '=' - tmp_units = 'MW/m3' - end if - - end subroutine constraint_eqn_003 - - subroutine constraint_eqn_004(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Global power balance equation for electrons - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Global power balance equation for electrons - !! This is a consistency equation - !! N.B. This constraint is currently NOT RECOMMENDED for use. - !! #=# physics - !! #=#=# consistency - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! i_rad_loss : input integer : switch for radiation loss term usage in power balance (see User Guide): - !! i_plasma_ignited : input integer : switch for ignition assumption: - !! pden_electron_transport_loss_mw : input real : electron transport power per volume (MW/m3) - !! pden_plasma_rad_mw : input real : total radiation power per volume (MW/m3) - !! pden_plasma_core_rad_mw : input real : total core radiation power per volume (MW/m3) - !! f_alpha_plasma : input real : fraction of alpha power deposited in plasma - !! f_pden_alpha_electron_mw : input real : alpha power per volume to electrons (MW/m3) - !! piepv : input real : ion/electron equilibration power per volume (MW/m3) - !! p_hcd_injected_electrons_mw : input real : auxiliary injected power to electrons (MW) - !! vol_plasma : input real : plasma volume (m3) - use physics_variables, only: i_rad_loss, i_plasma_ignited, pden_electron_transport_loss_mw, pden_plasma_core_rad_mw, f_alpha_plasma, & - f_pden_alpha_electron_mw, piepv, vol_plasma, pden_plasma_rad_mw - use current_drive_variables, only: p_hcd_injected_electrons_mw - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - ! pscaling : Local real : total transport power per volume (MW/m3) - real(dp) :: pscaling - real(dp) :: pnumerator, pdenom - pscaling = pden_electron_transport_loss_mw - ! Total power lost is scaling power plus radiation: - if (i_rad_loss == 0) then - pnumerator = pscaling + pden_plasma_rad_mw - else if (i_rad_loss == 1) then - pnumerator = pscaling + pden_plasma_core_rad_mw - else - pnumerator = pscaling - end if - - ! if plasma not ignited include injected power - if (i_plasma_ignited == 0) then - pdenom = f_alpha_plasma*f_pden_alpha_electron_mw + piepv + p_hcd_injected_electrons_mw/vol_plasma - else - ! if plasma ignited - pdenom = f_alpha_plasma*f_pden_alpha_electron_mw + piepv - end if - - tmp_cc = 1.0D0 - pnumerator / pdenom - tmp_con = pdenom * (1.0D0 - tmp_cc) - tmp_err = pdenom * tmp_cc - tmp_symbol = '=' - tmp_units = 'MW/m3' - - end subroutine constraint_eqn_004 - - subroutine constraint_eqn_005(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for density upper limit - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for density upper limit - !! #=# physics - !! #=#=# fdene, dnelimt - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! i_density_limit : input integer : switch for density limit to enforce (constraint equation 5): - !! fdene : input real : f-value for density limit - !! dene : input real : electron density (/m3) - !! dnelimt : input real : density limit (/m3) - !! dnla : input real : line averaged electron density (m-3) - use physics_variables, only: i_density_limit, dnelimt, dnla, dene - use constraint_variables, only: fdene - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - ! Apply Greenwald limit to line-averaged density - if (i_density_limit == 7) then - tmp_cc = dnla/dnelimt - 1.0D0 * fdene - tmp_con = fdene * dnelimt - tmp_err = fdene * dnelimt - dnla - tmp_symbol = '<' - tmp_units = '/m3' - else - tmp_cc = dene/dnelimt - 1.0D0 * fdene - tmp_con = dnelimt * (1.0D0 - tmp_cc) - tmp_err = dene * tmp_cc - tmp_symbol = '<' - tmp_units = '/m3' - end if - - end subroutine constraint_eqn_005 - - subroutine constraint_eqn_006(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for epsilon beta-poloidal upper limit - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for epsilon beta-poloidal upper limit - !! #=# physics - !! #=#=# fbeta_poloidal_eps, beta_poloidal_eps_max - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! fbeta_poloidal_eps : input real : f-value for epsilon beta-poloidal - !! beta_poloidal_eps_max : input real : maximum (eps*beta_poloidal) - !! eps : input real : inverse aspect ratio - !! beta_poloidal : input real : poloidal beta - use physics_variables, only: beta_poloidal_eps_max, eps, beta_poloidal - use constraint_variables, only: fbeta_poloidal_eps - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = (eps*beta_poloidal)/beta_poloidal_eps_max - 1.0D0 * fbeta_poloidal_eps - tmp_con = beta_poloidal_eps_max * (1.0D0 - tmp_cc) - tmp_err = (eps*beta_poloidal) * tmp_cc - tmp_symbol = '<' - tmp_units = '' - - end subroutine constraint_eqn_006 - - subroutine constraint_eqn_007(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for hot beam ion density - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for hot beam ion density - !! This is a consistency equation (NBI) - !! #=# physics - !! #=#=# consistency - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! i_plasma_ignited : input integer : switch for ignition assumption: - !! Obviously, i_plasma_ignited must be zero if current drive is required. - !! If i_plasma_ignited=1, any auxiliary power is assumed to be used only - !! during plasma start-up, and is excluded from all steady-state - !! power balance calculations. - !! beam_density_out : input real : hot beam ion density from calculation (/m3) - !! nd_beam_ions : input real : hot beam ion density, variable (/m3) - use physics_variables, only: i_plasma_ignited, beam_density_out, nd_beam_ions - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - ! Do not assume plasma ignition: - if (i_plasma_ignited == 0) then - tmp_cc = 1.0D0 - beam_density_out/nd_beam_ions - tmp_con = nd_beam_ions * (1.0D0 - tmp_cc) - tmp_err = nd_beam_ions * tmp_cc - tmp_symbol = '=' - tmp_units = '/m3' - else - tmp_cc = 0 - tmp_con = 0 - tmp_err = 0 - tmp_symbol = '' - tmp_units = '' - call report_error(1) - end if - - end subroutine constraint_eqn_007 - - subroutine constraint_eqn_008(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for neutron wall load upper limit - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for neutron wall load upper limit - !! #=# physics - !! #=#=# fwalld, walalw - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! fwalld : input real : f-value for maximum wall load - !! walalw : input real : allowable wall-load (MW/m2) - !! pflux_fw_neutron_mw : input real : average neutron wall load (MW/m2) - use constraint_variables, only: fwalld, walalw - use physics_variables, only: pflux_fw_neutron_mw - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = pflux_fw_neutron_mw/walalw - 1.0D0 * fwalld - tmp_con = fwalld * walalw - tmp_err = fwalld * walalw - pflux_fw_neutron_mw - tmp_symbol = '<' - tmp_units = 'MW/m2' - - end subroutine constraint_eqn_008 - - subroutine constraint_eqn_009(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for fusion power upper limit - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for fusion power upper limit - !! #=# physics - !! #=#=# fp_fusion_total_max_mw, p_fusion_total_max_mw - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! fp_fusion_total_max_mw : input real : f-value for maximum fusion power - !! p_fusion_total_max_mw : input real : maximum fusion power (MW) - !! p_fusion_total_mw : input real : fusion power (MW) - use constraint_variables, only: fp_fusion_total_max_mw, p_fusion_total_max_mw - use physics_variables, only: p_fusion_total_mw - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = p_fusion_total_mw/p_fusion_total_max_mw - 1.0D0 * fp_fusion_total_max_mw - tmp_con = p_fusion_total_max_mw * (1.0D0 - tmp_cc) - tmp_err = p_fusion_total_mw * tmp_cc - tmp_symbol = '<' - tmp_units = 'MW' - - end subroutine constraint_eqn_009 - - subroutine constraint_eqn_010(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! author: P B Lloyd, CCFE, Culham Science Centre - !! Equation for field at TF coil - !! This is a consistency equation - !! (do not use for stellarators) - !! #=# tfcoil - !! #=#=# consistency - !! rmajor | plasma major radius (m) - !! bt | toroidal field on axis (T) - !! r_b_tf_inboard_peak | radius of maximum toroidal field (m) - !! b_tf_inboard_peak | peak field at toroidal field coil (T) - - !! This constraint is depreciated - - implicit none - - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - !! Constraints output - - ! This constraint is depreciated - call report_error(236) - - tmp_con = 1.0D0 - tmp_err = 0.0D0 - tmp_symbol = '=' - tmp_units = '' - - end subroutine constraint_eqn_010 - - subroutine constraint_eqn_011(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for radial build - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for radial build - !! (This is a consistency equation.) - !! #=# build - !! #=#=# consistency - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! rbld : input real : sum of thicknesses to the major radius (m) - !! rmajor : input real : plasma major radius (m) - use build_variables, only: rbld - use physics_variables, only: rmajor - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = 1.0D0 - rbld/rmajor - tmp_con = rmajor * (1.0D0 - tmp_cc) - tmp_err = rmajor * tmp_cc - tmp_symbol = '=' - tmp_units = 'm' - - end subroutine constraint_eqn_011 - - subroutine constraint_eqn_012(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for volt-second capability lower limit - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for volt-second capability lower limit - !! #=# pfcoil - !! #=#=# fvs, vs_plasma_total_required - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! vs_plasma_total_required : input real : total V-s needed (Wb) - !! vs_plasma_total_required (lower limit) is positive; vs_cs_pf_total_pulse (available) is negative - !! fvs : input real : f-value for flux-swing (V-s) requirement (STEADY STATE) - !! vs_cs_pf_total_pulse : input real : total flux swing for pulse (Wb) - use physics_variables, only: vs_plasma_total_required - use constraint_variables, only: fvs - use pfcoil_variables, only: vs_cs_pf_total_pulse - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - !! vs_cs_pf_total_pulse is negative, requires sign change - tmp_cc = 1.0D0 - fvs * (-vs_cs_pf_total_pulse)/vs_plasma_total_required - tmp_con = vs_plasma_total_required * (1.0D0 - tmp_cc) - tmp_err = vs_plasma_total_required * tmp_cc - tmp_symbol = '>' - tmp_units = 'V.sec' - - end subroutine constraint_eqn_012 - - subroutine constraint_eqn_013(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for burn time lower limit - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for burn time lower limit - !! #=# times - !! #=#=# ft_burn, t_burn_min - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! ft_burn : input real : f-value for minimum burn time - !! t_burn : input real : burn time (s) (calculated if i_pulsed_plant=1) - !! t_burn_min : input real : minimum burn time (s) - use constraint_variables, only: ft_burn,t_burn_min - use times_variables, only: t_burn - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = 1.0D0 - ft_burn * t_burn/t_burn_min - tmp_con = t_burn_min / ft_burn - tmp_err = t_burn_min / ft_burn - t_burn - tmp_symbol = '>' - tmp_units = 'sec' - - end subroutine constraint_eqn_013 - - subroutine constraint_eqn_014(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation to fix number of NBI decay lengths to plasma centre - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation to fix number of NBI decay lengths to plasma centre - !! This is a consistency equation - !! #=# current_drive - !! #=#=# consistency - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! n_beam_decay_lengths_core : input real : neutral beam e-decay lengths to plasma centre - !! tbeamin : input real : permitted neutral beam e-decay lengths to plasma centre - use current_drive_variables, only: n_beam_decay_lengths_core, tbeamin - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = 1.0D0 - n_beam_decay_lengths_core/tbeamin - tmp_con = tbeamin * (1.0D0 - tmp_cc) - tmp_err = tbeamin * tmp_cc - tmp_symbol = '=' - tmp_units = '' - - end subroutine constraint_eqn_014 - - subroutine constraint_eqn_015(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for L-H power threshold limit - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for L-H power threshold limit - !! #=# physics - !! #=#=# fl_h_threshold, p_l_h_threshold_mw - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! fl_h_threshold : input real : f-value for L-H power threshold - !! p_l_h_threshold_mw : input real : L-H mode power threshold (MW) - !! p_plasma_separatrix_mw : input real : power to conducted to the divertor region (MW) - use constraint_variables, only: fl_h_threshold - use physics_variables, only: p_l_h_threshold_mw, p_plasma_separatrix_mw - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = 1.0D0 - fl_h_threshold * p_plasma_separatrix_mw / p_l_h_threshold_mw - tmp_con = p_l_h_threshold_mw - tmp_err = p_l_h_threshold_mw - p_plasma_separatrix_mw / fl_h_threshold - if (fl_h_threshold > 1.0D0) then - tmp_symbol = '>' - else - tmp_symbol = '<' - end if - tmp_units = 'MW' - - end subroutine constraint_eqn_015 - - subroutine constraint_eqn_016(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for net electric power lower limit - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for net electric power lower limit - !! #=# heat_transport - !! #=#=# fpnetel, pnetelin - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! fpnetel : input real : f-value for net electric power - !! p_plant_electric_net_mw : input real : net electric power (MW) - !! pnetelin : input real : required net electric power (MW) - use constraint_variables, only: fpnetel, pnetelin - use heat_transport_variables, only: p_plant_electric_net_mw - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = 1.0D0 - fpnetel * p_plant_electric_net_mw / pnetelin - tmp_con = pnetelin - tmp_err = p_plant_electric_net_mw - pnetelin / fpnetel - tmp_symbol = '>' - tmp_units = 'MW' - - end subroutine constraint_eqn_016 - - subroutine constraint_eqn_017(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for radiation power upper limit - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for radiation power upper limit - !! #=# physics - !! #=#=# fradpwr, pradmaxpv - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! f_alpha_plasma : input real : fraction of alpha power deposited in plasma - !! p_hcd_injected_total_mw : input real : total auxiliary injected power (MW) - !! vol_plasma : input real : plasma volume (m3) - !! pden_alpha_total_mw : input real : alpha power per volume (MW/m3) - !! pden_non_alpha_charged_mw : input real : non-alpha charged particle fusion power per volume (MW/m3) - !! pden_plasma_ohmic_mw : input real : ohmic heating power per volume (MW/m3) - !! fradpwr : input real : f-value for core radiation power limit - !! pden_plasma_rad_mw : input real : total radiation power per volume (MW/m3) - use physics_variables, only: f_alpha_plasma, vol_plasma, pden_alpha_total_mw, pden_non_alpha_charged_mw, pden_plasma_ohmic_mw, pden_plasma_rad_mw - use current_drive_variables, only: p_hcd_injected_total_mw - use constraint_variables, only: fradpwr - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - real(dp) :: pradmaxpv - !! Maximum possible power/vol_plasma that can be radiated (local) - - pradmaxpv = p_hcd_injected_total_mw/vol_plasma + pden_alpha_total_mw*f_alpha_plasma + pden_non_alpha_charged_mw + pden_plasma_ohmic_mw - tmp_cc = pden_plasma_rad_mw/pradmaxpv - 1.0D0 * fradpwr - tmp_con = pradmaxpv * (1.0D0 - tmp_cc) - tmp_err = pden_plasma_rad_mw * tmp_cc - tmp_symbol = '<' - tmp_units = 'MW/m3' - - end subroutine constraint_eqn_017 - - subroutine constraint_eqn_018(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for divertor heat load upper limit - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for divertor heat load upper limit - !! #=# divertor - !! #=#=# fpflux_div_heat_load_mw, pflux_div_heat_load_max_mw - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! fpflux_div_heat_load_mw : input real : f-value for divertor heat load - !! pflux_div_heat_load_max_mw : input real : heat load limit (MW/m2) - !! pflux_div_heat_load_mw : input real : divertor heat load (MW/m2) - use constraint_variables, only: fpflux_div_heat_load_mw - use divertor_variables, only: pflux_div_heat_load_max_mw, pflux_div_heat_load_mw - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = pflux_div_heat_load_mw/pflux_div_heat_load_max_mw - 1.0D0 * fpflux_div_heat_load_mw - tmp_con = pflux_div_heat_load_max_mw * (1.0D0 - tmp_cc) - tmp_err = pflux_div_heat_load_mw * tmp_cc - tmp_symbol = '<' - tmp_units = 'MW/m2' - - end subroutine constraint_eqn_018 - - subroutine constraint_eqn_019(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for MVA (power) upper limit: resistive TF coil set - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for MVA upper limit - !! #=# tfcoil - !! #=#=# fmva, mvalim - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! p_cp_resistive_mw : input real : peak resistive TF coil inboard leg power (total) (MW) - !! p_tf_leg_resistive_mw : input real : TF coil outboard leg resistive power (total) (MW) - !! fmva : input real : f-value for maximum MVA - !! mvalim : input real : MVA limit for resistive TF coil set (total) (MW) - use tfcoil_variables, only: p_cp_resistive_mw, p_tf_leg_resistive_mw - use constraint_variables, only: fmva, mvalim - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - ! totmva : local real : total MVA in TF coil (MW) - real(dp) :: totmva - - totmva = p_cp_resistive_mw + p_tf_leg_resistive_mw - tmp_cc = totmva/mvalim - 1.0D0 * fmva - tmp_con = mvalim * (1.0D0 - tmp_cc) - tmp_err = totmva * tmp_cc - tmp_symbol = '<' - tmp_units = 'MVA' - - end subroutine constraint_eqn_019 - - subroutine constraint_eqn_020(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for neutral beam tangency radius upper limit - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for neutral beam tangency radius upper limit - !! #=# current_drive - !! #=#=# fportsz, rtanmax - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! fportsz : input real : f-value for neutral beam tangency radius limit - !! rtanmax : input real : maximum tangency radius for centreline of beam (m) - !! rtanbeam : input real : neutral beam centreline tangency radius (m) - use constraint_variables, only: fportsz - use current_drive_variables, only: rtanmax, rtanbeam - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = rtanbeam/rtanmax - 1.0D0 * fportsz - tmp_con = rtanmax * (1.0D0 - tmp_cc) - tmp_err = rtanbeam * tmp_cc - tmp_symbol = '<' - tmp_units = 'm' - - end subroutine constraint_eqn_020 - - subroutine constraint_eqn_021(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for minor radius lower limit - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for minor radius lower limit - !! #=# physics - !! #=#=# frminor, aplasmin - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! frminor : input real : f-value for minor radius limit - !! rminor : input real : plasma minor radius (m) - !! aplasmin : input real : minimum minor radius (m) - use constraint_variables, only: frminor - use physics_variables, only: rminor - use build_variables, only: aplasmin - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = 1.0D0 - frminor * rminor/aplasmin - tmp_con = aplasmin * (1.0D0 - tmp_cc) - tmp_err = aplasmin * tmp_cc - tmp_symbol = '>' - tmp_units = '' - - end subroutine constraint_eqn_021 - - subroutine constraint_eqn_022(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for divertor collision/connection length ratio upper limit - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for divertor collision/connection length ratio upper limit - !! #=# divertor - !! #=#=# fdivcol, rlenmax - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! fdivcol : input real : f-value for divertor collisionality - !! rlenmax : input real : maximum value for length ratio (rlclolcn) - !! rlclolcn : input real : ratio of collision length / connection length - - implicit none - - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - !! Constraints output - - ! This constraint is depreciated - call report_error(289) - - tmp_con = 1.0D0 - tmp_err = 0.0D0 - tmp_symbol = '=' - tmp_units = '' - - end subroutine constraint_eqn_022 - - subroutine constraint_eqn_023(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for conducting shell radius / rminor upper limit - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for conducting shell radius / rminor upper limit - !! #=# physics - !! #=#=#fr_conducting_wall, f_r_conducting_wall - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! rminor : input real : plasma minor radius (m) - !! dr_fw_plasma_gap_outboard : input real : gap between plasma and first wall, outboard side (m) - !! dr_fw_outboard : input real : outboard first wall thickness, initial estimate (m) - !! dr_blkt_outboard : input real : outboard blanket thickness (m) - !!fr_conducting_wall : input real : f-value for conducting wall radius / rminor limit - !! f_r_conducting_wall : input real : maximum ratio of conducting wall distance to plasma minor radius for vertical stability - use physics_variables, only: rminor, f_r_conducting_wall - use build_variables, only: dr_fw_plasma_gap_outboard, dr_fw_outboard, dr_blkt_outboard - use constraint_variables, only:fr_conducting_wall - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - ! rcw : local real : conducting shell radius (m) - real(dp) :: rcw - - rcw = rminor + dr_fw_plasma_gap_outboard + dr_fw_outboard + dr_blkt_outboard - tmp_cc = rcw / (f_r_conducting_wall * rminor) - 1.0D0 * fr_conducting_wall - tmp_con = f_r_conducting_wall*rminor * (1.0D0 - tmp_cc) - tmp_err = rcw * tmp_cc - tmp_symbol = '<' - tmp_units = 'm' - - end subroutine constraint_eqn_023 - - subroutine constraint_eqn_024(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for beta upper limit - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for beta upper limit - !! #=# physics - !! #=#=# fbeta_max, beta_max - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! i_beta_component : input integer : switch for beta limit scaling (constraint equation 24): - !! istell : input integer : switch for stellarator option (set via device.dat): - !! fbeta_max : input real : f-value for beta limit - !! beta_max : input real : allowable beta - !! beta : input real : total plasma beta (calculated if ipedestal =3) - !! beta_fast_alpha : input real : fast alpha beta component - !! beta_beam : input real : neutral beam beta component - !! bt : input real : toroidal field - !! btot : input real : total field - use physics_variables, only: i_beta_component, beta_max, beta, beta_beam, beta_fast_alpha, bt, btot - use stellarator_variables, only: istell - use constraint_variables, only: fbeta_max - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - ! Include all beta components: relevant for both tokamaks and stellarators - if ((i_beta_component == 0).or.(istell /= 0)) then - tmp_cc = beta/beta_max - 1.0D0 * fbeta_max - tmp_con = beta_max - tmp_err = beta_max - beta / fbeta_max - tmp_symbol = '<' - tmp_units = '' - ! Here, the beta limit applies to only the thermal component, not the fast alpha or neutral beam parts - else if (i_beta_component == 1) then - tmp_cc = (beta-beta_fast_alpha-beta_beam)/beta_max - 1.0D0 * fbeta_max - tmp_con = beta_max - tmp_err = beta_max - (beta-beta_fast_alpha-beta_beam) / fbeta_max - tmp_symbol = '<' - tmp_units = '' - ! Beta limit applies to thermal + neutral beam: components of the total beta, i.e. excludes alphas - else if (i_beta_component == 2) then - tmp_cc = (beta-beta_fast_alpha)/beta_max - 1.0D0 * fbeta_max - tmp_con = beta_max * (1.0D0 - tmp_cc) - tmp_err = (beta-beta_fast_alpha) * tmp_cc - tmp_symbol = '<' - tmp_units = '' - ! Beta limit applies to toroidal beta - else if (i_beta_component == 3) then - tmp_cc = (beta*(btot/bt)**2)/beta_max - 1.0D0 * fbeta_max - tmp_con = beta_max - tmp_err = beta_max - (beta*(btot/bt)**2) / fbeta_max - tmp_symbol = '<' - tmp_units = '' - end if - - end subroutine constraint_eqn_024 - - subroutine constraint_eqn_025(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for peak toroidal field upper limit - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for peak toroidal field upper limit - !! #=# tfcoil - !! #=#=# fpeakb, bmxlim - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! fpeakb : input real : f-value for maximum toroidal field - !! bmxlim : input real : maximum peak toroidal field (T) - !! b_tf_inboard_peak : input real : mean peak field at TF coil (T) - use constraint_variables, only: fpeakb, bmxlim - use tfcoil_variables, only: b_tf_inboard_peak - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = b_tf_inboard_peak/bmxlim - 1.0D0 * fpeakb - tmp_con = bmxlim * (1.0D0 - tmp_cc) - tmp_err = b_tf_inboard_peak * tmp_cc - tmp_symbol = '<' - tmp_units = 'T' - - end subroutine constraint_eqn_025 - - subroutine constraint_eqn_026(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for Central Solenoid current density upper limit at EOF - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for Central Solenoid current density upper limit at EOF - !! #=# pfcoil - !! #=#=# fjohc, j_cs_critical_flat_top_end - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! fjohc : input real : f-value for central solenoid current at end-of-flattop - !! j_cs_critical_flat_top_end : input real : allowable central solenoid current density at end of flat-top (A/m2) - !! j_cs_flat_top_end : input real : central solenoid overall current density at end of flat-top (A/m2) - use constraint_variables, only: fjohc - use pfcoil_variables, only: j_cs_critical_flat_top_end, j_cs_flat_top_end - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = j_cs_flat_top_end/j_cs_critical_flat_top_end - 1.0D0 * fjohc - tmp_con = j_cs_critical_flat_top_end - tmp_err = j_cs_critical_flat_top_end - j_cs_flat_top_end / fjohc - tmp_symbol = '<' - tmp_units = 'A/m2' - - end subroutine constraint_eqn_026 - - subroutine constraint_eqn_027(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for Central Solenoid current density upper limit at BOP - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for Central Solenoid current density upper limit at BOP - !! #=# pfcoil - !! #=#=# fjohc0, j_cs_critical_pulse_start - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! fjohc0 : input real : f-value for central solenoid current at beginning of pulse - !! j_cs_critical_pulse_start : input real : allowable central solenoid current density at beginning of pulse (A/m2) - !! j_cs_pulse_start : input real : central solenoid overall current density at beginning of pulse (A/m2) - use constraint_variables, only: fjohc0 - use pfcoil_variables, only: j_cs_critical_pulse_start, j_cs_pulse_start - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = j_cs_pulse_start/j_cs_critical_pulse_start - 1.0D0 * fjohc0 - tmp_con = j_cs_critical_pulse_start - tmp_err = j_cs_critical_pulse_start - j_cs_pulse_start / fjohc0 - tmp_symbol = '<' - tmp_units = 'A/m2' - - end subroutine constraint_eqn_027 - - subroutine constraint_eqn_028(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for fusion gain (big Q) lower limit - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for fusion gain (big Q) lower limit - !! #=# physics - !! #=#=# fqval, bigqmin - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! fqval : input real : pf-value for Q - !! bigq : input real : Fusion gain; P_fusion / (P_injection + P_ohmic) - !! bigqmin : input real : minimum fusion gain Q - !! i_plasma_ignited : input integer : switch for ignition assumption: - !! Obviously, i_plasma_ignited must be zero if current drive is required. - !! If i_plasma_ignited=1, any auxiliary power is assumed to be used only - !! during plasma start-up, and is excluded from all steady-state - !! power balance calculations. - use constraint_variables, only: fqval, bigqmin - use current_drive_variables, only: bigq - use physics_variables, only: i_plasma_ignited - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - ! if plasma is not ignited ... - if (i_plasma_ignited == 0) then - tmp_cc = 1.0D0 - fqval * bigq/bigqmin - tmp_con = bigqmin * (1.0D0 - tmp_cc) - tmp_err = bigqmin * tmp_cc - tmp_symbol = '>' - tmp_units = '' - ! if plasma is ignited report error - else - tmp_cc = 0 - tmp_con = 0 - tmp_err = 0 - tmp_symbol = '' - tmp_units = '' - call report_error(4) - end if - - end subroutine constraint_eqn_028 - - subroutine constraint_eqn_029(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for inboard major radius: This is a consistency equation - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for inboard major radius: This is a consistency equation - !! #=# build - !! #=#=# consistency - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! rmajor : input real : plasma major radius (m) (iteration variable 3) - !! rminor : input real : plasma minor radius (m) - !! rinboard : input real : plasma inboard radius (m) - use physics_variables, only: rmajor, rminor - use build_variables, only: rinboard - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = 1.0D0 - (rmajor - rminor) / rinboard - tmp_con = rinboard * (1.0D0 - tmp_cc) - tmp_err = rinboard * tmp_cc - tmp_symbol = '=' - tmp_units = 'm' - - end subroutine constraint_eqn_029 - - subroutine constraint_eqn_030(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for injection power upper limit - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for injection power upper limit - !! #=# current_drive - !! #=#=# fpinj, p_hcd_injected_max - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! p_hcd_injected_total_mw : input real : total auxiliary injected power (MW) - !! fpinj : input real : f-value for injection power - !! p_hcd_injected_max : input real : Maximum allowable value for injected power (MW) - use current_drive_variables, only: p_hcd_injected_total_mw, p_hcd_injected_max - use constraint_variables, only: fpinj - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = p_hcd_injected_total_mw/p_hcd_injected_max - 1.0D0 * fpinj - tmp_con = p_hcd_injected_max - tmp_err = p_hcd_injected_max - p_hcd_injected_total_mw / fpinj - tmp_symbol = '<' - tmp_units = 'MW' - - end subroutine constraint_eqn_030 - - subroutine constraint_eqn_031(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for TF coil case stress upper limit (SCTF) - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for TF coil case stress upper limit (SCTF) - !! #=# tfcoil - !! #=#=# fstrcase, sig_tf_case_max - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! fstrcase : input real : f-value for TF coil case stress - !! sig_tf_case_max : input real : Allowable maximum shear stress in TF coil case (Tresca criterion) (Pa) - !! sig_tf_case : input real : Constrained stress in TF coil case (Pa) - use constraint_variables, only: fstrcase - use tfcoil_variables, only: sig_tf_case_max, sig_tf_case - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = sig_tf_case/sig_tf_case_max - 1.0D0 * fstrcase - tmp_con = sig_tf_case_max - tmp_err = sig_tf_case_max - sig_tf_case / fstrcase - tmp_symbol = '<' - tmp_units = 'Pa' - - end subroutine constraint_eqn_031 - - subroutine constraint_eqn_032(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for TF coil conduit stress upper limit (SCTF) - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for TF coil conduit stress upper limit (SCTF) - !! #=# tfcoil - !! #=#=# fstrcond, sig_tf_wp_max - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! fstrcond : input real : f-value for TF coil conduit stress - !! sig_tf_wp_max : input real : Allowable maximum shear stress in TF coil conduit (Tresca criterion) (Pa) - !! sig_tf_wp : input real : Constrained stress in TF conductor conduit (Pa) - use constraint_variables, only: fstrcond - use tfcoil_variables, only: sig_tf_wp_max, sig_tf_wp - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = sig_tf_wp/sig_tf_wp_max - 1.0D0 * fstrcond - tmp_con = sig_tf_wp_max - tmp_err = sig_tf_wp_max - sig_tf_wp / fstrcond - tmp_symbol = '<' - tmp_units = 'Pa' - - end subroutine constraint_eqn_032 - - subroutine constraint_eqn_033(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for TF coil operating/critical J upper limit (SCTF) - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for TF coil operating/critical J upper limit (SCTF) - !! #=# tfcoil - !! #=#=# fiooic, jwdgcrt - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! fiooic : input real : f-value for TF coil operating current / critical - !! jwdgcrt : input real : critical current density for winding pack (A/m2) - !! j_tf_wp : input real : winding pack current density (A/m2) - use constraint_variables, only: fiooic - use tfcoil_variables, only: jwdgcrt, j_tf_wp - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - if (fiooic > 0.7D0) call report_error(285) - tmp_cc = j_tf_wp/jwdgcrt - 1.0D0 * fiooic - tmp_con = jwdgcrt * (1.0D0 - tmp_cc) - tmp_err = j_tf_wp * tmp_cc - tmp_symbol = '<' - tmp_units = 'A/m2' - - end subroutine constraint_eqn_033 - - subroutine constraint_eqn_034(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for TF coil dump voltage upper limit (SCTF) - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for TF coil dump voltage upper limit (SCTF) - !! #=# tfcoil - !! #=#=# fvdump, vdalw - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! fvdump : input real : f-value for dump voltage - !! vdalw : input real : max voltage across TF coil during quench (kV) - !! vtfskv : input real : voltage across a TF coil during quench (kV) - use constraint_variables, only: fvdump - use tfcoil_variables, only: vdalw, vtfskv - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = vtfskv/vdalw - 1.0D0 * fvdump - tmp_con = vdalw - tmp_err = vdalw - vtfskv - tmp_symbol = '<' - tmp_units = 'V' - - end subroutine constraint_eqn_034 - - subroutine constraint_eqn_035(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for TF coil J_wp/J_prot upper limit (SCTF) - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for TF coil J_wp/J_prot upper limit (SCTF) - !! #=# tfcoil - !! #=#=# fjprot, jwdgpro - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! fjprot : input real : f-value for TF coil winding pack current density - !! jwdgpro : input real : allowable TF coil winding pack current density, for dump temperature rise protection (A/m2) - !! j_tf_wp : input real : winding pack current density (A/m2) - use constraint_variables, only: fjprot - use tfcoil_variables, only: jwdgpro, j_tf_wp - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = j_tf_wp/jwdgpro - 1.0D0 * fjprot - tmp_con = jwdgpro - tmp_err = j_tf_wp - jwdgpro - tmp_symbol = '<' - tmp_units = 'A/m2' - - end subroutine constraint_eqn_035 - - subroutine constraint_eqn_036(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for TF coil s/c temperature margin lower limit (SCTF) - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for TF coil s/c temperature margin lower limit (SCTF) - !! #=# tfcoil - !! #=#=# ftmargtf, tmargmin_tf - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! ftmargtf : input real : f-value for TF coil temperature margin - !! tmargtf : input real : TF coil temperature margin (K) - !! tmargmin_tf : input real : minimum allowable temperature margin : TF coils (K) - use constraint_variables, only: ftmargtf - use tfcoil_variables, only: tmargtf, tmargmin_tf - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = 1.0D0 - ftmargtf * tmargtf/tmargmin_tf - tmp_con = tmargmin_tf - tmp_err = tmargmin_tf - tmargtf - tmp_symbol = '>' - tmp_units = 'K' - - end subroutine constraint_eqn_036 - - subroutine constraint_eqn_037(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for current drive gamma upper limit - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for current drive gamma upper limit - !! #=# current_drive - !! #=#=# fgamcd, gammax - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! fgamcd : input real : f-value for current drive gamma - !! gammax : input real : maximum current drive gamma - !! eta_cd_norm_hcd_primary : input real : normalised current drive efficiency (1.0e20 A/W-m2) - use constraint_variables, only: fgamcd, gammax - use current_drive_variables, only: eta_cd_norm_hcd_primary - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = eta_cd_norm_hcd_primary/gammax - 1.0D0 * fgamcd - tmp_con = gammax * (1.0D0 - tmp_cc) - tmp_err = eta_cd_norm_hcd_primary * tmp_cc - tmp_symbol = '<' - tmp_units = '1E20 A/Wm2' - - end subroutine constraint_eqn_037 - - subroutine constraint_eqn_038(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! TODO Remove - !! Obsolete - !! author: P B Lloyd, CCFE, Culham Science Centre - !! Obsolete - !! #=# empty - !! #=#=# empty - implicit none - ! Dummy formal arguments, for compliance with interface - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = 0 - tmp_con = 0 - tmp_err = 0 - tmp_symbol = '' - tmp_units = '' - - end subroutine constraint_eqn_038 - - subroutine constraint_eqn_039(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for first wall temperature upper limit - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for first wall temperature upper limit - !! #=# fwbs - !! #=#=# ftpeak, temp_fw_max - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! ftpeak : input real : f-value for first wall peak temperature - !! temp_fw_max : input real : maximum temperature of first wall material (K) (i_thermal_electric_conversion>1) - !! temp_fw_peak : input real : peak first wall temperature (K) - use constraint_variables, only: ftpeak - use fwbs_variables, only: temp_fw_max, temp_fw_peak - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - ! If the temperature peak == 0 then report an error - if (temp_fw_peak < 1.0D0) call report_error(5) - tmp_cc = temp_fw_peak/temp_fw_max - 1.0D0 * ftpeak - tmp_con = temp_fw_max * (1.0D0 - tmp_cc) - tmp_err = temp_fw_peak * tmp_cc - tmp_symbol = '<' - tmp_units = 'K' - - end subroutine constraint_eqn_039 - - subroutine constraint_eqn_040(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for auxiliary power lower limit - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for auxiliary power lower limit - !! #=# current_drive - !! #=#=# fauxmn, auxmin - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! fauxmn : input real : f-value for minimum auxiliary power - !! p_hcd_injected_total_mw : input real : total auxiliary injected power (MW) - !! auxmin : input real : minimum auxiliary power (MW) - use constraint_variables, only: fauxmn, auxmin - use current_drive_variables, only: p_hcd_injected_total_mw - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = 1.0D0 - fauxmn * p_hcd_injected_total_mw/auxmin - tmp_con = auxmin * (1.0D0 - tmp_cc) - tmp_err = auxmin * tmp_cc - tmp_symbol = '>' - tmp_units = 'MW' - - end subroutine constraint_eqn_040 - - subroutine constraint_eqn_041(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for plasma current ramp-up time lower limit - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for plasma current ramp-up time lower limit - !! #=# times - !! #=#=# ft_current_ramp_up, t_current_ramp_up_min - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! ft_current_ramp_up : input real : f-value for plasma current ramp-up time - !! t_current_ramp_up : input real : plasma current ramp-up time for current initiation (s) - !! t_current_ramp_up_min : input real : minimum plasma current ramp-up time (s) - use constraint_variables, only: ft_current_ramp_up, t_current_ramp_up_min - use times_variables, only: t_current_ramp_up - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = 1.0D0 - ft_current_ramp_up * t_current_ramp_up/t_current_ramp_up_min - tmp_con = t_current_ramp_up_min * (1.0D0 - tmp_cc) - tmp_err = t_current_ramp_up_min * tmp_cc - tmp_symbol = '>' - tmp_units = 'sec' - - end subroutine constraint_eqn_041 - - subroutine constraint_eqn_042(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for cycle time lower limit - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for cycle time lower limit - !! #=# times - !! #=#=# ftcycl, tcycmn - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! ftcycl : input real : f-value for cycle time - !! t_cycle : input real : full cycle time (s) - !! tcycmn : input real : minimum cycle time (s) - use constraint_variables, only: ftcycl, tcycmn - use times_variables, only: t_cycle - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - ! if the minimum cycle time == 0 report an error - if (tcycmn < 1.0D0) call report_error(6) - tmp_cc = 1.0D0 - ftcycl * t_cycle/tcycmn - tmp_con = tcycmn * (1.0D0 - tmp_cc) - tmp_err = tcycmn * tmp_cc - tmp_symbol = '>' - tmp_units = 'sec' - - end subroutine constraint_eqn_042 - - subroutine constraint_eqn_043(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for average centrepost temperature: This is a consistency equation (TART) - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for average centrepost temperature: This is a consistency equation (TART) - !! #=# tfcoil - !! #=#=# consistency - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! temp_cp_average : input real : average temp of TF coil inboard leg conductor (C)e - !! tcpav2 : input real : centrepost average temperature (C) (for consistency) - !! itart : input integer : switch for spherical tokamak (ST) models: - use tfcoil_variables, only: temp_cp_average, tcpav2 - use physics_variables, only: itart - use tfcoil_variables, only: i_tf_sup - - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - ! if the machine isn't a ST then report error - if (itart == 0) call report_error(7) - - ! For some reasons these lines are needed to make VMCON CONVERGE .... - if ( i_tf_sup == 0 ) then ! Copper case - temp_cp_average = temp_cp_average - 273.15D0 - tcpav2 = tcpav2 - 273.15D0 - end if - - tmp_cc = 1.0D0 - temp_cp_average/tcpav2 - tmp_con = tcpav2 * (1.0D0 - tmp_cc) - tmp_err = tcpav2 * tmp_cc - tmp_symbol = '=' - tmp_units = 'deg C' - - ! For some reasons these lines are needed to make VMCON CONVERGE .... - if ( i_tf_sup == 0 ) then ! Copper case - temp_cp_average = temp_cp_average + 273.15D0 - tcpav2 = tcpav2 + 273.15D0 - end if - - - - end subroutine constraint_eqn_043 - - subroutine constraint_eqn_044(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for centrepost temperature upper limit (TART) - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for centrepost temperature upper limit (TART) - !! #=# tfcoil - !! #=#=# fptemp, ptempalw - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! fptemp : input real : f-value for peak centrepost temperature - !! ptempalw : input real : maximum peak centrepost temperature (K) - !! tcpmax : input real : peak centrepost temperature (K) - !! itart : input integer : switch for spherical tokamak (ST) models: - use constraint_variables, only: fptemp - use tfcoil_variables, only: ptempalw, tcpmax - use physics_variables, only: itart - use tfcoil_variables, only: i_tf_sup - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - ! if the machine isn't a ST then report error - if (itart == 0) call report_error(8) - - ! For some reasons these lines are needed to make VMCON CONVERGE .... - if ( i_tf_sup == 0 ) then ! Copper case - ptempalw = ptempalw - 273.15D0 - tcpmax = tcpmax - 273.15D0 - end if - - tmp_cc = tcpmax/ptempalw - 1.0D0 * fptemp - tmp_con = ptempalw * (1.0D0 - tmp_cc) - tmp_err = tcpmax * tmp_cc - tmp_symbol = '<' - tmp_units = 'deg C' - - ! For some reasons these lines are needed to make VMCON CONVERGE .... - if ( i_tf_sup == 0 ) then ! Copper case - ptempalw = ptempalw + 273.15D0 - tcpmax = tcpmax + 273.15D0 - end if - - end subroutine constraint_eqn_044 - - subroutine constraint_eqn_045(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for edge safety factor lower limit (TART) - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for edge safety factor lower limit (TART) - !! #=# tfcoil - !! #=#=# fq, q95_min - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! fq : input real : f-value for edge safety factor - !! q95 : safety factor 'near' plasma edge - !! (unless i_plasma_current = 2 (ST current scaling), in which case q = mean edge safety factor qbar) - !! q95_min : input real : lower limit for edge safety factor - !! itart : input integer : switch for spherical tokamak (ST) models: - use constraint_variables, only: fq - use physics_variables, only: q95, q95_min, itart - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - ! if the machine isn't a ST then report error - if (itart == 0) call report_error(9) - tmp_cc = 1.0D0 - fq * q95/q95_min - tmp_con = q95_min * (1.0D0 - tmp_cc) - tmp_err = q95_min * tmp_cc - tmp_symbol = '<' - tmp_units = '' - - end subroutine constraint_eqn_045 - - subroutine constraint_eqn_046(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for Ip/Irod upper limit (TART) - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for Ip/Irod upper limit (TART) - !! #=# tfcoil - !! #=#=# fipir, cratmx - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! eps : input real : inverse aspect ratio - !! fipir : input real : f-value for Ip/Irod upper limit - !! c_tf_total : input real : total (summed) current in TF coils (A) - !! plasma_current : input real : plasma current (A) - !! itart : input integer : switch for spherical tokamak (ST) models: - use physics_variables, only: eps, plasma_current, itart - use constraint_variables, only: fipir - use tfcoil_variables, only: c_tf_total - implicit none - ! cratmx : local real : maximum ratio of plasma current to centrepost current - real(dp) :: cratmx - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - ! if the machine isn't a ST then report error - if (itart == 0) call report_error(10) - cratmx = 1.0D0 + 4.91D0*(eps-0.62D0) - tmp_cc = (plasma_current / c_tf_total) / cratmx - 1.0D0 * fipir - tmp_con = cratmx * (1.0D0 - tmp_cc) - tmp_err = plasma_current/c_tf_total * tmp_cc - tmp_symbol = '<' - tmp_units = '' - - end subroutine constraint_eqn_046 - - subroutine constraint_eqn_047(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! TODO Remove - !! Issue #508 Remove RFP option: Relevant only to reversed field pinch devices - !! author: P B Lloyd, CCFE, Culham Science Centre - !! Issue #508 Remove RFP option: Relevant only to reversed field pinch devices - !! Equation for TF coil toroidal thickness upper limit - !! #=# empty - !! #=#=# empty - implicit none - ! Dummy formal arguments, just to comply with the subroutine interface - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = 0 - tmp_con = 0 - tmp_err = 0 - tmp_symbol = '' - tmp_units = '' - - end subroutine constraint_eqn_047 - - subroutine constraint_eqn_048(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for poloidal beta upper limit - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for poloidal beta upper limit - !! #=# physics - !! #=#=# fbeta_poloidal, beta_poloidal_max - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! fbeta_poloidal : input real : rf-value for poloidal beta - !! beta_poloidal_max : input real : maximum poloidal beta - !! beta_poloidal : input real : poloidal beta - use constraint_variables, only: fbeta_poloidal, beta_poloidal_max - use physics_variables, only: beta_poloidal - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = beta_poloidal/beta_poloidal_max - 1.0D0 * fbeta_poloidal - tmp_con = beta_poloidal_max * (1.0D0 - tmp_cc) - tmp_err = beta_poloidal * tmp_cc - tmp_symbol = '<' - tmp_units = '' - - end subroutine constraint_eqn_048 - - subroutine constraint_eqn_049(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! TODO Remove - !! Issue #508 Remove IFE option: Equation for repetition rate upper limit - !! author: P B Lloyd, CCFE, Culham Science Centre - !! Issue #508 Remove IFE option: Equation for repetition rate upper limit - !! #=# empty - !! #=#=# empty - ! Dummy formal arguments, just to comply with the subroutine interface - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = 0 - tmp_con = 0 - tmp_err = 0 - tmp_symbol = '' - tmp_units = '' - - end subroutine constraint_eqn_049 - - subroutine constraint_eqn_050(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! IFE option: Equation for repetition rate upper limit - !! author: P B Lloyd, CCFE, Culham Science Centre - !! author: S I Muldrew, CCFE, Culham Science Centre - !! IFE option: Equation for repetition rate upper limit - !! #=# IFE - !! #=#=# frrmax, rrmax - use ife_variables, only: frrmax, ife, rrmax, reprat - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - if (ife /= 1) then - call report_error(12) - end if - - tmp_cc = reprat/rrmax - 1.0D0 * frrmax - tmp_con = rrmax * (1.0D0 - tmp_cc) - tmp_err = reprat * tmp_cc - tmp_symbol = '<' - tmp_units = 'Hz' - - end subroutine constraint_eqn_050 - - subroutine constraint_eqn_051(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation to enforce startup flux = available startup flux - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation to enforce startup flux = available startup flux - !! #=# pfcoil - !! #=#=# consistency - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! vs_plasma_res_ramp : input real : resistive losses in startup V-s (Wb) - !! vs_plasma_ind_ramp : input real : internal and external plasma inductance V-s (Wb)) - !! vs_cs_pf_total_ramp : input real : total flux swing for startup (Wb) - use physics_variables, only: vs_plasma_res_ramp, vs_plasma_ind_ramp - use pfcoil_variables, only: vs_cs_pf_total_ramp, fvs_cs_pf_total_ramp - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = 1.0D0 - fvs_cs_pf_total_ramp * abs((vs_plasma_res_ramp+vs_plasma_ind_ramp) / vs_cs_pf_total_ramp) - tmp_con = vs_cs_pf_total_ramp * (1.0D0 - tmp_cc) - tmp_err = vs_cs_pf_total_ramp * tmp_cc - tmp_symbol = '=' - tmp_units = 'V.s' - - end subroutine constraint_eqn_051 - - subroutine constraint_eqn_052(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for tritium breeding ratio lower limit - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for tritium breeding ratio lower limit - !! #=# fwbs - !! #=#=# ftbr, tbrmin - !! ? TODO should this only be for certain blanket models ? - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! ftbr : input real : f-value for minimum tritium breeding ratio - !! tbr : input real : tritium breeding ratio (i_blanket_type=2,3 (KIT HCPB/HCLL)) - !! tbrmin : input real : minimum tritium breeding ratio (If i_blanket_type=1, tbrmin=minimum 5-year time-averaged tritium breeding ratio) - use constraint_variables, only: ftbr, tbrmin - use fwbs_variables, only: tbr - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = 1.0D0 - ftbr * tbr/tbrmin - tmp_con = tbrmin * (1.0D0 - tmp_cc) - tmp_err = tbrmin * tmp_cc - tmp_symbol = '>' - tmp_units = '' - - end subroutine constraint_eqn_052 - - subroutine constraint_eqn_053(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for fast neutron fluence on TF coil upper limit - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for fast neutron fluence on TF coil upper limit - !! #=# fwbs - !! #=#=# fflutf, nflutfmax - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! fflutf : input real : f-value for maximum TF coil nuclear heating - !! nflutfmax : input real : max fast neutron fluence on TF coil (n/m2) - !! nflutf : input real : peak fast neutron fluence on TF coil superconductor (n/m2) - use constraint_variables, only: fflutf, nflutfmax - use fwbs_variables, only: nflutf - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = nflutf/nflutfmax - 1.0D0 * fflutf - tmp_con = nflutfmax * (1.0D0 - tmp_cc) - tmp_err = nflutf * tmp_cc - tmp_symbol = '<' - tmp_units = 'neutron/m2' - - end subroutine constraint_eqn_053 - - subroutine constraint_eqn_054(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for peak TF coil nuclear heating upper limit - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for peak TF coil nuclear heating upper limit - !! #=# fwbs - !! #=#=# fptfnuc, ptfnucmax - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! fptfnuc : input real : f-value for maximum TF coil nuclear heating - !! ptfnucmax : input real : maximum nuclear heating in TF coil (MW/m3) - !! ptfnucpm3 : input real : nuclear heating in the TF coil (MW/m3) (blktmodel>0) - use constraint_variables, only: fptfnuc, ptfnucmax - use fwbs_variables, only: ptfnucpm3 - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = ptfnucpm3/ptfnucmax - 1.0D0 * fptfnuc - tmp_con = ptfnucmax * (1.0D0 - tmp_cc) - tmp_err = ptfnucpm3 * tmp_cc - tmp_symbol = '<' - tmp_units = 'MW/m3' - - end subroutine constraint_eqn_054 - - subroutine constraint_eqn_055(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! TODO Remove - !! vvhemax is no longer calculated in PROCESS and this constraint is disabled - implicit none - - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - call report_error(173) - end subroutine constraint_eqn_055 - - subroutine constraint_eqn_056(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for power through separatrix / major radius upper limit - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for power through separatrix / major radius upper limit - !! #=# current_drive - !! #=#=# fnbshinef, f_p_beam_shine_through_max - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! fpsepr : input real : f-value for maximum Psep/R limit - !! pseprmax : input real : maximum ratio of power crossing the separatrix to plasma major radius (Psep/R) (MW/m) - !! p_plasma_separatrix_mw : input real : power to be conducted to the divertor region (MW) - !! rmajor : input real : plasma major radius (m) - use constraint_variables, only: fpsepr, pseprmax - use physics_variables, only: p_plasma_separatrix_mw, rmajor - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = (p_plasma_separatrix_mw/rmajor)/pseprmax - 1.0D0 * fpsepr - tmp_con = pseprmax * (1.0D0 - tmp_cc) - tmp_err = (p_plasma_separatrix_mw/rmajor) * tmp_cc - tmp_symbol = '<' - tmp_units = 'MW/m' - - end subroutine constraint_eqn_056 - - subroutine constraint_eqn_057(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! TODO Remove - !! Obsolete - !! author: P B Lloyd, CCFE, Culham Science Centre - !! Obsolete - !! #=# empty - !! #=#=# empty - ! Dummy formal arguments, just to comply with the subroutine interface - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = 0 - tmp_con = 0 - tmp_err = 0 - tmp_symbol = '' - tmp_units = '' - - end subroutine constraint_eqn_057 - - subroutine constraint_eqn_058(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! TODO Remove - !! Obsolete - !! author: P B Lloyd, CCFE, Culham Science Centre - !! Obsolete - !! #=# empty - !! #=#=# empty - ! Dummy formal arguments, just to comply with the subroutine interface - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = 0 - tmp_con = 0 - tmp_err = 0 - tmp_symbol = '' - tmp_units = '' - - end subroutine constraint_eqn_058 - - subroutine constraint_eqn_059(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for neutral beam shine-through fraction upper limit - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for neutral beam shine-through fraction upper limit - !! #=# current_drive - !! #=#=# fnbshinef, f_p_beam_shine_through_max - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! fnbshinef : input real : f-value for maximum neutral beam shine-through fraction - !! f_p_beam_shine_through_max : input real : maximum neutral beam shine-through fraction - !! f_p_beam_shine_through : input real : neutral beam shine-through fraction - use constraint_variables, only: fnbshinef, f_p_beam_shine_through_max - use current_drive_variables, only: f_p_beam_shine_through - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - tmp_cc = f_p_beam_shine_through/f_p_beam_shine_through_max - 1.0D0 * fnbshinef - tmp_con = f_p_beam_shine_through_max * (1.0D0 - tmp_cc) - tmp_err = f_p_beam_shine_through * tmp_cc - tmp_symbol = '<' - tmp_units = '' - end subroutine constraint_eqn_059 - - subroutine constraint_eqn_060(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for Central Solenoid s/c temperature margin lower limit - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for Central Solenoid s/c temperature margin lower limit - !! #=# tfcoil - !! #=#=# ftmargoh, tmargmin_cs - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! ftmargoh : input real : f-value for central solenoid temperature margin - !! temp_cs_margin : input real : Central solenoid temperature margin (K) - !! tmargmin_cs : input real : Minimum allowable temperature margin : CS (K) - use constraint_variables, only: ftmargoh - use pfcoil_variables, only: temp_cs_margin - use tfcoil_variables, only: tmargmin_cs - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = 1.0D0 - ftmargoh * temp_cs_margin/tmargmin_cs - tmp_con = tmargmin_cs - tmp_err = tmargmin_cs - temp_cs_margin - tmp_symbol = '>' - tmp_units = 'K' - - end subroutine constraint_eqn_060 - - subroutine constraint_eqn_061(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for availability lower limit - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for availability limit - !! #=# cost - !! #=#=# favail, avail_min - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! favail : input real : F-value for minimum availability - !! cfactr : input real : Total plant availability fraction - !! avail_min : input real : Minimum availability - use cost_variables, only: favail, cfactr, avail_min - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = 1.0D0 - favail * cfactr / avail_min - tmp_con = avail_min * (1.0D0 - tmp_cc) - tmp_err = cfactr * tmp_cc - tmp_symbol = '>' - tmp_units = '' - - end subroutine constraint_eqn_061 - - subroutine constraint_eqn_062(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Lower limit on f_alpha_energy_confinement the ratio of alpha particle to energy confinement times - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Lower limit on f_alpha_energy_confinement the ratio of alpha particle to energy confinement times - !! #=# physics - !! #=#=# falpha_energy_confinement, f_alpha_energy_confinement_min - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! falpha_energy_confinement : input real : f-value for lower limit on f_alpha_energy_confinement the ratio of alpha particle to energy confinement - !! t_alpha_confinement : input real : alpha particle confinement time (s) - !! t_energy_confinement : input real : global thermal energy confinement time (sec) - !! f_alpha_energy_confinement_min : input real : Lower limit on f_alpha_energy_confinement the ratio of alpha particle to energy confinement times - use constraint_variables, only: falpha_energy_confinement, f_alpha_energy_confinement_min - use physics_variables, only: t_alpha_confinement, t_energy_confinement - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = 1.0D0 - falpha_energy_confinement * (t_alpha_confinement / t_energy_confinement) / f_alpha_energy_confinement_min - tmp_con = f_alpha_energy_confinement_min - tmp_err = (t_alpha_confinement / t_energy_confinement) * tmp_cc - tmp_symbol = '>' - tmp_units = '' - - end subroutine constraint_eqn_062 - - subroutine constraint_eqn_063(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Upper limit on niterpump (vacuum_model = simple) - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Upper limit on niterpump (vacuum_model = simple) - !! #=# vacuum - !! #=#=# fniterpump, tfno - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! fniterpump : input real : f-value for constraint that number of pumps < tfno - !! tfno : input real : number of TF coils (default = 50 for stellarators) - !! niterpump : input real : number of high vacuum pumps (real number), each with the throughput - use constraint_variables, only: fniterpump - use tfcoil_variables, only: n_tf_coils - use vacuum_variables, only: niterpump - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = niterpump/n_tf_coils - 1.0D0 * fniterpump - tmp_con = n_tf_coils - tmp_err = n_tf_coils * tmp_cc - tmp_symbol = '<' - tmp_units = '' - - end subroutine constraint_eqn_063 - - subroutine constraint_eqn_064(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Upper limit on Zeff - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Upper limit on Zeff - !! #=# physics - !! #=#=# fzeffmax, zeffmax - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! fzeffmax : input real : f-value for maximum zeff - !! zeffmax : input real : maximum value for Zeff - !! zeff : input real : plasma effective charge - use constraint_variables, only: fzeffmax, zeffmax - use physics_variables, only: zeff - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = zeff/zeffmax - 1.0D0 * fzeffmax - tmp_con = zeffmax - tmp_err = zeffmax * tmp_cc - tmp_symbol = '<' - tmp_units = '' - - end subroutine constraint_eqn_064 - - subroutine constraint_eqn_065(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Upper limit on stress of the vacuum vessel that occurs when the TF coil quenches. - !! author: Timothy Nunn, UKAEA - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! fmaxvvstress : input real : f-value for constraint on maximum VV stress - !! max_vv_stress : input real : Maximum permitted stress of the VV (Pa) - !! vv_stress_quench : input real : Stress of the VV (Pa) - use constraint_variables, only: fmaxvvstress - use tfcoil_variables, only: max_vv_stress - use sctfcoil_module, only: vv_stress_quench - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = vv_stress_quench / max_vv_stress - 1.0d0 * fmaxvvstress - tmp_con = max_vv_stress - tmp_err = max_vv_stress * tmp_cc - tmp_symbol = '<' - tmp_units = 'Pa' - - end subroutine constraint_eqn_065 - - subroutine constraint_eqn_066(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Upper limit on rate of change of energy in poloidal field - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Limit on rate of change of energy in poloidal field - !! #=# pfcoil - !! #=#=# fpoloidalpower, maxpoloidalpower - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! fpoloidalpower : input real : f-value for constraint on rate of change of energy in poloidal field - !! maxpoloidalpower : input real : Maximum permitted absolute rate of change of stored energy in poloidal field (MW) - !! peakpoloidalpower : input real : Peak absolute rate of change of stored energy in poloidal field (MW) (11/01/16) - use constraint_variables, only: fpoloidalpower - use pf_power_variables, only: maxpoloidalpower, peakpoloidalpower - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = peakpoloidalpower / maxpoloidalpower - 1.0d0 * fpoloidalpower - tmp_con = maxpoloidalpower - tmp_err = maxpoloidalpower * tmp_cc - tmp_symbol = '<' - tmp_units = 'MW' - - end subroutine constraint_eqn_066 - - subroutine constraint_eqn_067(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Simple upper limit on radiation wall load - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Simple upper limit on radiation wall load - !! #=# physics - !! #=#=# fradwall, pflux_fw_rad_max - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! fradwall : input real : f-value for upper limit on radiation wall load - !! pflux_fw_rad_max : input real : Maximum permitted radiation wall load (MW/m^2) - !! pflux_fw_rad_max_mw : input real : Peak radiation wall load (MW/m^2) - use constraint_variables, only: fradwall, pflux_fw_rad_max, pflux_fw_rad_max_mw - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = pflux_fw_rad_max_mw / pflux_fw_rad_max - 1.0d0 * fradwall - tmp_con = pflux_fw_rad_max - tmp_err = pflux_fw_rad_max * tmp_cc - tmp_symbol = '<' - tmp_units = 'MW/m^2' - - end subroutine constraint_eqn_067 - - subroutine constraint_eqn_068(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Upper limit on Psep scaling (PsepB/qAR) - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! New Psep scaling (PsepB/qAR) - !! Issue #464 - !! #=# physics - !! #=#=# fpsepbqar, psepbqarmax - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! fpsepbqar : input real : f-value for upper limit on psepbqar, maximum Psep*Bt/qAR limit - !! psepbqarmax : input real : maximum permitted value of ratio of Psep*Bt/qAR (MWT/m) - !! p_plasma_separatrix_mw : input real : Power to conducted to the divertor region (MW) - !! bt : input real : toroidal field on axis (T) (iteration variable 2) - !! q95 : input real : safety factor q at 95% flux surface - !! aspect : input real : aspect ratio (iteration variable 1) - !! rmajor : input real : plasma major radius (m) (iteration variable 3) - !! i_q95_fixed : input int : Switch that allows for fixing q95 only in this constraint. - !! q95_fixed : input real : fixed safety factor q at 95% flux surface - use constraint_variables, only: fpsepbqar, psepbqarmax, i_q95_fixed, q95_fixed - use physics_variables, only: p_plasma_separatrix_mw, bt, q95, aspect, rmajor - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - if (i_q95_fixed == 1) then - tmp_cc = ((p_plasma_separatrix_mw*bt)/(q95_fixed*aspect*rmajor)) / psepbqarmax - 1.0d0 * fpsepbqar - tmp_err = (p_plasma_separatrix_mw*bt)/(q95_fixed*aspect*rmajor) - psepbqarmax - else - tmp_cc = ((p_plasma_separatrix_mw*bt)/(q95*aspect*rmajor)) / psepbqarmax - 1.0d0 * fpsepbqar - tmp_err = (p_plasma_separatrix_mw*bt)/(q95*aspect*rmajor) - psepbqarmax - end if - tmp_con = psepbqarmax - tmp_symbol = '<' - tmp_units = 'MWT/m' - - end subroutine constraint_eqn_068 - - subroutine constraint_eqn_069(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! TODO Remove - !! Ensure separatrix power is less than value from Kallenbach divertor - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Ensure separatrix power is less than value from Kallenbach divertor - !! #=# divertor_kallenbach - !! #=#=# consistency, psep_kallenbach - !! fpsep has been removed from the equation. - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! psep_kallenbach : input real : Power conducted through the separatrix, as calculated by the divertor model [W] - !! p_plasma_separatrix_mw : input real : power to conducted to the divertor region (MW) - ! use div_kal_vars, only: psep_kallenbach - use physics_variables, only: p_plasma_separatrix_mw - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - ! From Kallenbach model, should be reserved if the model is going to be added back - - end subroutine constraint_eqn_069 - - subroutine constraint_eqn_070(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! TODO Remove - !! Separatrix density consistency - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Separatrix density consistency - !! #=# divertor_kallenbach - !! #=#=# consistency - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! teomp : input real : Separatrix temperature calculated by the Kallenbach divertor model [eV] - !! tesep : input real : Electron temperature at separatrix [keV] - ! use div_kal_vars, only: teomp - use physics_variables, only: tesep - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - ! From Kallenbach model, should be reserved if the model is going to be added back - - end subroutine constraint_eqn_070 - - subroutine constraint_eqn_071(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! TODO Remove - !! Separatrix density consistency - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Separatrix density consistency - !! #=# divertor_kallenbach - !! #=#=# consistency - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! neomp : input real : Mean SOL density at OMP calculated by the Kallenbach divertor model [m-3] - !! nesep : input real : electron density at separatrix [m-3] (ipedestal=1,2, calculated if 3) - !! neratio : input real : Ratio of mean SOL density at OMP to separatrix density at OMP (iteration variable 121) - ! use div_kal_vars, only: neomp, neratio - use physics_variables, only: nesep - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - ! From Kallenbach model, should be reserved if the model is going to be added back - - end subroutine constraint_eqn_071 - - subroutine constraint_eqn_072(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Upper limit on central Solenoid Tresca yield stress - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Central Solenoid Tresca yield criterion - !! #=# pfcoil - !! #=#=# foh_stress, alstroh - !! In the case if the bucked and wedged option ( i_tf_bucking >= 2 ) the constrained - !! stress is the largest the largest stress of the - !! - CS stress at maximum current (conservative as the TF inward pressure is not taken - !! into account) - !! - CS stress at flux swing (no current in CS) from the TF inward pressure - !! This allow to cover the 2 worst stress scenario in the bucked and wedged design - !! Otherwise (free standing TF), the stress limits are only set by the CS stress at max current - !! Reverse the sign so it works as an inequality constraint (tmp_cc > 0) - !! This will have no effect if it is used as an equality constraint because it will be squared. - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! foh_stress : input real : f-value for Tresca yield criterion in Central Solenoid - !! alstroh : input real : allowable hoop stress in Central Solenoid structural material (Pa) - !! s_shear_cs_peak : input real : Maximum shear stress coils/central solenoid (Pa) - !! sig_tf_cs_bucked : input real : Maximum shear stress in CS case at flux swing (no current in CS) - !! can be significant for the bucked and weged design - !! i_tf_bucking : input integer : switch for TF structure design - use constraint_variables, only: foh_stress - use pfcoil_variables, only: alstroh, s_shear_cs_peak - use tfcoil_variables, only: sig_tf_cs_bucked, i_tf_bucking - use build_variables, only: i_tf_inside_cs - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - ! bucked and wedged desing (see subroutine comment) - if ( i_tf_bucking >= 2 .and. i_tf_inside_cs == 0 ) then - tmp_cc = max(s_shear_cs_peak, sig_tf_cs_bucked) / alstroh - 1.0d0 * foh_stress - tmp_err = alstroh - max(s_shear_cs_peak, sig_tf_cs_bucked) - ! Free standing CS - else - tmp_cc = s_shear_cs_peak / alstroh - 1.0d0 * foh_stress - tmp_err = alstroh - s_shear_cs_peak - end if - - tmp_con = alstroh - tmp_symbol = '<' - tmp_units = 'Pa' - - end subroutine constraint_eqn_072 - - subroutine constraint_eqn_073(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Lower limit to ensure separatrix power is greater than the L-H power + auxiliary power - !! Related to constraint 15 - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Ensure separatrix power is greater than the L-H power + auxiliary power - !! #=# physics - !! #=#=# fplhsep, p_plasma_separatrix_mw - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! fplhsep : input real : F-value for Psep >= Plh + Paux : for consistency of two values of separatrix power - !! p_l_h_threshold_mw : input real : L-H mode power threshold (MW) - !! p_plasma_separatrix_mw : input real : power to be conducted to the divertor region (MW) - !! p_hcd_injected_total_mw : inout real : total auxiliary injected power (MW) - use physics_variables, only: fplhsep, p_l_h_threshold_mw, p_plasma_separatrix_mw - use current_drive_variables, only: p_hcd_injected_total_mw - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = 1.0d0 - fplhsep * p_plasma_separatrix_mw / (p_l_h_threshold_mw+p_hcd_injected_total_mw) - tmp_con = p_plasma_separatrix_mw - tmp_err = p_plasma_separatrix_mw * tmp_cc - tmp_symbol = '>' - tmp_units = 'MW' - - end subroutine constraint_eqn_073 - - subroutine constraint_eqn_074(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Upper limit to ensure TF coil quench temperature < tmax_croco - !! ONLY used for croco HTS coil - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Ensure TF coil quench temperature < tmax_croco ONLY used for croco HTS coil - !! #=# physics - !! #=#=# fcqt, tmax_croco - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! fcqt : input real : f-value: TF coil quench temparature remains below tmax_croco - !! croco_quench_temperature : input real : CroCo strand: Actual temp reached during a quench (K) - !! tmax_croco : input real : CroCo strand: maximum permitted temp during a quench (K) - use constraint_variables, only: fcqt - use tfcoil_variables, only: croco_quench_temperature, tmax_croco - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = croco_quench_temperature / tmax_croco - 1.0d0 * fcqt - tmp_con = croco_quench_temperature - tmp_err = croco_quench_temperature * tmp_cc - tmp_symbol = '<' - tmp_units = 'K' - - end subroutine constraint_eqn_074 - - subroutine constraint_eqn_075(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Upper limit to ensure that TF coil current / copper area < Maximum value - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Ensure that TF coil current / copper area < Maximum value - !! ONLY used for croco HTS coil - !! #=# physics - !! #=#=# f_coppera_m2, copperA_m2_max - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! copperA_m2 : input real : TF coil current / copper area (A/m2) - !! copperA_m2_max : input real : Maximum TF coil current / copper area (A/m2) - !! f_coppera_m2 : input real : f-value for TF coil current / copper area < copperA_m2_max - use rebco_variables, only: copperA_m2, copperA_m2_max, f_coppera_m2 - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = copperA_m2 / copperA_m2_max - 1.0d0 * f_coppera_m2 - tmp_con = copperA_m2 - tmp_err = copperA_m2 * tmp_cc - tmp_symbol = '<' - tmp_units = 'A/m2' - - end subroutine constraint_eqn_075 - - subroutine constraint_eqn_076(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Upper limit for Eich critical separatrix density model: Added for issue 558 - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Eich critical separatrix density model - !! Added for issue 558 with ref to http://iopscience.iop.org/article/10.1088/1741-4326/aaa340/pdf - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! alpha_crit : output real : critical ballooning parameter value - !! nesep_crit : output real : critical electron density at separatrix [m-3] - !! kappa : input real : plasma separatrix elongation (calculated if i_plasma_geometry = 1-5, 7 or 9) - !! triang : input real : plasma separatrix triangularity (calculated if i_plasma_geometry = 1, 3-5 or 7) - !! aspect : input real : aspect ratio (iteration variable 1) - !! p_plasma_separatrix_mw : input real : power to conducted to the divertor region (MW) - !! dlimit(7) : input real array : density limit (/m3) as calculated using various models - !! fnesep : input real : f-value for Eich critical separatrix density - use physics_variables, only: alpha_crit, nesep_crit, kappa, triang, & - aspect, p_plasma_separatrix_mw, dlimit, nesep - use constraint_variables, only: fnesep - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - alpha_crit = (kappa ** 1.2D0) * (1.0D0 + 1.5D0 * triang) - nesep_crit = 5.9D0 * alpha_crit * (aspect ** (-2.0D0/7.0D0)) * & - (((1.0D0 + (kappa ** 2.0D0)) / 2.0D0) ** (-6.0D0/7.0D0)) & - * ((p_plasma_separatrix_mw* 1.0D6) ** (-11.0D0/70.0D0)) * dlimit(7) - tmp_cc = nesep / nesep_crit - 1.0D0 * fnesep - tmp_con = nesep - tmp_err = nesep * tmp_cc - tmp_symbol = '<' - tmp_units = 'm-3' - - end subroutine constraint_eqn_076 - - subroutine constraint_eqn_077(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for maximum TF current per turn upper limit - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; residual error in physical units; output string; units string - !! Equation for maximum TF current per turn upper limit - !! #=# tfcoil - !! #=#=# fcpttf, c_tf_turn, cpttf_max - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! fcpttf : input : f-value for TF coil current per turn - !! cpttf_max : input : allowable TF coil current per turn [A/turn] - !! c_tf_turn : input : TF coil current per turn [A/turn] - use constraint_variables, only: fcpttf - use tfcoil_variables, only: cpttf_max, c_tf_turn - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = c_tf_turn / cpttf_max - 1.0D0 * fcpttf - tmp_con = cpttf_max - tmp_err = cpttf_max * tmp_cc - tmp_symbol = '<' - tmp_units = 'A/turn' - - end subroutine constraint_eqn_077 - - subroutine constraint_eqn_078(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for Reinke criterion, divertor impurity fraction lower limit - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; residual error in physical units; output string; units string - !! Equation for Reinke criterion, divertor impurity fraction lower limit - !! #=# divertor - !! #=#=# freinke, fzactual, fzmin - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present; - !! and con will be printed out only if present. Thesw conditions were missing. - !! freinke : input : f-value for Reinke criterion (itv 147) - !! fzmin : input : minimum impurity fraction from Reinke model - !! fzactual : input : actual impurity fraction - use constraint_variables, only: freinke - use reinke_variables, only: fzactual, fzmin - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - ! write(*,*) 'freinke, fzact, fzmin = ', freinke, ', ', fzactual, ', ', fzmin - ! 1.0, 0.0, value - tmp_cc = 1.0D0 - freinke * fzactual/fzmin - !The following two pre-existing lines are not understood: - !KE note - cc is always 1, code never enters IF statement... - tmp_con = fzmin * (1.0D0 - tmp_cc) - tmp_err = fzmin * tmp_cc - tmp_symbol = '>' - tmp_units = '' - ! write(*,*) 'cc, con = ', tmp_cc, ', ', tmp_con - ! write(*,*) 'freinke, fzactual, fzmin = ', freinke, ', ', fzactual, ', ', fzmin - - end subroutine constraint_eqn_078 - - subroutine constraint_eqn_079(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for maximum CS field - !! author: P B Lloyd, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; residual error in physical units; output string; units string - !! Equation for maximum CS field - !! #=# pfcoil - !! #=#=# fb_cs_limit_max, b_cs_peak_flat_top_end, b_cs_peak_pulse_start, b_cs_limit_max - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! fb_cs_limit_max : input : F-value for CS mmax field (cons. 79, itvar 149) - !! b_cs_limit_max : input : Central solenoid max field limit [T] - !! b_cs_peak_pulse_start : input : maximum field in central solenoid at beginning of pulse (T) - !! b_cs_peak_flat_top_end : input real : maximum field in central solenoid at end of flat-top (EoF) (T) - !! (Note: original code has "b_cs_peak_flat_top_end/b_cs_peak_pulse_start | peak CS field [T]".) - use pfcoil_variables, only: fb_cs_limit_max, b_cs_limit_max, b_cs_peak_pulse_start, b_cs_peak_flat_top_end - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = max(b_cs_peak_flat_top_end, b_cs_peak_pulse_start) / b_cs_limit_max - 1.0D0 * fb_cs_limit_max - tmp_con = b_cs_limit_max - tmp_err = max(b_cs_peak_flat_top_end, b_cs_peak_pulse_start) * tmp_cc - tmp_symbol = '<' - tmp_units = 'A/turn' - - end subroutine constraint_eqn_079 - - subroutine constraint_eqn_080(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for p_plasma_separatrix_mw lower limit - !! author: J Morris, Culham Science Centre - !! args : output structure : residual error; constraint value; residual error in physical units; - !! output string; units string - !! Lower limit p_plasma_separatrix_mw - !! #=# physics - !! #=#=# fp_plasma_separatrix_min_mw, p_plasma_separatrix_mw - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! fp_plasma_separatrix_min_mw : input : F-value for lower limit on p_plasma_separatrix_mw (cons. 80, itvar 153) - !! p_plasma_separatrix_min_mw : input : Minimum power crossing separatrix p_plasma_separatrix_mw [MW] - !! p_plasma_separatrix_mw : input : Power crossing separatrix [MW] - use physics_variables, only: fp_plasma_separatrix_min_mw, p_plasma_separatrix_mw - use constraint_variables, only : p_plasma_separatrix_min_mw - implicit none - - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - tmp_cc = 1.0D0 - fp_plasma_separatrix_min_mw * p_plasma_separatrix_mw / p_plasma_separatrix_min_mw - tmp_con = p_plasma_separatrix_min_mw - tmp_err = p_plasma_separatrix_mw * tmp_cc - tmp_symbol = '>' - tmp_units = 'MW' - - end subroutine constraint_eqn_080 - - subroutine constraint_eqn_081(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Lower limit to ensure central density is larger that the pedestal one - !! author: S Kahn, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Lower limit ne0 > neped - !! !#=# physics - !! !#=#=# ne0, neped - !! Logic change during pre-factoring: err, symbol, units will be - !! assigned only if present. - !! fne0 : input : F-value for constraint on ne0 > neped - !! ne0 : input : Central electron density [m-3] - !! neped : input : Electron density at pedestal [m-3] - use physics_variables, only: ne0, fne0, neped - implicit none - - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - tmp_cc = 1.0D0 - fne0 * ne0/neped - tmp_con = fne0 - tmp_err = fne0 * tmp_cc - tmp_symbol = '>' - tmp_units = '/m3' - - end subroutine constraint_eqn_081 - - subroutine constraint_eqn_082(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for toroidal consistency of stellarator build - !! author: J Lion, IPP Greifswald - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! toroidalgap > dx_tf_inboard_out_toroidal - !! #=# tfcoil - !! #=#=# dx_tf_inboard_out_toroidal, ftoroidalgap - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! ftoroidalgap : input real : f-value for constraint toroidalgap > dx_tf_inboard_out_toroidal - !! toroidalgap : input real : minimal gap between two stellarator coils - !! dx_tf_inboard_out_toroidal : input real : total toroidal width of a tf coil - use tfcoil_variables, only: dx_tf_inboard_out_toroidal,ftoroidalgap,toroidalgap - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = 1.0D0 - ftoroidalgap * toroidalgap/dx_tf_inboard_out_toroidal - tmp_con = toroidalgap - tmp_err = toroidalgap - dx_tf_inboard_out_toroidal/ftoroidalgap - tmp_symbol = '<' - tmp_units = 'm' - - end subroutine constraint_eqn_082 - - subroutine constraint_eqn_083(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for radial consistency of stellarator build - !! author: J Lion, IPP Greifswald - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! available_radial_space > required_radial_space - !! #=# build - !! #=#=# required_radial_space, f_avspace - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! f_avspace : input real : f-value for constraint available_radial_space > required_radial_space - !! available_radial_space : input real : avaible space in radial direction as given by each s.-configuration - !! required_radial_space : input real : required space in radial direction - use build_variables, only: available_radial_space, required_radial_space, f_avspace - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = 1.0D0 - f_avspace * available_radial_space/required_radial_space - tmp_con = available_radial_space * (1.0D0 - tmp_cc) - tmp_err = required_radial_space * tmp_cc - tmp_symbol = '<' - tmp_units = 'm' - end subroutine constraint_eqn_083 - - subroutine constraint_eqn_084(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for the lower limit of beta - !! author: J Lion, IPP Greifswald - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! (beta-beta_fast_alpha) > beta_min - !! #=# physics - !! #=#=# beta_fast_alpha, beta, fbeta_min - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! fbeta_min : input real : f-value for constraint beta-beta_fast_alpha > beta_min - !! beta_min : input real : Lower limit for beta - !! beta : input real : plasma beta - - use physics_variables, only: beta_min, beta - use constraint_variables, only: fbeta_min - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - - tmp_cc = 1.0D0 - fbeta_min * (beta)/beta_min - tmp_con = beta_min * (1.0D0 - tmp_cc) - tmp_err = (beta) * tmp_cc - tmp_symbol = '>' - tmp_units = '' - - - end subroutine constraint_eqn_084 - - subroutine constraint_eqn_085(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equality constraint for the centerpost (CP) lifetime - !! Author : S Kahn - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Depending on the chosen option : i_cp_lifetime - !! - 0 : The CP full power year lifelime is set by the user (cplife_input) - !! - 1 : The CP lifelime is equal to the divertor one - !! - 2 : The CP lifetime is equal to the breeding blankets one - !! - 3 : The CP lifetime is equal to the plant one - !! #=# availability - !! #=#=# consistency - !! Logic change during pre-factoring: err, symbol, units will be assigned - !! only if present. - !! cplife : input real : calculated CP full power year lifetime (years) - !! life_blkt_fpy : input real : calculated first wall/blanket power year lifetime (years) - !! divlife : input real : calculated divertor power year lifetime (years) - !! i_cp_lifetime : input integer : switch chosing which plant element the CP - !! the CP lifetime must equate - use cost_variables, only : cplife, divlife, cplife_input, & - tlife, i_cp_lifetime - use fwbs_variables, only : life_blkt_fpy - - implicit none - - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - !! Constraints output - - - ! The CP lifetime is equal to the the divertor one - if ( i_cp_lifetime == 0 ) then - tmp_cc = 1.0D0 - cplife/cplife_input - - else if ( i_cp_lifetime == 1 ) then - tmp_cc = 1.0D0 - cplife/divlife - - ! The CP lifetime is equal to the tritium breeding blankets / FW one - else if ( i_cp_lifetime == 2 ) then - tmp_cc = 1.0D0 - cplife/life_blkt_fpy - - ! The CP lifetime is equal to the - else if ( i_cp_lifetime == 3 ) then - tmp_cc = 1.0D0 - cplife/tlife - end if - - tmp_con = divlife * (1.0D0 - tmp_cc) - tmp_err = divlife * tmp_cc - tmp_symbol = '=' - tmp_units = 'years' - - end subroutine constraint_eqn_085 - - subroutine constraint_eqn_086(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Upper limit on the turn edge length in the TF winding pack - !! Author : S Kahn - !! args : output structure : residual error; constraint value; - !! residual error in physical units; - !! t_turn_tf : input real : TF coil turn edge length including turn insulation [m] - !! f_t_turn_tf : input real : f-value for TF turn edge length constraint - !! t_turn_tf_max : input real : TF turn edge length including turn insulation upper limit [m] - use tfcoil_variables, only : t_turn_tf, f_t_turn_tf, t_turn_tf_max - - implicit none - - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - !! Constraints output - tmp_cc = t_turn_tf / t_turn_tf_max - 1.0D0 * f_t_turn_tf - tmp_con = t_turn_tf_max * (1.0D0 - tmp_cc) - tmp_err = t_turn_tf_max * tmp_cc - tmp_symbol = '<' - tmp_units = 'm' - - end subroutine constraint_eqn_086 - - - subroutine constraint_eqn_087(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for TF coil cryogenic power upper limit - !! author: S. Kahn, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! p_cryo_plant_electric_mw : input real : cryogenic plant power (MW) - !! f_crypmw : input real : f-value for maximum cryogenic plant power - !! p_cryo_plant_electric_max_mw : input real : Maximum cryogenic plant power (MW) - - use heat_transport_variables, only: p_cryo_plant_electric_mw, p_cryo_plant_electric_max_mw, f_crypmw - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = p_cryo_plant_electric_mw / p_cryo_plant_electric_max_mw - 1.0D0 * f_crypmw - tmp_con = p_cryo_plant_electric_max_mw * (1.0D0 - tmp_cc) - tmp_err = p_cryo_plant_electric_mw * tmp_cc - tmp_symbol = '<' - tmp_units = 'MW' - end subroutine constraint_eqn_087 - - subroutine constraint_eqn_088(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for TF coil vertical strain upper limit (absolute value) - !! author: CPS Swanson, PPPL, USA - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! Equation for TF coil vertical strain upper limit (absolute value) - !! #=# tfcoil - !! #=#=# fstr_wp, str_wp_max - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! fstr_wp : input real : f-value for TF coil strain - !! str_wp_max : input real : Allowable maximum TF coil vertical strain - !! str_wp : input real : Constrained TF coil vertical strain - use constraint_variables, only: fstr_wp - use tfcoil_variables, only: str_wp_max, str_wp - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = abs(str_wp) / str_wp_max - 1.0D0 * fstr_wp - tmp_con = str_wp_max - tmp_err = str_wp_max - abs(str_wp) / fstr_wp - tmp_symbol = '<' - tmp_units = '' - end subroutine constraint_eqn_088 - - subroutine constraint_eqn_089(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Upper limit to ensure that the Central Solenoid [OH] coil current / copper area < Maximum value - !! author: G Turkington, CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! #=# physics - !! #=#=# f_copperaoh_m2, copperaoh_m2_max - !! and hence also optional here. - !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. - !! copperaoh_m2 : input real : CS coil current at EOF / copper area [A/m2] - !! copperaoh_m2_max : input real : maximum coil current / copper area [A/m2] - !! f_copperaoh_m2 : input real : f-value for CS coil current / copper area - use rebco_variables, only: copperaoh_m2, copperaoh_m2_max, f_copperaoh_m2 - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - tmp_cc = copperaoh_m2 / copperaoh_m2_max - 1.0d0 * f_copperaoh_m2 - tmp_con = copperaoh_m2 - tmp_err = copperaoh_m2 * tmp_cc - tmp_symbol = '<' - tmp_units = 'A/m2' - - end subroutine constraint_eqn_089 - - subroutine constraint_eqn_090(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Lower limit for CS coil stress load cycles - !! author: A. Pearce, G Turkington CCFE, Culham Science Centre - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! fncycle : input real : f-value for constraint n_cycle > n_cycle_min - !! n_cycle : input real : Allowable number of cycles for CS - !! n_cycle_min : input real : Minimum required cycles for CS - use CS_fatigue_variables, only: n_cycle, n_cycle_min, bkt_life_csf - use constraint_variables, only: fncycle - use cost_variables, only: ibkt_life, bktcycles - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - !! Switch to relay the calculated fw/blanket lifetime cycles as the minimum required CS stress cycles. - !! bkt_life_cycle = 1 turns on the relay. Otherwise the models run independently. - if (ibkt_life == 1 .and. bkt_life_csf == 1 ) then - n_cycle_min = bktcycles - end if - - tmp_cc = 1.0D0 - fncycle * n_cycle / n_cycle_min - tmp_con = n_cycle_min * (1.0D0 - tmp_cc) - tmp_err = n_cycle * tmp_cc - tmp_symbol = '>' - tmp_units = '' - - end subroutine constraint_eqn_090 - - subroutine constraint_eqn_091(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Lower limit to ensure ECRH te is greater than required te for ignition - !! at lower values for n and B. Or if the design point is ECRH heatable (if i_plasma_ignited==0) - !! stellarators only (but in principle usable also for tokamaks). - !! author: J Lion, IPP Greifswald - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! powerht_local > powerscaling - !! #=# physics - !! #=#=# fecrh_ignition, powerht_local, powerscaling - !! fecrh_ignition : input real : f-value for constraint powerht_local > powerscaling - !! max_gyrotron_frequency : input real : Max. av. gyrotron frequency - !! te0_ecrh_achievable : input real : Max. achievable electron temperature at ignition point - use constraint_variables, only: fecrh_ignition - use stellarator_variables, only: max_gyrotron_frequency, te0_ecrh_achievable, powerscaling_constraint, powerht_constraint - use physics_variables, only: i_plasma_ignited - use current_drive_variables, only: p_hcd_primary_extra_heat_mw - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - ! Achievable ECRH te needs to be larger than needed te for igntion - if(i_plasma_ignited==0) then - tmp_cc = 1.0D0 - fecrh_ignition* (powerht_constraint+p_hcd_primary_extra_heat_mw)/powerscaling_constraint - else - tmp_cc = 1.0D0 - fecrh_ignition* powerht_constraint/powerscaling_constraint - endif - - tmp_con = powerscaling_constraint * (1.0D0 - tmp_cc) - tmp_err = powerht_constraint * tmp_cc - tmp_symbol = '<' - tmp_units = 'MW' - end subroutine constraint_eqn_091 - - subroutine constraint_eqn_092(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) - !! Equation for checking is D/T ratio is consistent, and sums to 1. - !! author: G Turkington, UKAEA - !! args : output structure : residual error; constraint value; - !! residual error in physical units; output string; units string - !! f_deuterium : input : fraction of deuterium ions - !! f_tritium : input : fraction of tritium ions - !! f_helium3 : input : fraction of helium-3 ions - use physics_variables, only: f_deuterium, f_tritium, f_helium3 - implicit none - real(dp), intent(out) :: tmp_cc - real(dp), intent(out) :: tmp_con - real(dp), intent(out) :: tmp_err - character(len=1), intent(out) :: tmp_symbol - character(len=10), intent(out) :: tmp_units - - - ! Iterate over f_tritium and calculate f_deuterium - f_deuterium = 1.0D0 - (f_tritium + f_helium3) - tmp_cc = 1.0D0 - (f_deuterium + f_tritium + f_helium3) - tmp_con = 1.0D0 - tmp_err = tmp_con * tmp_cc - tmp_symbol = '=' - tmp_units = 'fraction' - - end subroutine constraint_eqn_092 - - -end module constraints diff --git a/source/fortran/constraint_variables.f90 b/source/fortran/constraint_variables.f90 index 31824dccb3..5d90917269 100644 --- a/source/fortran/constraint_variables.f90 +++ b/source/fortran/constraint_variables.f90 @@ -310,99 +310,4 @@ module constraint_variables real(dp) :: fcqt !! TF coil quench temparature remains below tmax_croco !! (`constraint equation 74`, `iteration variable 141`) - - contains - - subroutine init_constraint_variables - !! Initialise module variables - implicit none - - auxmin = 0.1D0 - beta_poloidal_max = 0.19D0 - bigqmin = 10.0D0 - bmxlim = 12.0D0 - fauxmn = 1.0D0 - fbeta_poloidal_eps = 1.0D0 - fbeta_poloidal = 1.0D0 - fbeta_max = 1.0D0 - fbeta_min = 1.0D0 - fcpttf = 1.0D0 - fr_conducting_wall = 1.0D0 - fdene = 1.0D0 - fdtmp = 1.0D0 - fflutf = 1.0D0 - fp_fusion_total_max_mw = 1.0D0 - fgamcd = 1.0D0 - fpflux_div_heat_load_mw = 1.0D0 - fiooic = 0.5D0 - fipir = 1.0D0 - q95_fixed = 3.0D0 - fjohc = 1.0D0 - fjohc0 = 1.0D0 - fjprot = 1.0D0 - fl_h_threshold = 1.0D0 - fmva = 1.0D0 - fnbshinef = 1.0D0 - fncycle = 1.0D0 - fnesep = 1.0D0 - foh_stress = 1.0D0 - fpeakb = 1.0D0 - fpinj = 1.0D0 - fpnetel = 1.0D0 - fportsz = 1.0D0 - fpsepbqar = 1.0D0 - fpsepr = 1.0D0 - fptemp = 1.0D0 - fptfnuc = 1.0D0 - fq = 1.0D0 - fqval = 1.0D0 - fradpwr = 0.99D0 - fradwall = 1.0D0 - freinke = 1.0D0 - frminor = 1.0D0 - fstrcase = 1.0D0 - fstrcond = 1.0D0 - fstr_wp = 1.0D0 - fmaxvvstress = 1.0D0 - ftbr = 1.0D0 - ft_burn = 1.0D0 - ftcycl = 1.0D0 - ftmargoh = 1.0D0 - ftmargtf = 1.0D0 - ft_current_ramp_up = 1.0D0 - ftpeak = 1.0D0 - fvdump = 1.0D0 - fvs = 1.0D0 - fvvhe = 1.0D0 - fwalld = 1.0D0 - fzeffmax = 1.0D0 - gammax = 2.0D0 - i_q95_fixed = 0 - pflux_fw_rad_max = 1.0D0 - mvalim = 40.0D0 - f_p_beam_shine_through_max = 1.0D-3 - nflutfmax = 1.0D23 - p_plasma_separatrix_min_mw = 150.0D0 - f_fw_rad_max = 3.33D0 - pflux_fw_rad_max_mw = 0.0D0 - pnetelin = 1.0D3 - p_fusion_total_max_mw = 1.5D3 - psepbqarmax = 9.5D0 - pseprmax = 25.0D0 - ptfnucmax = 1.0D-3 - tbrmin = 1.1D0 - t_burn_min = 1.0D0 - tcycmn = 0.0D0 - t_current_ramp_up_min = 1.0D0 - vvhealw = 1.0D0 - walalw = 1.0D0 - f_alpha_energy_confinement_min = 5.0D0 - falpha_energy_confinement = 1.0D0 - fniterpump = 1.0D0 - zeffmax = 3.6D0 - fpoloidalpower = 1.0D0 - fpsep = 1.0D0 - fcqt = 1.0D0 - fecrh_ignition = 1.0D0 - end subroutine init_constraint_variables end module constraint_variables