From afc95e2ee45e2c97246889631dfe948e8d277054 Mon Sep 17 00:00:00 2001 From: kripner-personal Date: Mon, 17 Apr 2023 10:55:24 +0200 Subject: [PATCH 01/10] First version of OMAS/IMAS import. --- freegs/omasio.py | 123 ++++++++++++++++++++++++++++++++++++++++++++ freegs/test_omas.py | 17 ++++++ requirements.txt | 1 + 3 files changed, 141 insertions(+) create mode 100644 freegs/omasio.py create mode 100644 freegs/test_omas.py diff --git a/freegs/omasio.py b/freegs/omasio.py new file mode 100644 index 0000000..a0b2898 --- /dev/null +++ b/freegs/omasio.py @@ -0,0 +1,123 @@ +# IDS schema view: https://gafusion.github.io/omas/schema.html + +import omas +import numpy as np +from .machine import ShapedCoil, FilamentCoil + +# OMAS coil types: +# multi-element coil is always translated to FilamentCoil (only support rectangular geometry) +# 'outline', 'rectangle' --> ShapeCoil +# oblique, arcs of circle, thick line and annulus are not supported at the moment +OMAS_COIL_TYPES = {1: 'outline', + 2: 'rectangle', + 3: 'oblique', + 4: 'arcs of circle', + 5: 'annulus', + 6: 'thick line'} + + +def _identify_name(ods: omas.ODS) -> str | None: + if "name" in ods: + return ods["name"] + elif "identifier" in ods: + return ods["identifier"] + else: + return None + + +def _identify_geometry_type(ods: omas.ODS) -> str | None: + """ Identify coil geometry type from OMAS data. + :param ods: pf_active.coil.element data + :return: coil geometry type + """ + geometry_type_idx = ods["geometry"]["geometry_type"] if "geometry" in ods else None + geometry_type = OMAS_COIL_TYPES[geometry_type_idx] if geometry_type_idx else None + + if not geometry_type: + if "outline" in ods["geometry"]: + geometry_type = "outline" + elif "rectangle" in ods["geometry"]: + geometry_type = "rectangle" + # The IDS definition of annulus is unclear to me. + # elif "annulus" in ods["geometry"]: + # geometry_type = "annulus" + else: + geometry_type = None + + return geometry_type + + +def _load_omas_coil(ods: omas.ODS): + # Identify coil name + coil_name = _identify_name(ods) + + # Read current + if "current" in ods: + if len(ods["current"]["data"]) > 1: + print( + f"Warning: multiple coil currents found. Using first one for time: " + f"{ods['pf_active.coil'][idx]['current']['time'][0]}") + + coil_current = ods["current"]["data"][0] + else: + coil_current = 0.0 + + # Multicoil or simple coil? + if len(ods["element"]) > 1: + r_filaments = ods["element"][:]["geometry"]["rectangle"]["r"] + z_filaments = ods["element"][:]["geometry"]["rectangle"]["z"] + turns = ods["element"][:]["turns_with_sign"] + if not np.all(np.abs(turns) == 1): + raise ValueError("Multicoil with non-unit turns is not supported yet.") + + # TODO: check if turns are interpreted correctly + coil = (coil_name, FilamentCoil(r_filaments, z_filaments, coil_current, turns=len(r_filaments))) + else: + print("Simple coil") + if not coil_name: + coil_name = _identify_name(ods["element"][0]) + + if not coil_name: + raise ValueError(f"Coil name not identified: \n {ods['pf_active.coil'][idx]}") + + geometry_type = _identify_geometry_type(ods["element"][0]) + + print(f"Coil name: {coil_name}") + + # Read turns: + if "turns" in ods["element"][0]: + turns = ods["element"][0]["turns_with_sign"] + else: + turns = 1 + + # Init FreeGS coil + if geometry_type == "outline": + outline_r = ods["element"][0]["geometry"]["outline"]["r"] + outline_z = ods["element"][0]["geometry"]["outline"]["z"] + coil = (coil_name, ShapedCoil(zip(outline_r, outline_z)), coil_current, turns) + elif geometry_type == "rectangle": + r_centre = ods["element"][0]["geometry"]["rectangle"]["r"] + z_centre = ods["element"][0]["geometry"]["rectangle"]["z"] + width = ods["element"][0]["geometry"]["rectangle"]["width"] + height = ods["element"][0]["geometry"]["rectangle"]["height"] + coil_r = [r_centre - width / 2, r_centre + width / 2, r_centre + width / 2, r_centre - width / 2] + coil_z = [z_centre - height / 2, z_centre - height / 2, z_centre + height / 2, z_centre + height / 2] + coil = (coil_name, ShapedCoil(zip(coil_r, coil_z)), coil_current, turns) + else: + raise ValueError(f"Coil geometry type {geometry_type} not supported yet.") + return coil + + +def _load_omas_coils(ods: omas.ODS): + coils = [_load_omas_coil(ods["pf_active.coil"][idx]) for idx in ods["pf_active.coil"]] + return coils + + +def main(): + ods = omas.ods_sample() + _load_omas_coils(ods) + + +if __name__ == '__main__': + np.linspace(10, 20) + main() diff --git a/freegs/test_omas.py b/freegs/test_omas.py new file mode 100644 index 0000000..8096ae0 --- /dev/null +++ b/freegs/test_omas.py @@ -0,0 +1,17 @@ +import omas + +from freegs.omasio import _load_omas_coils + + +# ODS feature with pf_active data +def coils_ods(): + ods = omas.ods_sample() + # TODO add multi-element coil + return ods + + +def test_omas_coils(): + ods = coils_ods() + + coils = _load_omas_coils(ods) + print(coils) diff --git a/requirements.txt b/requirements.txt index c79fce3..cb9dece 100644 --- a/requirements.txt +++ b/requirements.txt @@ -4,3 +4,4 @@ matplotlib>=2.0.0 h5py>=2.10.0 Shapely>=1.7.1 importlib-metadata<4.3,>=1.1.0 +omas>=0.56.0 From db937a11949dd217d39fd39886c0d75037677eab Mon Sep 17 00:00:00 2001 From: kripner-personal Date: Mon, 17 Apr 2023 11:09:21 +0200 Subject: [PATCH 02/10] Adjust test. Fix types. --- freegs/omasio.py | 27 +++++++++++---------------- freegs/test_omas.py | 7 ++++++- 2 files changed, 17 insertions(+), 17 deletions(-) diff --git a/freegs/omasio.py b/freegs/omasio.py index a0b2898..b941c50 100644 --- a/freegs/omasio.py +++ b/freegs/omasio.py @@ -1,7 +1,10 @@ # IDS schema view: https://gafusion.github.io/omas/schema.html +from typing import Tuple, Type, List import omas import numpy as np + +from freegs.coil import Coil from .machine import ShapedCoil, FilamentCoil # OMAS coil types: @@ -47,7 +50,9 @@ def _identify_geometry_type(ods: omas.ODS) -> str | None: return geometry_type -def _load_omas_coil(ods: omas.ODS): +def _load_omas_coil(ods: omas.ODS) -> Tuple[str, Type[Coil]]: + """ Load coil from OMAS data. + :param ods: 'pf_active.coil.:' data""" # Identify coil name coil_name = _identify_name(ods) @@ -56,7 +61,7 @@ def _load_omas_coil(ods: omas.ODS): if len(ods["current"]["data"]) > 1: print( f"Warning: multiple coil currents found. Using first one for time: " - f"{ods['pf_active.coil'][idx]['current']['time'][0]}") + f"{ods['current']['time'][0]}") coil_current = ods["current"]["data"][0] else: @@ -78,7 +83,7 @@ def _load_omas_coil(ods: omas.ODS): coil_name = _identify_name(ods["element"][0]) if not coil_name: - raise ValueError(f"Coil name not identified: \n {ods['pf_active.coil'][idx]}") + raise ValueError(f"Coil name not identified: \n {ods}") geometry_type = _identify_geometry_type(ods["element"][0]) @@ -94,7 +99,7 @@ def _load_omas_coil(ods: omas.ODS): if geometry_type == "outline": outline_r = ods["element"][0]["geometry"]["outline"]["r"] outline_z = ods["element"][0]["geometry"]["outline"]["z"] - coil = (coil_name, ShapedCoil(zip(outline_r, outline_z)), coil_current, turns) + coil = (coil_name, ShapedCoil(list(zip(outline_r, outline_z))), coil_current, turns) elif geometry_type == "rectangle": r_centre = ods["element"][0]["geometry"]["rectangle"]["r"] z_centre = ods["element"][0]["geometry"]["rectangle"]["z"] @@ -102,22 +107,12 @@ def _load_omas_coil(ods: omas.ODS): height = ods["element"][0]["geometry"]["rectangle"]["height"] coil_r = [r_centre - width / 2, r_centre + width / 2, r_centre + width / 2, r_centre - width / 2] coil_z = [z_centre - height / 2, z_centre - height / 2, z_centre + height / 2, z_centre + height / 2] - coil = (coil_name, ShapedCoil(zip(coil_r, coil_z)), coil_current, turns) + coil = (coil_name, ShapedCoil(list(zip(coil_r, coil_z))), coil_current, turns) else: raise ValueError(f"Coil geometry type {geometry_type} not supported yet.") return coil -def _load_omas_coils(ods: omas.ODS): +def _load_omas_coils(ods: omas.ODS) -> List[Tuple[str, Type[Coil]]]: coils = [_load_omas_coil(ods["pf_active.coil"][idx]) for idx in ods["pf_active.coil"]] return coils - - -def main(): - ods = omas.ods_sample() - _load_omas_coils(ods) - - -if __name__ == '__main__': - np.linspace(10, 20) - main() diff --git a/freegs/test_omas.py b/freegs/test_omas.py index 8096ae0..59428b3 100644 --- a/freegs/test_omas.py +++ b/freegs/test_omas.py @@ -14,4 +14,9 @@ def test_omas_coils(): ods = coils_ods() coils = _load_omas_coils(ods) - print(coils) + + coil_names = [coil[0] for coil in coils] + + assert "samp0" in coil_names + assert "samp1" in coil_names + assert "samp2" in coil_names From b0d44351018205434aee83d5763452c0bf09e860 Mon Sep 17 00:00:00 2001 From: kripner-personal Date: Tue, 18 Apr 2023 18:08:56 +0200 Subject: [PATCH 03/10] add dependency to pyproject.toml --- pyproject.toml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/pyproject.toml b/pyproject.toml index beff0ff..3698e6e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -32,6 +32,8 @@ dependencies = [ "Shapely>=1.7.1", "importlib-metadata<4.3,>=1.1.0", "freeqdsk>=0.1.0", + "omas>=0.56.0", + ] dynamic = ["version"] From e5b1dd911b21e23de12c5184e41c8cb18d021faa Mon Sep 17 00:00:00 2001 From: kripner-personal Date: Wed, 19 Apr 2023 10:12:50 +0200 Subject: [PATCH 04/10] fix typing and code compatibility --- freegs/omasio.py | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/freegs/omasio.py b/freegs/omasio.py index b941c50..5ff0ace 100644 --- a/freegs/omasio.py +++ b/freegs/omasio.py @@ -1,4 +1,6 @@ # IDS schema view: https://gafusion.github.io/omas/schema.html +from __future__ import annotations + from typing import Tuple, Type, List import omas @@ -50,7 +52,7 @@ def _identify_geometry_type(ods: omas.ODS) -> str | None: return geometry_type -def _load_omas_coil(ods: omas.ODS) -> Tuple[str, Type[Coil]]: +def _load_omas_coil(ods: omas.ODS) -> Tuple[str, Coil]: """ Load coil from OMAS data. :param ods: 'pf_active.coil.:' data""" # Identify coil name @@ -76,7 +78,7 @@ def _load_omas_coil(ods: omas.ODS) -> Tuple[str, Type[Coil]]: raise ValueError("Multicoil with non-unit turns is not supported yet.") # TODO: check if turns are interpreted correctly - coil = (coil_name, FilamentCoil(r_filaments, z_filaments, coil_current, turns=len(r_filaments))) + coil = FilamentCoil(r_filaments, z_filaments, coil_current, turns=len(r_filaments)) else: print("Simple coil") if not coil_name: @@ -99,7 +101,7 @@ def _load_omas_coil(ods: omas.ODS) -> Tuple[str, Type[Coil]]: if geometry_type == "outline": outline_r = ods["element"][0]["geometry"]["outline"]["r"] outline_z = ods["element"][0]["geometry"]["outline"]["z"] - coil = (coil_name, ShapedCoil(list(zip(outline_r, outline_z))), coil_current, turns) + coil = ShapedCoil(list(zip(outline_r, outline_z)), coil_current, turns) elif geometry_type == "rectangle": r_centre = ods["element"][0]["geometry"]["rectangle"]["r"] z_centre = ods["element"][0]["geometry"]["rectangle"]["z"] @@ -107,12 +109,17 @@ def _load_omas_coil(ods: omas.ODS) -> Tuple[str, Type[Coil]]: height = ods["element"][0]["geometry"]["rectangle"]["height"] coil_r = [r_centre - width / 2, r_centre + width / 2, r_centre + width / 2, r_centre - width / 2] coil_z = [z_centre - height / 2, z_centre - height / 2, z_centre + height / 2, z_centre + height / 2] - coil = (coil_name, ShapedCoil(list(zip(coil_r, coil_z))), coil_current, turns) + coil = ShapedCoil(list(zip(coil_r, coil_z)), coil_current, turns) else: raise ValueError(f"Coil geometry type {geometry_type} not supported yet.") - return coil + + if coil_name is None: + raise ValueError(f"Coil name not identified: \n {ods}") + + coil_tuple = (coil_name, coil) + return coil_tuple -def _load_omas_coils(ods: omas.ODS) -> List[Tuple[str, Type[Coil]]]: +def _load_omas_coils(ods: omas.ODS) -> List[Tuple[str, Coil]]: coils = [_load_omas_coil(ods["pf_active.coil"][idx]) for idx in ods["pf_active.coil"]] return coils From 1f5ea35771832cc5796a3c817ddeb1e7a67192c7 Mon Sep 17 00:00:00 2001 From: Lukas Kripner Date: Wed, 19 Apr 2023 13:49:41 +0200 Subject: [PATCH 05/10] Circuit loader. --- freegs/omasio.py | 109 ++++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 99 insertions(+), 10 deletions(-) diff --git a/freegs/omasio.py b/freegs/omasio.py index 5ff0ace..bd03d67 100644 --- a/freegs/omasio.py +++ b/freegs/omasio.py @@ -1,13 +1,13 @@ # IDS schema view: https://gafusion.github.io/omas/schema.html from __future__ import annotations -from typing import Tuple, Type, List +from typing import Tuple, Type, List, Dict import omas import numpy as np from freegs.coil import Coil -from .machine import ShapedCoil, FilamentCoil +from .machine import ShapedCoil, FilamentCoil, Circuit # OMAS coil types: # multi-element coil is always translated to FilamentCoil (only support rectangular geometry) @@ -22,6 +22,10 @@ def _identify_name(ods: omas.ODS) -> str | None: + """ Identify coil name from OMAS data. + Primary identifier is 'name', secondary is 'identifier'. + Note: Not sure what is an intended difference between 'name' and 'identifier'. + """ if "name" in ods: return ods["name"] elif "identifier" in ods: @@ -30,6 +34,16 @@ def _identify_name(ods: omas.ODS) -> str | None: return None +def _load_omas_power_supplies(ods: omas.ODS) -> List[Dict]: + """ Load power supplies names from OMAS data. + :param ods: 'pf_active.power_supply.:' data""" + # FreeGS does not have PowerSupply class, so we return data as list of dicts + # The voltages can be loaded in the same way as currents to support equilibrium evolution. + power_supplies_names = [{"name": _identify_name(ods[idx]), + "current": _load_omas_current(ods[idx])} for idx in ods] + return power_supplies_names + + def _identify_geometry_type(ods: omas.ODS) -> str | None: """ Identify coil geometry type from OMAS data. :param ods: pf_active.coil.element data @@ -52,22 +66,97 @@ def _identify_geometry_type(ods: omas.ODS) -> str | None: return geometry_type -def _load_omas_coil(ods: omas.ODS) -> Tuple[str, Coil]: - """ Load coil from OMAS data. - :param ods: 'pf_active.coil.:' data""" - # Identify coil name - coil_name = _identify_name(ods) +def _load_omas_current(ods: omas.ODS) -> float: + """ Load current from OMAS data. + :param ods: any IDS substructure which can contain current structure data""" # Read current if "current" in ods: if len(ods["current"]["data"]) > 1: print( - f"Warning: multiple coil currents found. Using first one for time: " + f"Warning: multiple circuit currents found. Using first one for time: " f"{ods['current']['time'][0]}") - coil_current = ods["current"]["data"][0] + circuit_current = ods["current"]["data"][0] else: - coil_current = 0.0 + circuit_current = 0.0 + + return circuit_current + + +def _circuit_connection_to_linear(circuit_structure: np.ndarray, n_supplies: int) -> Tuple[Tuple, Tuple]: + num_rows, num_cols = circuit_structure.shape + + supply_idx = set() + coil_idx = set() + + # Loop through each node in the circuit + for i in range(num_rows): + # Loop through each supply or coil side connected to the node + for j in range(num_cols): + if circuit_structure[i, j] == 1: + # Determine if the connection is to a supply or coil + if j < 2 * n_supplies: + index = j // 2 + supply_idx.add(index) + else: + index = (j - 2 * n_supplies) // 2 + coil_idx.add(index) + + return tuple(supply_idx), tuple(coil_idx) + + +def _load_omas_circuit(ods: omas.ODS, coils: List[Tuple[str, Coil]], power_supplies: List[Dict]) -> Tuple[str, Circuit]: + """ Load circuit from OMAS data. + :param ods: 'pf_active.circuit.:' data""" + # Identify circuit name + + # IDS circuit description can be found here: https://gafusion.github.io/omas/schema/schema_pf%20active.html. + + # Get linear circuit (coil and supply) structure + supply_idx, coil_idx = _circuit_connection_to_linear(ods["connections"], len(power_supplies)) + + if len(supply_idx) == 0 or len(coil_idx) == 0: + raise ValueError(f"Invalid circuit structure. No supplies or coils found for circuit {ods}.") + + if len(supply_idx) > 1: + print(f"Warning: multiple supplies found for circuit {ods}.") + + # Construct circuit name. First from circuit name, then from supply name, then from coil names. + circuit_name = _identify_name(ods) + if not circuit_name: + if power_supplies[supply_idx[0]]["name"]: + supply_names = [power_supplies[supply_idx[idx]]["name"] for idx in supply_idx if + power_supplies[supply_idx[idx]]["name"]] + circuit_name = "+".join(supply_names) + elif coils[coil_idx[0]][0]: + coil_names = [coil[0] for coil in coils if coil[0]] + circuit_name = "+".join(coil_names) + else: + raise ValueError(f"Unable to identify circuit name for circuit {ods}.") + + circuit_current = _load_omas_current(ods) + + # Init FreeGS circuit + # TODO: Recognize correctly the multiplier for the circuit current + circuit = Circuit([(coils[idx][0], coils[idx[1]], 1.0) for idx in coil_idx], circuit_current) + return circuit_name, circuit + + +def _load_omas_circuits(ods: omas.ODS) -> List[Tuple[str, Circuit]]: + coils = _load_omas_coils(ods) + power_supplies = _load_omas_power_supplies(ods) + return [_load_omas_circuit(ods["pf_active.circuit"][idx], coils, power_supplies) for idx in ods["pf_active.circuit"]] + + +def _load_omas_coil(ods: omas.ODS) -> Tuple[str, Coil]: + """ Load coil from OMAS data. + :param ods: 'pf_active.coil.:' data""" + # Identify coil name + coil_name = _identify_name(ods) + + # Read current + coil_current = _load_omas_current(ods) # Multicoil or simple coil? if len(ods["element"]) > 1: From 06517e87b807c206be60aabe9b8b3a95633fb003 Mon Sep 17 00:00:00 2001 From: Lukas Kripner Date: Wed, 19 Apr 2023 14:54:43 +0200 Subject: [PATCH 06/10] Load limiter and complete machine. --- freegs/omasio.py | 49 +++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 46 insertions(+), 3 deletions(-) diff --git a/freegs/omasio.py b/freegs/omasio.py index bd03d67..0691277 100644 --- a/freegs/omasio.py +++ b/freegs/omasio.py @@ -7,7 +7,7 @@ import numpy as np from freegs.coil import Coil -from .machine import ShapedCoil, FilamentCoil, Circuit +from .machine import ShapedCoil, FilamentCoil, Circuit, Wall, Machine # OMAS coil types: # multi-element coil is always translated to FilamentCoil (only support rectangular geometry) @@ -143,10 +143,13 @@ def _load_omas_circuit(ods: omas.ODS, coils: List[Tuple[str, Coil]], power_suppl return circuit_name, circuit -def _load_omas_circuits(ods: omas.ODS) -> List[Tuple[str, Circuit]]: +def _load_omas_circuits(ods: omas.ODS) -> List[Tuple[str, Circuit]] | None: coils = _load_omas_coils(ods) power_supplies = _load_omas_power_supplies(ods) - return [_load_omas_circuit(ods["pf_active.circuit"][idx], coils, power_supplies) for idx in ods["pf_active.circuit"]] + if "pf_active.circuit" not in ods: + return None + return [_load_omas_circuit(ods["pf_active.circuit"][idx], coils, power_supplies) for idx in + ods["pf_active.circuit"]] def _load_omas_coil(ods: omas.ODS) -> Tuple[str, Coil]: @@ -212,3 +215,43 @@ def _load_omas_coil(ods: omas.ODS) -> Tuple[str, Coil]: def _load_omas_coils(ods: omas.ODS) -> List[Tuple[str, Coil]]: coils = [_load_omas_coil(ods["pf_active.coil"][idx]) for idx in ods["pf_active.coil"]] return coils + + +def _load_omas_limiter(ods: omas.ODS): + # Look for the limiter in the ODS. + + # Try to find limiter description in the OMAS data. + # If not found, raise. + if "limiter" not in ods["wall.description_2d.:"]: + return None + + # Identify the first occurrence of limiter id ods + limiter_id = next((idx for idx in ods["wall.description_2d"] if "limiter" in ods["wall.description_2d"][idx])) + limiter_ods = ods["wall.description_2d"][limiter_id] + + # Currently is supported only continuous limiter + if limiter_ods["type.index"] != 0: + raise ValueError("Only continuous limiter is supported yet.") + + # Create FreeGS wall object + limiter = Wall(limiter_ods["unit.0.outline.r"], limiter_ods["unit.0.outline.z"]) + + return limiter + + +def load_omas_machine(ods: omas.ODS): + """ Load machine from OMAS data. + :param ods: OMAS data""" + + # Load circuits + coils = _load_omas_circuits(ods) + if not coils: + coils = _load_omas_coils(ods) + + # Load limiter + limiter = _load_omas_limiter(ods) + + # Load machine + machine = Machine(coils, limiter) + + return machine From bb620f4f111385011d19bb3ae8f05170dd65900e Mon Sep 17 00:00:00 2001 From: Lukas Kripner Date: Thu, 20 Apr 2023 12:02:12 +0200 Subject: [PATCH 07/10] Add first wall plotting. --- freegs/machine.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/freegs/machine.py b/freegs/machine.py index 3062634..9e34efe 100644 --- a/freegs/machine.py +++ b/freegs/machine.py @@ -858,6 +858,10 @@ def plot(self, axis=None, show=True): """ for label, coil in self.coils: axis = coil.plot(axis=axis, show=False) + + if self.wall is not None: + axis.plot(self.wall.R, self.wall.Z, 'k-') + if show: import matplotlib.pyplot as plt From 1a8bb9b23335cdefba96946b3d00e799ca0752b5 Mon Sep 17 00:00:00 2001 From: Lukas Kripner Date: Thu, 20 Apr 2023 12:02:53 +0200 Subject: [PATCH 08/10] Add test coil and test reconstruction. --- freegs/omasio.py | 6 +---- freegs/test_omas.py | 64 ++++++++++++++++++++++++++++++++++++++++----- 2 files changed, 58 insertions(+), 12 deletions(-) diff --git a/freegs/omasio.py b/freegs/omasio.py index 0691277..d77509f 100644 --- a/freegs/omasio.py +++ b/freegs/omasio.py @@ -76,7 +76,6 @@ def _load_omas_current(ods: omas.ODS) -> float: print( f"Warning: multiple circuit currents found. Using first one for time: " f"{ods['current']['time'][0]}") - circuit_current = ods["current"]["data"][0] else: circuit_current = 0.0 @@ -172,7 +171,6 @@ def _load_omas_coil(ods: omas.ODS) -> Tuple[str, Coil]: # TODO: check if turns are interpreted correctly coil = FilamentCoil(r_filaments, z_filaments, coil_current, turns=len(r_filaments)) else: - print("Simple coil") if not coil_name: coil_name = _identify_name(ods["element"][0]) @@ -181,8 +179,6 @@ def _load_omas_coil(ods: omas.ODS) -> Tuple[str, Coil]: geometry_type = _identify_geometry_type(ods["element"][0]) - print(f"Coil name: {coil_name}") - # Read turns: if "turns" in ods["element"][0]: turns = ods["element"][0]["turns_with_sign"] @@ -227,7 +223,7 @@ def _load_omas_limiter(ods: omas.ODS): # Identify the first occurrence of limiter id ods limiter_id = next((idx for idx in ods["wall.description_2d"] if "limiter" in ods["wall.description_2d"][idx])) - limiter_ods = ods["wall.description_2d"][limiter_id] + limiter_ods = ods["wall.description_2d"][limiter_id]["limiter"] # Currently is supported only continuous limiter if limiter_ods["type.index"] != 0: diff --git a/freegs/test_omas.py b/freegs/test_omas.py index 59428b3..7b10b2a 100644 --- a/freegs/test_omas.py +++ b/freegs/test_omas.py @@ -1,22 +1,72 @@ +import matplotlib.pyplot as plt +import numpy as np import omas +import pytest -from freegs.omasio import _load_omas_coils +import freegs +from freegs.omasio import _load_omas_coils, load_omas_machine # ODS feature with pf_active data -def coils_ods(): +@pytest.fixture +def example_ods(): ods = omas.ods_sample() - # TODO add multi-element coil - return ods + ods["pf_active.coil.+.element"][0]["geometry"]["rectangle"]["r"] = 1.2 + ods["pf_active.coil.-1.element"][0]["geometry"]["rectangle"]["z"] = -1.5 + ods["pf_active.coil.-1.element"][0]["geometry"]["rectangle"]["width"] = 0.1 + ods["pf_active.coil.-1.element"][0]["geometry"]["rectangle"]["height"] = 0.1 + # set geometry_type = "rectangle" + ods["pf_active.coil.-1.element"][0]["geometry"]["geometry_type"] = 2 + ods["pf_active.coil.-1.name"] = "test_coil" + + return ods -def test_omas_coils(): - ods = coils_ods() - coils = _load_omas_coils(ods) +def test_omas_coils(example_ods): + coils = _load_omas_coils(example_ods) coil_names = [coil[0] for coil in coils] assert "samp0" in coil_names assert "samp1" in coil_names assert "samp2" in coil_names + + +def test_omas_machine(example_ods): + machine = load_omas_machine(example_ods) + + +# Run only on request, since it takes a while to run +@pytest.mark.skip(reason="Run only on request") +def test_omas_test_reconstruction(example_ods): + machine = load_omas_machine(example_ods) + + eq = freegs.Equilibrium(machine, + Rmin=0.5, Rmax=3.0, + Zmin=-1.6, Zmax=1.5, + nx=65, ny=65) + + ip = example_ods["equilibrium"]["time_slice"][0]["global_quantities"]["ip"] + p_ax = example_ods["equilibrium"]["time_slice"][0]["profiles_1d"]["pressure"][0] + btor = example_ods["equilibrium"]["time_slice"][0]["global_quantities"]["magnetic_axis"]["b_field_tor"] + rax = example_ods["equilibrium"]["time_slice"][0]["global_quantities"]["magnetic_axis"]["r"] + f0 = btor * rax + + profiles = freegs.jtor.ConstrainPaxisIp(eq, p_ax, ip, f0) + + # Try some circular equilibrium + r_lcfs = example_ods["equilibrium"]["time_slice"][0]["boundary"]["outline"]["r"] + z_lcfs = example_ods["equilibrium"]["time_slice"][0]["boundary"]["outline"]["z"] + + isoflux = [(r_lcfs[0], z_lcfs[0], r, z) for r, z in zip(r_lcfs[1:], z_lcfs[1:])] + constraints = freegs.control.constrain(isoflux=isoflux) + + freegs.solve(eq, profiles, constraints) + + ax = plt.gca() + ax.set_aspect("equal") + eq.plot(axis=ax, show=False) + machine.plot(axis=ax, show=False) + constraints.plot(axis=ax, show=False) + plt.show() From 085c3fd2a46a435abc3b3ded1250b6079f79afcd Mon Sep 17 00:00:00 2001 From: Lukas Kripner Date: Thu, 20 Apr 2023 17:30:33 +0200 Subject: [PATCH 09/10] Add circuit test --- freegs/omasio.py | 31 ++++++++++-------- freegs/test_omas.py | 80 ++++++++++++++++++++++++++++++++++++++------- 2 files changed, 86 insertions(+), 25 deletions(-) diff --git a/freegs/omasio.py b/freegs/omasio.py index d77509f..37fc85d 100644 --- a/freegs/omasio.py +++ b/freegs/omasio.py @@ -34,9 +34,9 @@ def _identify_name(ods: omas.ODS) -> str | None: return None -def _load_omas_power_supplies(ods: omas.ODS) -> List[Dict]: +def _load_omas_supplies(ods: omas.ODS) -> List[Dict]: """ Load power supplies names from OMAS data. - :param ods: 'pf_active.power_supply.:' data""" + :param ods: 'pf_active.supply.:' data""" # FreeGS does not have PowerSupply class, so we return data as list of dicts # The voltages can be loaded in the same way as currents to support equilibrium evolution. power_supplies_names = [{"name": _identify_name(ods[idx]), @@ -84,6 +84,7 @@ def _load_omas_current(ods: omas.ODS) -> float: def _circuit_connection_to_linear(circuit_structure: np.ndarray, n_supplies: int) -> Tuple[Tuple, Tuple]: + # TODO This function does not support non-linear circuits and connections which generate opposite currents. num_rows, num_cols = circuit_structure.shape supply_idx = set() @@ -92,6 +93,7 @@ def _circuit_connection_to_linear(circuit_structure: np.ndarray, n_supplies: int # Loop through each node in the circuit for i in range(num_rows): # Loop through each supply or coil side connected to the node + node_count = 0 for j in range(num_cols): if circuit_structure[i, j] == 1: # Determine if the connection is to a supply or coil @@ -101,6 +103,9 @@ def _circuit_connection_to_linear(circuit_structure: np.ndarray, n_supplies: int else: index = (j - 2 * n_supplies) // 2 coil_idx.add(index) + node_count += 1 + if node_count > 2: + raise ValueError(f"Non-linear circuit structure. Node {i} has more than 2 connections.") return tuple(supply_idx), tuple(coil_idx) @@ -138,13 +143,13 @@ def _load_omas_circuit(ods: omas.ODS, coils: List[Tuple[str, Coil]], power_suppl # Init FreeGS circuit # TODO: Recognize correctly the multiplier for the circuit current - circuit = Circuit([(coils[idx][0], coils[idx[1]], 1.0) for idx in coil_idx], circuit_current) + circuit = Circuit([(coils[idx][0], coils[idx][1], 1.0) for idx in coil_idx], circuit_current) return circuit_name, circuit -def _load_omas_circuits(ods: omas.ODS) -> List[Tuple[str, Circuit]] | None: - coils = _load_omas_coils(ods) - power_supplies = _load_omas_power_supplies(ods) +def load_omas_circuits(ods: omas.ODS) -> List[Tuple[str, Circuit]] | None: + coils = load_omas_coils(ods) + power_supplies = _load_omas_supplies(ods["pf_active.supply"]) if "pf_active.circuit" not in ods: return None return [_load_omas_circuit(ods["pf_active.circuit"][idx], coils, power_supplies) for idx in @@ -162,14 +167,14 @@ def _load_omas_coil(ods: omas.ODS) -> Tuple[str, Coil]: # Multicoil or simple coil? if len(ods["element"]) > 1: - r_filaments = ods["element"][:]["geometry"]["rectangle"]["r"] - z_filaments = ods["element"][:]["geometry"]["rectangle"]["z"] - turns = ods["element"][:]["turns_with_sign"] + r_filaments = ods["element.:.geometry.rectangle.r"] + z_filaments = ods["element.:.geometry.rectangle.z"] + turns = ods["element.:.turns_with_sign"] if not np.all(np.abs(turns) == 1): raise ValueError("Multicoil with non-unit turns is not supported yet.") # TODO: check if turns are interpreted correctly - coil = FilamentCoil(r_filaments, z_filaments, coil_current, turns=len(r_filaments)) + coil = FilamentCoil(r_filaments, z_filaments, current=coil_current, turns=len(r_filaments)) else: if not coil_name: coil_name = _identify_name(ods["element"][0]) @@ -208,7 +213,7 @@ def _load_omas_coil(ods: omas.ODS) -> Tuple[str, Coil]: return coil_tuple -def _load_omas_coils(ods: omas.ODS) -> List[Tuple[str, Coil]]: +def load_omas_coils(ods: omas.ODS) -> List[Tuple[str, Coil]]: coils = [_load_omas_coil(ods["pf_active.coil"][idx]) for idx in ods["pf_active.coil"]] return coils @@ -240,9 +245,9 @@ def load_omas_machine(ods: omas.ODS): :param ods: OMAS data""" # Load circuits - coils = _load_omas_circuits(ods) + coils = load_omas_circuits(ods) if not coils: - coils = _load_omas_coils(ods) + coils = load_omas_coils(ods) # Load limiter limiter = _load_omas_limiter(ods) diff --git a/freegs/test_omas.py b/freegs/test_omas.py index 7b10b2a..a231bf8 100644 --- a/freegs/test_omas.py +++ b/freegs/test_omas.py @@ -4,7 +4,7 @@ import pytest import freegs -from freegs.omasio import _load_omas_coils, load_omas_machine +from freegs.omasio import load_omas_coils, load_omas_machine, load_omas_circuits # ODS feature with pf_active data @@ -12,19 +12,65 @@ def example_ods(): ods = omas.ods_sample() - ods["pf_active.coil.+.element"][0]["geometry"]["rectangle"]["r"] = 1.2 - ods["pf_active.coil.-1.element"][0]["geometry"]["rectangle"]["z"] = -1.5 - ods["pf_active.coil.-1.element"][0]["geometry"]["rectangle"]["width"] = 0.1 - ods["pf_active.coil.-1.element"][0]["geometry"]["rectangle"]["height"] = 0.1 - # set geometry_type = "rectangle" - ods["pf_active.coil.-1.element"][0]["geometry"]["geometry_type"] = 2 - ods["pf_active.coil.-1.name"] = "test_coil" + r_coil = 1.3 + z_coil = -1.5 + dr_coil = 0.1 + dz_coil = 0.1 + n_elements = 9 + + # Add a coil with 9 filaments: + ods["pf_active.coil.+.name"] = "multicoil" + + for i in range(n_elements): + icol = i % 3 + irow = i // 3 + rc = r_coil - 3 * dr_coil + irow * dr_coil + zc = z_coil - 3 * dz_coil + icol * dz_coil + + ods["pf_active.coil.-1.element.+.geometry"]["rectangle"]["r"] = rc + ods["pf_active.coil.-1.element.-1.geometry"]["rectangle"]["z"] = zc + ods["pf_active.coil.-1.element.-1.geometry"]["rectangle"]["width"] = 0.1 + ods["pf_active.coil.-1.element.-1.geometry"]["rectangle"]["height"] = 0.1 + ods["pf_active.coil.-1.element.-1.geometry"]["geometry_type"] = 2 + + # For the consistency with `example_circuit_ods` fixture: + ods["pf_active.coil.-1.element.0.identifier"] = "sampl_circuit" + + return ods + + +@pytest.fixture +def example_circuit_ods(example_ods): + ods = example_ods.copy() + + n_coils = len(ods["pf_active.coil"]) + + # Generate Power Supplies: + for coil_idx in ods["pf_active.coil"]: + coil = ods["pf_active.coil"][coil_idx] + ods["pf_active.supply.+.name"] = "PS_" + coil["element.0.identifier"] + + for coil_idx in ods["pf_active.coil"]: + coil = ods["pf_active.coil"][coil_idx] + # Add a circuit: + ods["pf_active.circuit.+.name"] = "C_" + coil["element.0.identifier"] + + # Connect both poles in a loop + connections = np.zeros((2, 4 * n_coils), dtype=int) + # Connection on two poles + connections[0, 2 * coil_idx] = 1 + connections[0, 2 * coil_idx + 2 * n_coils] = 1 + # Connection of the second two poles + connections[1, 2 * coil_idx + 1] = 1 + connections[1, 2 * coil_idx + 2 * n_coils + 1] = 1 + + ods["pf_active.circuit.-1.connections"] = connections return ods def test_omas_coils(example_ods): - coils = _load_omas_coils(example_ods) + coils = load_omas_coils(example_ods) coil_names = [coil[0] for coil in coils] @@ -33,18 +79,28 @@ def test_omas_coils(example_ods): assert "samp2" in coil_names +def test_omas_circuits(example_circuit_ods): + circuits = load_omas_circuits(example_circuit_ods) + + circuit_names = [circuit[0] for circuit in circuits] + + assert "C_samp0" in circuit_names + assert "C_samp1" in circuit_names + assert "C_samp2" in circuit_names + + def test_omas_machine(example_ods): machine = load_omas_machine(example_ods) # Run only on request, since it takes a while to run -@pytest.mark.skip(reason="Run only on request") +# @pytest.mark.skip(reason="Run only on request") def test_omas_test_reconstruction(example_ods): machine = load_omas_machine(example_ods) eq = freegs.Equilibrium(machine, - Rmin=0.5, Rmax=3.0, - Zmin=-1.6, Zmax=1.5, + Rmin=0.6, Rmax=2.8, + Zmin=-2.0, Zmax=1.8, nx=65, ny=65) ip = example_ods["equilibrium"]["time_slice"][0]["global_quantities"]["ip"] From 72da4acd7c7e6a378f6585ce66b3b7c5a8b5257e Mon Sep 17 00:00:00 2001 From: Lukas Kripner Date: Thu, 20 Apr 2023 17:34:57 +0200 Subject: [PATCH 10/10] Add circuit test --- freegs/test_omas.py | 44 +++++++++++++++++++++++++++----------------- 1 file changed, 27 insertions(+), 17 deletions(-) diff --git a/freegs/test_omas.py b/freegs/test_omas.py index a231bf8..3fde7c4 100644 --- a/freegs/test_omas.py +++ b/freegs/test_omas.py @@ -93,36 +93,46 @@ def test_omas_machine(example_ods): machine = load_omas_machine(example_ods) -# Run only on request, since it takes a while to run -# @pytest.mark.skip(reason="Run only on request") -def test_omas_test_reconstruction(example_ods): - machine = load_omas_machine(example_ods) - +def _test_freegs_inverse(machine, ods): + # Inverse problem eq = freegs.Equilibrium(machine, Rmin=0.6, Rmax=2.8, Zmin=-2.0, Zmax=1.8, nx=65, ny=65) - ip = example_ods["equilibrium"]["time_slice"][0]["global_quantities"]["ip"] - p_ax = example_ods["equilibrium"]["time_slice"][0]["profiles_1d"]["pressure"][0] - btor = example_ods["equilibrium"]["time_slice"][0]["global_quantities"]["magnetic_axis"]["b_field_tor"] - rax = example_ods["equilibrium"]["time_slice"][0]["global_quantities"]["magnetic_axis"]["r"] + ip = ods["equilibrium"]["time_slice"][0]["global_quantities"]["ip"] + p_ax = ods["equilibrium"]["time_slice"][0]["profiles_1d"]["pressure"][0] + btor = ods["equilibrium"]["time_slice"][0]["global_quantities"]["magnetic_axis"]["b_field_tor"] + rax = ods["equilibrium"]["time_slice"][0]["global_quantities"]["magnetic_axis"]["r"] f0 = btor * rax profiles = freegs.jtor.ConstrainPaxisIp(eq, p_ax, ip, f0) # Try some circular equilibrium - r_lcfs = example_ods["equilibrium"]["time_slice"][0]["boundary"]["outline"]["r"] - z_lcfs = example_ods["equilibrium"]["time_slice"][0]["boundary"]["outline"]["z"] + r_lcfs = ods["equilibrium"]["time_slice"][0]["boundary"]["outline"]["r"] + z_lcfs = ods["equilibrium"]["time_slice"][0]["boundary"]["outline"]["z"] isoflux = [(r_lcfs[0], z_lcfs[0], r, z) for r, z in zip(r_lcfs[1:], z_lcfs[1:])] constraints = freegs.control.constrain(isoflux=isoflux) freegs.solve(eq, profiles, constraints) - ax = plt.gca() - ax.set_aspect("equal") - eq.plot(axis=ax, show=False) - machine.plot(axis=ax, show=False) - constraints.plot(axis=ax, show=False) - plt.show() + # ax = plt.gca() + # ax.set_aspect("equal") + # eq.plot(axis=ax, show=False) + # machine.plot(axis=ax, show=False) + # constraints.plot(axis=ax, show=False) + # plt.show() + + +# Run only on request, since it takes a while to run +@pytest.mark.skip(reason="Run only on request") +def test_omas_test_reconstruction(example_ods): + machine = load_omas_machine(example_ods) + _test_freegs_inverse(machine, example_ods) + + +@pytest.mark.skip(reason="Run only on request") +def test_omas_w_circuits_test_reconstruction(example_circuit_ods): + machine = load_omas_machine(example_circuit_ods) + _test_freegs_inverse(machine, example_circuit_ods)