Update PROCESS v3 (OS) mappings for EU-DEMO (#2749)

* first pass mappings

* long_name not long name

* populate PROCESS params

* unit issue

* finish process params descriptions

* 💩 Enable really bad units

* mapping fixes

* start on special ase mechanics

* shield + vv shenanigans

* finish applying hacks

* add q_0 to test params

* add mslimit to test data

* add sigma cs

* add missing params to test data

* move q out of OUT mapping

* add proper TODOO

* ✨ Allow different in and out names for external codes (#2750)

* ✨ Allow different in and out names for external codes

* ✅ Add tests

* 🚧 A horrible hack for value=None unit=crap

* 🍻 Actually it might be ok

* manually turn shldith send back on

---------

Co-authored-by: james <james.cook1@ukaea.uk>
Co-authored-by: je-cook <81617086+je-cook@users.noreply.github.com>

bluemira/base/parameter_frame/_frame.py

@@ -405,6 +405,20 @@ def _validate_parameter_field(field, member_type: Type) -> Tuple[Type, ...]:
 def _validate_units(param_data: Dict, value_type: Iterable[Type]):
     try:
         quantity = pint.Quantity(param_data["value"], param_data["unit"])
+    except ValueError:
+        try:
+            quantity = pint.Quantity(f'{param_data["value"]}*{param_data["unit"]}')
+        except pint.errors.PintError as pe:
+            if param_data["value"] is None:
+                quantity = pint.Quantity(
+                    1 if param_data["unit"] in (None, "") else param_data["unit"]
+                )
+                param_data["source"] = f"{param_data.get('source', '')}\nMAD UNIT 🤯 😭:"
+            else:
+                raise ValueError("Unit conversion failed") from pe
+        else:
+            param_data["value"] = quantity.magnitude
+        param_data["unit"] = quantity.units
     except KeyError as ke:
         raise ValueError("Parameters need a value and a unit") from ke
     except TypeError:
@@ -434,6 +448,9 @@ def _validate_units(param_data: Dict, value_type: Iterable[Type]):
 
     param_data["unit"] = f"{unit:~P}"
 
+    if "MAD UNIT" in param_data.get("source", ""):
+        param_data["source"] += f"{quantity.magnitude}{param_data['unit']}"
+
 
 def _remake_units(dimensionality: Union[Dict, pint.util.UnitsContainer]) -> pint.Unit:
     """Reconstruct unit from its dimensionality"""

bluemira/codes/interface.py

@@ -236,9 +236,11 @@ def _map_external_outputs_to_bluemira_params(
         for bm_name, mapping in self.params.mappings.items():
             if not (mapping.recv or recv_all):
                 continue
-            output_value = self._get_output_or_raise(external_outputs, mapping.name)
+            output_value = self._get_output_or_raise(external_outputs, mapping.out_name)
             if mapping.unit is None:
-                bluemira_warn(f"{mapping.name} from code {self._name} has no known unit")
+                bluemira_warn(
+                    f"{mapping.out_name} from code {self._name} has no known unit"
+                )
                 value = output_value
             elif output_value is None:
                 value = output_value

bluemira/codes/params.py

@@ -132,6 +132,7 @@ class ParameterMapping:
     """
 
     name: str
+    out_name: Optional[str] = None
     send: bool = True
     recv: bool = True
     unit: Optional[str] = None
@@ -142,14 +143,17 @@ def __post_init__(self):
         """
         Freeze the dataclass
         """
-        self._frozen = ("name", "unit", "_frozen")
+        if self.out_name is None:
+            self.out_name = self.name
+        self._frozen = ("name", "out_name", "unit", "_frozen")
 
     def to_dict(self) -> Dict:
         """
         Convert this object to a dictionary with attributes as values.
         """
         return {
             "name": self.name,
+            "out_name": self.out_name,
             "send": self.send,
             "recv": self.recv,
             "unit": self.unit,
@@ -181,10 +185,10 @@ def __setattr__(self, attr: str, value: Union[bool, str]):
             Value of attribute
         """
         if (
-            attr not in ["send", "recv", "name", "unit", "_frozen"]
+            attr not in ["send", "recv", "name", "out_name", "unit", "_frozen"]
             or attr in self._frozen
         ):
-            raise KeyError(f"{attr} cannot be set for a {self.__class__.__name__}")
+            raise KeyError(f"{attr} cannot be set for a {type(self).__name__}")
         if attr in ["send", "recv"] and not isinstance(value, bool):
             raise ValueError(f"{attr} must be a bool")
         super().__setattr__(attr, value)

bluemira/codes/process/mapping.py

@@ -25,28 +25,33 @@
 from bluemira.codes.utilities import create_mapping
 
 IN_mappings = {
-    "P_el_net": ("pnetelin", "MW"),
     "n_TF": ("n_tf", "dimensionless"),
     "TF_ripple_limit": ("ripmax", "%"),
     "C_Ejima": ("gamma", "dimensionless"),
     "e_nbi": ("enbeam", "keV"),
     "P_hcd_ss": ("pinjalw", "MW"),
     "eta_nb": ("etanbi", "dimensionless"),
     "e_mult": ("emult", "dimensionless"),
-    "tk_sh_out": ("shldoth", "m"),
-    "tk_sh_top": ("shldtth", "m"),
-    "tk_sh_bot": ("shldlth", "m"),
-    "tk_vv_out": ("d_vv_out", "m"),
-    "tk_vv_top": ("d_vv_top", "m"),
-    "tk_vv_bot": ("d_vv_bot", "m"),
     "tk_cr_vv": ("ddwex", "m"),
     "tk_tf_front_ib": ("casthi", "m"),
     "tk_tf_side": ("casths", "m"),
     "PsepB_qAR_max": ("psepbqarmax", "MW.T/m"),
+    "q_0": ("q0", "dimensionless"),
+    "m_s_limit": ("m_s_limit", "dimensionless"),
+    "delta": ("triang", "dimensionless"),
+    "sigma_tf_case_max": ("sig_tf_case_max", "Pa"),
+    "sigma_tf_wp_max": ("sig_tf_wp_max", "Pa"),
+    "sigma_cs_wp_max": ("alstroh", "Pa"),
+    "H_star": ("hfact", "dimensionless"),
+    "bb_pump_eta_el": ("etahtp", "dimensionless"),
+    "bb_pump_eta_isen": ("etaiso", "dimensionless"),
+    "bb_t_inlet": ("inlet_temp", "K"),
+    "bb_t_outlet": ("outlet_temp", "K"),
+    "eta_ecrh": ("etaech", "dimensionless"),
+    "gamma_ecrh": ("gamma_ecrh", "1e20 A/W/m^2"),
 }
 
 OUT_mappings = {
-    "P_el_net_process": ("pnetelmw", "MW"),
     "R_0": ("rmajor", "m"),
     "B_0": ("bt", "T"),
     "kappa_95": ("kappa95", "dimensionless"),
@@ -101,36 +106,47 @@
 
 IO_mappings = {
     "A": ("aspect", "dimensionless"),
+    "tau_flattop": (("tbrnmn", "tburn"), "s"),
+    "P_el_net": (("pnetelin", "pnetelmw"), "MW"),
     "tk_bb_ib": ("blnkith", "m"),
     "tk_bb_ob": ("blnkoth", "m"),
-    "tk_sh_in": ("shldith", "m"),
     "tk_vv_in": ("d_vv_in", "m"),
     "tk_sol_ib": ("scrapli", "m"),
     "tk_sol_ob": ("scraplo", "m"),
-    # Thermal shield thickness is a constant for us
-    # "tk_ts": ("thshield_ob", "m"),
-    # "tk_ts": ("thshield_vb", "m"),
-    "tk_ts": ("thshield_ib", "m"),
     "g_cs_tf": ("gapoh", "m"),
     "g_ts_tf": ("tftsgap", "m"),
     "g_vv_bb": ("vvblgap", "m"),
 }
 
 NONE_mappings = {
-    "tau_flattop": ("tburn", "s"),
     "B_tf_peak": ("bmaxtfrp", "T"),
-    "q_95": ("q95", "dimensionless"),
     "T_e": ("te", "keV"),
     "Z_eff": ("zeff", "amu"),
     "V_p": ("vol", "m^3"),
     "l_i": ("rli", "dimensionless"),
     "f_ni": ("faccd", "dimensionless"),
     "tk_tf_outboard": ("tfthko", "m"),
-    "sigma_tf_case_max": ("sig_tf_case_max", "Pa"),
-    "sigma_tf_wp_max": ("sig_tf_wp_max", "Pa"),
     "h_cp_top": ("h_cp_top", "m"),
     "h_tf_max_in": ("hmax", "m"),
     "r_tf_inboard_out": ("r_tf_inboard_out", "m"),
+    # The following mappings are not 1:1
+    "tk_sh_in": ("shldith", "m"),
+    "tk_sh_out": ("shldoth", "m"),
+    "tk_sh_top": ("shldtth", "m"),
+    "tk_sh_bot": ("shldlth", "m"),
+    "tk_vv_out": ("d_vv_out", "m"),
+    "tk_vv_top": ("d_vv_top", "m"),
+    "tk_vv_bot": ("d_vv_bot", "m"),
+    # Thermal shield thickness is a constant for us
+    "tk_ts": ("thshield_ib", "m"),
+    # "tk_ts": ("thshield_ob", "m"),
+    # "tk_ts": ("thshield_vb", "m"),
+    # TODO: q is not properly put in the MFILE output
+    # This should be ok OK most of the time as q_95 is input and then
+    # used as the lower bound of the q iteration variable, but this
+    # should be fixed as soon as PROCESS deal with this issue on
+    # their side
+    "q_95": ("q", "dimensionless"),
 }
 
 mappings = create_mapping(IN_mappings, OUT_mappings, IO_mappings, NONE_mappings)

bluemira/codes/process/params.py

@@ -64,6 +64,9 @@ class ProcessSolverParams(MappedParameterFrame):
     P_el_net: Parameter[float]
     """Net electrical power output [megawatt]."""
 
+    tau_flattop: Parameter[float]
+    """Flat-top duration [second]."""
+
     P_hcd_ss: Parameter[float]
     """Steady-state HCD power [megawatt]."""
 
@@ -100,6 +103,48 @@ class ProcessSolverParams(MappedParameterFrame):
     PsepB_qAR_max: Parameter[float]
     """Maximum PsepB/q95AR vale [MW.T/m]"""
 
+    q_0: Parameter[float]
+    """Plasma safety factor on axis [dimensionless]"""
+
+    q_95: Parameter[float]
+    """Plasma safety factor at the 95th percentile flux surface [dimensionless]"""
+
+    m_s_limit: Parameter[float]
+    """Margin to vertical stability [dimensionless]"""
+
+    delta: Parameter[float]
+    """Triangularity [dimensionless]"""
+
+    sigma_tf_case_max: Parameter[float]
+    """Maximum von Mises stress in the TF coil case nose [pascal]."""
+
+    sigma_tf_wp_max: Parameter[float]
+    """Maximum von Mises stress in the TF coil winding pack [pascal]."""
+
+    sigma_cs_wp_max: Parameter[float]
+    """Maximum von Mises stress in the CS coil winding pack [pascal]."""
+
+    H_star: Parameter[float]
+    """H factor (radiation corrected) [dimensionless]."""
+
+    bb_pump_eta_el: Parameter[float]
+    """Breeding blanket pumping electrical efficiency [dimensionless]"""
+
+    bb_pump_eta_isen: Parameter[float]
+    """Breeding blanket pumping isentropic efficiency [dimensionless]"""
+
+    bb_t_inlet: Parameter[float]
+    """Breeding blanket inlet temperature [K]"""
+
+    bb_t_outlet: Parameter[float]
+    """Breeding blanket outlet temperature [K]"""
+
+    eta_ecrh: Parameter[float]
+    """Electron cyclotron resonce heating wallplug efficiency [dimensionless]"""
+
+    gamma_ecrh: Parameter[float]
+    """Electron cyclotron resonce heating current drive efficiency [TODO: UNITS!]"""
+
     # Out parameters
     B_0: Parameter[float]
     """Toroidal field at R_0 [tesla]."""
@@ -116,7 +161,6 @@ class ProcessSolverParams(MappedParameterFrame):
     delta_95: Parameter[float]
     """95th percentile plasma triangularity [dimensionless]."""
 
-    delta: Parameter[float]
     """Last closed surface plasma triangularity [dimensionless]."""
 
     f_bs: Parameter[float]
@@ -125,9 +169,6 @@ class ProcessSolverParams(MappedParameterFrame):
     g_vv_ts: Parameter[float]
     """Gap between VV and TS [meter]."""
 
-    H_star: Parameter[float]
-    """H factor (radiation corrected) [dimensionless]."""
-
     I_p: Parameter[float]
     """Plasma current [megaampere]."""
 
@@ -143,9 +184,6 @@ class ProcessSolverParams(MappedParameterFrame):
     P_brehms: Parameter[float]
     """Bremsstrahlung [megawatt]."""
 
-    P_el_net_process: Parameter[float]
-    """Net electrical power output as provided by PROCESS [megawatt]."""
-
     P_fus_DD: Parameter[float]
     """D-D fusion power [megawatt]."""
 
@@ -307,24 +345,12 @@ class ProcessSolverParams(MappedParameterFrame):
     l_i: Parameter[float]
     """Normalised internal plasma inductance [dimensionless]."""
 
-    q_95: Parameter[float]
-    """Plasma safety factor [dimensionless]."""
-
     r_tf_inboard_out: Parameter[float]
     """Outboard Radius of the TF coil inboard leg tapered region [meter]."""
 
-    sigma_tf_case_max: Parameter[float]
-    """Maximum von Mises stress in the TF coil case nose [pascal]."""
-
-    sigma_tf_wp_max: Parameter[float]
-    """Maximum von Mises stress in the TF coil winding pack nose [pascal]."""
-
     T_e: Parameter[float]
     """Average plasma electron temperature [kiloelectron_volt]."""
 
-    tau_flattop: Parameter[float]
-    """Flat-top duration [second]."""
-
     tk_tf_outboard: Parameter[float]
     """TF coil outboard thickness [meter]."""
 

bluemira/codes/utilities.py

@@ -113,8 +113,14 @@ def create_mapping(
     ]:
         if puts is not None:
             for bm_key, (ec_key, unit) in puts.items():
+                if isinstance(ec_key, tuple):
+                    ec_in = ec_key[0]
+                    ec_out = ec_key[1]
+                else:
+                    ec_in = ec_out = ec_key
+
                 mappings[bm_key] = ParameterMapping(
-                    ec_key, send=sr["send"], recv=sr["recv"], unit=unit
+                    ec_in, ec_out, send=sr["send"], recv=sr["recv"], unit=unit
                 )
 
     return mappings

data/reactors/EU-DEMO/systems_code/mockPROCESS.json

@@ -5,7 +5,7 @@
   "H_star": 1.1,
   "I_p": 19.117,
   "P_brehms": 60.956,
-  "P_el_net_process": 500.0,
+  "P_el_net": 500.0,
   "P_fus": 1790.6,
   "P_fus_DD": 2.1206,
   "P_fus_DT": 1788.5,

eudemo/config/mockPROCESS.json

@@ -8,7 +8,6 @@
   "P_bd_in": 50.0,
   "P_brehms": 59.668,
   "P_el_net": 500,
-  "P_el_net_process": 500.0,
   "P_fus": 1994.6,
   "P_fus_DD": 2.4746,
   "P_fus_DT": 1992.1,