diff --git a/README.md b/README.md index 3c36865..c9f76be 100644 --- a/README.md +++ b/README.md @@ -68,12 +68,14 @@ Templates are ready-made configurations for common IDSs, and for the results of some kinds of codes. Start from one with `idscope -T NAME SOURCE ...` or File → New from template; `idscope --list-templates` prints their names. Templates for the results of a kind of code also go by the names of such codes: `-T nubdec` (NuBDeC) -is `nbi`, and `-T sonamu` and `-T freegs` (SONAMU, FreeGS) are `free_boundary`. +is `nbi`, `-T sonamu` and `-T freegs` (SONAMU, FreeGS) are `free_boundary`, and +`-T autofrd` (autofrd, which draws a machine) is `device`. | template | shows | |---|---| | `core_profiles` | Ip, central Te/Ti/ne, loop voltage, li and beta over time; Te/Ti, ne/ni, q, current densities and Zeff profiles; Te and ne against (t, x) | | `core_sources` | for every source (beams, waves, ohmic, ...): profiles of electron and ion heating, current density, particles and torque density; power, driven current and torque over time; the power at the selected time, as bars | +| `device` (`autofrd`) | the machine in the poloidal plane: the vacuum vessel (its layers, between their inner and outer contours), the wall (including parts only at some toroidal angles, such as ports), the passive conductors, the PF coils and the central solenoid, and the conductor of the TF coils; the currents of the coils and passive conductors, and the vacuum toroidal field, when the source has them | | `equilibrium` | profiles: q, pressure, current densities (j_phi, j_parallel), dp/dpsi, FF', elongation and triangularities | | `free_boundary` (`sonamu`, `freegs`) | the poloidal plane: flux surfaces, plasma boundary, X-points and those the solver was asked to hold, magnetic axis, coils and wall; the coil currents, as bars; p', FF', q and pressure against psi_n | | `nbi` (`nubdec`) | what the beams give the plasma: profiles of electron and ion heating, driven current and torque densities; fast ion density and pressures; where the ions of the beams are born (markers) and where the neutrons are (a map per reaction), in the poloidal plane; the neutron rates over time, D-D beam-thermal and thermal apart | @@ -87,7 +89,9 @@ voltages of a pf_active written by an equilibrium code, which only has currents, the totals of a core_sources written by a beam deposition code, which only has profiles. The output of a beam code has no equilibrium and no wall: to draw the births over the flux surfaces, give the data entry the code was run on as a second -source. +source. A machine description has the geometry and no currents; given the outputs +of several runs of autofrd, the `device` template draws each machine with a line +style of its own. A template saved with File → Save as becomes a configuration that can be edited. ![equilibrium template on the ITER 134173/106 training data](docs/template-equilibrium.png) @@ -100,6 +104,8 @@ A template saved with File → Save as becomes a configuration that can be edite ![free_boundary template comparing two synthetic sources](docs/template-free_boundary.png) +![device template comparing two synthetic sources](docs/template-device.png) + ## Using the scope | jScope | idscope | @@ -265,7 +271,12 @@ distribution_sources/source[:]/ggd/particles[0]/values a map per source, on a maps: the 2D geometry of coils and passive loops (outline, rectangle, parallelogram, arcs of circle, annulus or thick line), `r`/`z` outlines, `r`/`z` points, and the markers of a particle distribution, as dots (the columns - of their `positions` that their `coordinate_identifier` says are R and Z). The IDS + of their `positions` that their `coordinate_identifier` says are R and Z). Also + the layers of a vessel (`wall/description_2d[0]/vessel/unit[:]/annular`), filled + between their inner and outer contours, or around their centreline as thick as + their `thickness`; and conductors made of line segments, arcs of circle and full + circles (`tf/coil[:]/conductor[:]/elements`), drawn at their distance from the + axis and their height, since their points are in three dimensions. The IDS browser lists them as *shape*. - The `values` of a GGD quantity on a structured grid (grid identifier `structured_spaces`: two spaces of one dimension, e.g. R and Z) are a map, on the diff --git a/docs/template-device.png b/docs/template-device.png new file mode 100644 index 0000000..fecd625 Binary files /dev/null and b/docs/template-device.png differ diff --git a/src/idscope/data/extract.py b/src/idscope/data/extract.py index 291ae19..c5e872a 100644 --- a/src/idscope/data/extract.py +++ b/src/idscope/data/extract.py @@ -13,6 +13,8 @@ ``map`` (z vs x1, x2) ``equilibrium/time_slice/profiles_2d[0]/psi`` ``shape`` (in R, Z) ``pf_active/coil[:]/element[:]/geometry``, ``wall/description_2d[0]/limiter/unit[:]/outline``, + ``wall/description_2d[0]/vessel/unit[:]/annular``, + ``tf/coil[:]/conductor[:]/elements``, ``equilibrium/time_slice/boundary/x_point[:]`` ======================== ========================================================= @@ -47,6 +49,9 @@ GEOMETRY, OUTLINE, POINT, MARKERS = "geometry", "outline", "point", "markers" """Kinds of shapes: closed cross-sections, curves, points and the markers of a particle distribution (e.g. where the ions of a neutral beam are born).""" +ANNULAR, CONDUCTOR = "annular", "conductor" +"""Also kinds of shapes: a layer of a vessel between two contours, and a conductor +made of line segments and arcs of circle (e.g. a toroidal field coil).""" _R_INDEX, _Z_INDEX = 4, 3 """Indices of R and Z in the Data Dictionary's coordinate identifiers.""" @@ -369,8 +374,9 @@ def is_time_aos(meta) -> bool: def shape_kind(meta) -> Optional[str]: """The kind of shape a structure (or the elements of an array of structures) - describes: a Data Dictionary 2D geometry, an (r, z) curve, an (r, z) point or - the markers of a particle distribution.""" + describes: a Data Dictionary 2D geometry, an (r, z) curve, an (r, z) point, + the markers of a particle distribution, a layer of a vessel or the elements of + a conductor.""" if meta.data_type not in (IDSDataType.STRUCTURE, IDSDataType.STRUCT_ARRAY): return None children = {child.name: child for child in meta} @@ -378,6 +384,10 @@ def shape_kind(meta) -> Optional[str]: return GEOMETRY if {"positions", "coordinate_identifier"} <= set(children): return MARKERS + if {"outline_inner", "outline_outer"} <= set(children): + return ANNULAR + if {"types", "start_points", "end_points"} <= set(children): + return CONDUCTOR r, z = children.get("r"), children.get("z") if r is None or z is None: return None @@ -673,10 +683,16 @@ def walk(node, pos, key, labels): walk(ids, 0, (), ()) members = [] seen: Dict[str, int] = {} + steps = [plan.steps[pos].name for pos in plan.channels] for key in sorted(found): labels = found[key] - # Levels with a single element (e.g. coils of one element) add nothing - shown = [name for name, count in labels if count > 1] or [labels[-1][0]] + # Levels with a single element (e.g. coils of one element) add nothing; when + # all have one, the last level with a name tells the member (the coil "TF" + # rather than its one conductor, "conductor[0]") + named = [name for (name, _), step in zip(labels, steps) + if not name.startswith(f"{step}[")] + shown = [name for name, count in labels if count > 1] or \ + (named[-1:] or [labels[-1][0]]) name = "/".join(shown) # Elements can share a name, e.g. two core_sources both "nbi": nbi, nbi#2 seen[name] = seen.get(name, 0) + 1 @@ -1127,6 +1143,10 @@ def read_shape(node, shape: str) -> Optional[np.ndarray]: return geometry_outline(node) if shape == MARKERS: return marker_positions(node) + if shape == ANNULAR: + return annular_outline(node) + if shape == CONDUCTOR: + return conductor_line(node) r, z = _leaf_array(node.r), _leaf_array(node.z) if shape == POINT: if not (np.isfinite(r) and np.isfinite(z)): @@ -1157,14 +1177,17 @@ def marker_positions(markers) -> Optional[np.ndarray]: def _is_closed(outline, points: np.ndarray) -> bool: - """Whether a curve is closed: as its ``closed`` flag says (e.g. limiter units + """Whether a curve is closed: as its ``closed`` flag says, or that of the + structure around it (e.g. the contours of a vessel layer and the limiter units of DD 3), else when its ends are about one segment apart.""" - try: - from imas.util import get_parent + flag = getattr(outline, "closed", None) + if not (isinstance(flag, IDSPrimitive) and flag.has_value): + try: + from imas.util import get_parent - flag = getattr(get_parent(outline), "closed", None) - except Exception: # pragma: no cover - depends on the imas-python version - flag = None + flag = getattr(get_parent(outline), "closed", None) + except Exception: # pragma: no cover - depends on the imas-python version + flag = None if isinstance(flag, IDSPrimitive) and flag.has_value: return int(flag.value) == 1 finite = points[np.all(np.isfinite(points), axis=1)] @@ -1285,6 +1308,197 @@ def geometry_outline(geometry) -> Optional[np.ndarray]: return None +def _curve(node) -> Optional[Tuple[np.ndarray, bool]]: + """The finite (R, Z) points of an (r, z) curve, and whether it is closed; None + when it is not filled.""" + r, z = _leaf_array(node.r).ravel(), _leaf_array(node.z).ravel() + if r.size == 0 or r.shape != z.shape: + return None + points = np.column_stack([r, z]) + points = points[np.all(np.isfinite(points), axis=1)] + if len(points) < 2: + return None + return points, _is_closed(node, points) + + +def _widened(points: np.ndarray, thickness: np.ndarray, closed: bool): + """The two sides of a band around a curve: ``thickness`` is its width across + each segment (from point i to i + 1), or one width for all. None when it does + not fit the curve.""" + segments = (np.roll(points, -1, axis=0) - points) if closed else np.diff(points, axis=0) + thickness = thickness[np.isfinite(thickness)] + if thickness.size == 1: + thickness = np.full(len(segments), float(thickness[0])) + if thickness.size < len(segments): + return None + thickness = thickness[: len(segments)] + lengths = np.hypot(*segments.T) + lengths[lengths == 0] = 1.0 + normals = np.column_stack([-segments[:, 1], segments[:, 0]]) / lengths[:, None] + # At a point, between the segments on either side of it (the ends of an open + # curve have one), mitred at corners + if closed: + before, after = np.roll(normals, 1, axis=0), normals + widths = (np.roll(thickness, 1) + thickness) / 2.0 + else: + before, after = np.vstack([normals[:1], normals]), np.vstack([normals, normals[-1:]]) + widths = (np.concatenate([thickness[:1], thickness]) + + np.concatenate([thickness, thickness[-1:]])) / 2.0 + total = before + after + norms = np.hypot(*total.T) + norms[norms == 0] = 1.0 + mitre = np.minimum(2.0 / norms, 4.0) # 1 / cos of half the turn at the point + offset = total / norms[:, None] * (widths / 2.0 * mitre)[:, None] + return points + offset, points - offset + + +def annular_outline(annular) -> Optional[np.ndarray]: + """The contours of a layer of a vessel (``annular``): its inner and outer + outlines, else its centreline widened by its thickness. Closed contours are + separated by a row of NaN, so that the ring between them can be filled; open + ones make the band between them, as one closed outline. None when it is not + filled.""" + inner, outer = _curve(annular.outline_inner), _curve(annular.outline_outer) + if inner is not None and outer is not None: + (inner, closed), (outer, _) = inner, outer + else: + centre = _curve(annular.centreline) + if centre is None: + return None + points, closed = centre + sides = _widened(points, _leaf_array(annular.thickness).ravel(), closed) + if sides is None: + return None + inner, outer = sides + if closed: + return np.vstack([_closed(inner), np.full((1, 2), np.nan), _closed(outer)]) + # The band: along one side and back along the other + if np.hypot(*(inner[-1] - outer[-1])) < np.hypot(*(inner[-1] - outer[0])): + outer = outer[::-1] + return _closed(np.vstack([inner, outer])) + + +_TURN_POINTS = 64 +"""Points on a full circle; arcs have as many for the angle they turn by.""" + + +def conductor_line(elements) -> Optional[np.ndarray]: + """(R, Z) along the elements of a conductor, e.g. of a toroidal field coil: + line segments, arcs of circle through an intermediate point, and full circles + about a centre. The elements are in three dimensions (r, phi, z); in the + poloidal plane they are drawn at their distance from the axis and their height. + Rows of NaN separate elements that do not meet. None when it is not filled.""" + types = _leaf_array(elements.types).ravel() + count = types.size + if count == 0: + return None + + def points(node) -> np.ndarray: + """(x, y, z) of the point of every element; NaN where it is not given.""" + r, z = _leaf_array(node.r).ravel(), _leaf_array(node.z).ravel() + if r.size != count or z.size != count: + return np.full((count, 3), np.nan) + phi = _leaf_array(node.phi).ravel() if hasattr(node, "phi") else np.empty(0) + phi = np.where(np.isfinite(phi), phi, 0.0) if phi.size == count \ + else np.zeros(count) + return np.column_stack([r * np.cos(phi), r * np.sin(phi), z]) + + start, middle, end, centre = ( + points(getattr(elements, name)) + for name in ("start_points", "intermediate_points", "end_points", "centres") + ) + parts: List[np.ndarray] = [] + last = None + for i, kind in enumerate(types): + if kind == 1: + drawn = _segment(start[i], end[i]) + elif kind == 2: + drawn = _arc_through(start[i], middle[i], end[i], centre[i]) + elif kind == 3: + drawn = _full_circle(centre[i], start[i], middle[i]) + else: + drawn = None + if drawn is None: + continue + if last is not None and not np.allclose(drawn[0], last, rtol=0.0, atol=1e-6): + parts.append(np.full((1, 3), np.nan)) + parts.append(drawn) + last = drawn[-1] + if not parts: + return None + xyz = np.vstack(parts) + return np.column_stack([np.hypot(xyz[:, 0], xyz[:, 1]), xyz[:, 2]]) + + +def _finite(*points) -> bool: + return all(bool(np.all(np.isfinite(point))) for point in points) + + +def _segment(a: np.ndarray, b: np.ndarray) -> Optional[np.ndarray]: + """Points along a line segment: its ends, and points between them when it + crosses toroidal angles (it then curves in R).""" + if not _finite(a, b): + return None + across = abs(a[0] * b[1] - a[1] * b[0]) > 1e-9 * (1.0 + a[:2] @ a[:2] + b[:2] @ b[:2]) + t = np.linspace(0.0, 1.0, 17 if across else 2)[:, None] + return a + t * (b - a) + + +def _on_circle(centre, u, v, radius: float, sweep: float) -> np.ndarray: + n = max(2, int(np.ceil(abs(sweep) / (2.0 * np.pi) * _TURN_POINTS))) + 1 + angles = np.linspace(0.0, sweep, n)[:, None] + return centre + radius * (np.cos(angles) * u + np.sin(angles) * v) + + +def _arc_through(a, middle, b, centre) -> Optional[np.ndarray]: + """Points along the arc of circle from ``a`` through ``middle`` to ``b``; without + the intermediate point, about the ``centre`` the shorter way.""" + if not _finite(a, b): + return None + if not _finite(middle): + if not _finite(centre): + return _segment(a, b) + mean = (a - centre) + (b - centre) + if np.linalg.norm(mean) == 0: # half a circle, either way + return None + middle = centre + np.linalg.norm(a - centre) * mean / np.linalg.norm(mean) + ab, ac = middle - a, b - a + normal = np.cross(ab, ac) + area = normal @ normal + if area <= 1e-24 * (ab @ ab) * (ac @ ac) or area == 0: # on a line + return _segment(a, b) + # The centre of the circle through the three points + circle = a + (np.cross(normal, ab) * (ac @ ac) + np.cross(ac, normal) * (ab @ ab)) \ + / (2.0 * area) + radius = float(np.linalg.norm(a - circle)) + u = (a - circle) / radius + v = np.cross(normal / np.sqrt(area), u) + # About the normal of a, middle, b, going from a the positive way passes the + # intermediate point before b + sweep = np.arctan2((b - circle) @ v, (b - circle) @ u) % (2.0 * np.pi) + return _on_circle(circle, u, v, radius, sweep) + + +def _full_circle(centre, a, middle) -> Optional[np.ndarray]: + """Points around a full circle about ``centre`` through ``a``, in the plane of + ``middle`` (else in the poloidal plane of ``a``).""" + if not _finite(centre, a): + return None + radius = float(np.linalg.norm(a - centre)) + if radius == 0: + return None + u = (a - centre) / radius + normal = np.cross(a - centre, middle - centre) if _finite(middle) else np.zeros(3) + if np.linalg.norm(normal) <= 1e-12 * radius**2: + phi = np.arctan2(a[1], a[0]) + normal = np.array([-np.sin(phi), np.cos(phi), 0.0]) # toroidal direction + normal = normal / np.linalg.norm(normal) + v = np.cross(normal, u) + if np.linalg.norm(v) == 0: + return None + return _on_circle(centre, u, v / np.linalg.norm(v), radius, 2.0 * np.pi) + + def apply_x(data: SignalData, x_data: SignalData) -> SignalData: """Use the values of ``x_data`` as x coordinate of the profiles in ``data``. diff --git a/src/idscope/demo.py b/src/idscope/demo.py index 52fc6a2..4b36785 100644 --- a/src/idscope/demo.py +++ b/src/idscope/demo.py @@ -472,8 +472,15 @@ def process(source, name: str, energy: Optional[str] = None) -> None: return ds +def _d_shape(a: float, b: float, theta: np.ndarray): + """(R, Z) of a D-shaped contour of the demo machine: half-width ``a``, half-height + ``b``, centred on R = 1.8 m.""" + return 1.8 + a * np.cos(theta + 0.35 * np.sin(theta)), b * np.sin(theta) + + def wall(factory): - """wall: a D-shaped first wall around the plasma, as one limiter unit.""" + """wall: a D-shaped first wall around the plasma, as one limiter unit, and the + vacuum vessel around it, as one annular unit (its inner and outer contours).""" from imas.ids_defs import IDS_TIME_MODE_INDEPENDENT w = factory.wall() @@ -483,11 +490,107 @@ def wall(factory): unit = w.description_2d[0].limiter.unit[0] unit.name = "first wall" theta = np.linspace(0.0, 2.0 * np.pi, 121)[:-1] # closed: the first point is not repeated - unit.outline.r = 1.8 + 0.72 * np.cos(theta + 0.35 * np.sin(theta)) - unit.outline.z = 1.15 * np.sin(theta) + unit.outline.r, unit.outline.z = _d_shape(0.72, 1.15, theta) + w.description_2d[0].vessel.unit.resize(1) + vessel = w.description_2d[0].vessel.unit[0] + vessel.name = "vacuum vessel" + for contour, (a, b) in ((vessel.annular.outline_inner, (0.8, 1.25)), + (vessel.annular.outline_outer, (0.85, 1.32))): + contour.r, contour.z = _d_shape(a, b, theta) + if hasattr(contour, "closed"): # Data Dictionary 3 + contour.closed = 1 return w +def pf_passive(factory, scale: float = 1.0, n_time: int = N_TIME): + """pf_passive: two stabilizing plates inside the vessel, above and below the + midplane on the outboard side, with the currents induced in them: by the ramp-up + of the plasma current, and, in opposite directions, by the plasma moving up and + down.""" + from imas.ids_defs import IDS_TIME_MODE_HOMOGENEOUS + + pf = factory.pf_passive() + pf.ids_properties.homogeneous_time = IDS_TIME_MODE_HOMOGENEOUS + t = np.linspace(0.0, 10.0, n_time) + pf.time = t + pf.loop.resize(2) + for loop, name, sign in zip(pf.loop, ("upper plate", "lower plate"), (1.0, -1.0)): + loop.name = name + loop.element.resize(1) + geometry = loop.element[0].geometry + geometry.geometry_type = 1 # an outline: the plate follows the wall + theta = sign * np.linspace(0.3, 0.8, 11) + inner, outer = _d_shape(0.745, 1.19, theta), _d_shape(0.775, 1.22, theta) + geometry.outline.r = np.concatenate([inner[0], outer[0][::-1]]) + geometry.outline.z = np.concatenate([inner[1], outer[1][::-1]]) + loop.element[0].turns_with_sign = 1.0 + loop.current = scale * (-25.0e3 * np.exp(-(((t - 1.0) / 0.6) ** 2)) + + sign * 3.0e3 * np.sin(3.7 * t)) + return pf + + +TF_COILS = 16 +_S = np.sqrt(0.5) +TF_CONDUCTOR = [ + # name, type (1: line, 2: arc), start, intermediate, end, centre (R, Z) [m] + ("inner leg", 1, (0.9, -1.5), None, (0.9, 1.5), None), + ("upper inner arc", 2, (0.9, 1.5), (1.33 - 0.43 * _S, 1.5 + 0.43 * _S), (1.33, 1.93), + (1.33, 1.5)), + ("top", 1, (1.33, 1.93), None, (2.3, 1.93), None), + ("upper outer arc", 2, (2.3, 1.93), (2.3 + 0.83 * _S, 1.1 + 0.83 * _S), (3.13, 1.1), + (2.3, 1.1)), + ("outer leg", 1, (3.13, 1.1), None, (3.13, -1.1), None), + ("lower outer arc", 2, (3.13, -1.1), (2.3 + 0.83 * _S, -1.1 - 0.83 * _S), (2.3, -1.93), + (2.3, -1.1)), + ("bottom", 1, (2.3, -1.93), None, (1.33, -1.93), None), + ("lower inner arc", 2, (1.33, -1.93), (1.33 - 0.43 * _S, -1.5 - 0.43 * _S), (0.9, -1.5), + (1.33, -1.5)), +] +"""The center line of the conductor of a toroidal field coil of the demo machine: +a D of straight legs and arcs, around the vessel and the PF coils.""" + + +def tf(factory, scale: float = 1.0, n_time: int = N_TIME): + """tf: the toroidal field coils (one standing for all, :data:`TF_COILS` of + them), their current and the vacuum field they make.""" + from imas.ids_defs import EMPTY_FLOAT, IDS_TIME_MODE_HOMOGENEOUS + + coils = factory.tf() + coils.ids_properties.homogeneous_time = IDS_TIME_MODE_HOMOGENEOUS + t = np.linspace(0.0, 10.0, n_time) + coils.time = t + coils.r0 = 1.8 + coils.is_periodic = 1 + coils.coils_n = TF_COILS + # R B_phi = mu0 N I / (2 pi), for N I of all the coils together + current = scale * 11.25e3 * np.ones(len(t)) # per turn, of 100 turns + field = 2.0e-7 * TF_COILS * 100 * current + vacuum = (coils.b_field_phi_vacuum_r if hasattr(coils, "b_field_phi_vacuum_r") + else coils.b_field_tor_vacuum_r) # Data Dictionary 4, 3 + vacuum.data = field + coils.coil.resize(1) + coil = coils.coil[0] + coil.name = "TF" + coil.turns = 100.0 + coil.current.data = current + coil.conductor.resize(1) + elements = coil.conductor[0].elements + if hasattr(elements, "names"): # not in every version of the Data Dictionary + elements.names = [name for name, *_ in TF_CONDUCTOR] + elements.types = np.array([kind for _, kind, *_ in TF_CONDUCTOR], dtype=np.int32) + + def column(points, k): + return np.array([EMPTY_FLOAT if p is None else p[k] for p in points]) + + for name, index in (("start_points", 2), ("intermediate_points", 3), + ("end_points", 4), ("centres", 5)): + points = [row[index] for row in TF_CONDUCTOR] + node = getattr(elements, name) + node.r, node.z = column(points, 0), column(points, 1) + node.phi = np.array([EMPTY_FLOAT if p is None else 0.0 for p in points]) + return coils + + def summary(factory, scale: float = 1.0, n_time: int = N_TIME, time: Optional[float] = None): """summary: the plasma current and the neutron rates (``time``: write this @@ -526,8 +629,9 @@ def write_entry(uri: str, scale: float = 1.0, dd_version: Optional[str] = None, for ids in (core_profiles(factory, scale, rho_tor_norm=rho_tor_norm), equilibrium(factory, scale, rho_tor_norm=rho_tor_norm), magnetics(factory, scale), summary(factory, scale), - pf_active(factory, scale), wall(factory), - core_sources(factory, scale), distribution_sources(factory, scale)): + pf_active(factory, scale), wall(factory), pf_passive(factory, scale), + tf(factory, scale), core_sources(factory, scale), + distribution_sources(factory, scale)): entry.put(ids) diff --git a/src/idscope/gui/panel.py b/src/idscope/gui/panel.py index 602a69f..c63f7e8 100644 --- a/src/idscope/gui/panel.py +++ b/src/idscope/gui/panel.py @@ -12,7 +12,7 @@ from idscope.config import PanelConfig, SignalConfig from idscope.data import MAP, PROFILE, SHAPE, TRACE, Channels, MapSeries, ProfileSeries from idscope.data import Shapes, SignalData, parse_alternatives -from idscope.data.extract import GEOMETRY, MARKERS, POINT +from idscope.data.extract import ANNULAR, GEOMETRY, MARKERS, POINT from idscope.data.sources import DataSource from idscope.errors import IdscopeError, NoDataError from idscope.gui import style @@ -715,7 +715,7 @@ def _make_shape_item(self, data: Shapes, pen, source_index: int): elif data.shape == MARKERS: # often thousands: small dots, no outline item = pg.PlotDataItem(pen=None, symbol="o", symbolSize=3, symbolPen=None, symbolBrush=pg.mkBrush(pen.color())) - elif data.shape == GEOMETRY: + elif data.shape in (GEOMETRY, ANNULAR): # filled (a vessel layer: its ring) fill = QtGui.QColor(pen.color()) fill.setAlpha(60) item = pg.PlotDataItem(pen=pen, connect="finite", fillLevel="enclosed", diff --git a/src/idscope/templates/__init__.py b/src/idscope/templates/__init__.py index 67ce72e..3d3b266 100644 --- a/src/idscope/templates/__init__.py +++ b/src/idscope/templates/__init__.py @@ -20,6 +20,7 @@ "nubdec": "nbi", "sonamu": "free_boundary", "freegs": "free_boundary", + "autofrd": "device", } """Other names of templates: codes whose results they show (e.g. in tbd-fusion).""" diff --git a/src/idscope/templates/device.yaml b/src/idscope/templates/device.yaml new file mode 100644 index 0000000..65daaf8 --- /dev/null +++ b/src/idscope/templates/device.yaml @@ -0,0 +1,62 @@ +# idscope template: a machine: its wall and vacuum vessel, its coils and passive +# conductors (autofrd in tbd-fusion writes these: "-T autofrd" is this template). +# +# idscope -T device SOURCE [SOURCE ...] +# +# The poloidal cross-section shows the vacuum vessel (the layers between the inner and +# outer contours of its annular units, or its block elements), the plasma-facing +# components (the limiter units of the wall; one that is only at some toroidal angles, +# such as a port, is drawn like the others), the passive conductors (pf_passive), the +# PF coils and the central solenoid (pf_active), and the center line of the conductor +# of the toroidal field coils (tf). The currents and the vacuum field of the TF coils +# are shown when the source has them: a machine description alone has the geometry, +# and the other panels then note that these are not filled. Given several sources, +# e.g. the runs of autofrd for two placements of the coils, the line style tells them +# apart. +title: device +columns: + - width: 1.5 + panels: + - title: Poloidal cross-section + kind: map + equal_aspect: true + x_label: R [m] + y_label: Z [m] + signals: + - path: wall/description_2d[0]/vessel/unit[:]/annular | wall/description_2d[0]/vessel/unit[:]/element[:]/outline + label: vessel + color: "#a0a0a0" + - path: wall/description_2d[0]/limiter/unit[:]/outline + label: wall + color: "#4c4c4c" + - path: pf_passive/loop[:]/element[:]/geometry + label: passive + color: "#8c564b" + - path: pf_active/coil[:]/element[:]/geometry | pf_active/coil[:]/geometry + label: coils + color: "#1f77b4" + - path: tf/coil[:]/conductor[:]/elements + label: TF coil + color: "#d62728" + - panels: + - title: Coil currents + kind: profile + signals: + - path: pf_active/coil[:]/current/data + scale: 1.0e-3 + y_label: I [kA] + - title: Passive currents + kind: profile + signals: + - path: pf_passive/loop[:]/current + scale: 1.0e-3 + y_label: I [kA] + - title: TF coil current + signals: + - path: tf/coil[:]/current/data + scale: 1.0e-3 + y_label: I [kA] + - title: Vacuum toroidal field + signals: + - path: tf/b_field_phi_vacuum_r/data | tf/b_field_tor_vacuum_r/data + y_label: R B_phi [T m] diff --git a/tests/test_browse.py b/tests/test_browse.py index 0dbb456..80df3f4 100644 --- a/tests/test_browse.py +++ b/tests/test_browse.py @@ -39,3 +39,18 @@ def test_shape_items(source_a): assert find_item(eq, "time_slice/global_quantities/magnetic_axis").kind == "shape" wall = find_item(list_signals(source_a, "wall"), "description_2d/limiter/unit/outline") assert wall.signal_path == "wall/description_2d[:]/limiter/unit[:]/outline" + + +def test_machine_shape_items(source_a): + """A layer of the vessel and the conductor of a coil are shapes as a whole, not + their contours or points one by one.""" + wall = list_signals(source_a, "wall") + vessel = find_item(wall, "description_2d/vessel/unit/annular") + assert vessel.kind == "shape" + assert vessel.signal_path == "wall/description_2d[:]/vessel/unit[:]/annular" + assert find_item(wall, "description_2d/vessel/unit/annular/outline_inner") is None + tf = list_signals(source_a, "tf") + conductor = find_item(tf, "coil/conductor/elements") + assert conductor.kind == "shape" and not conductor.time_dependent + assert conductor.signal_path == "tf/coil[:]/conductor[:]/elements" + assert find_item(tf, "coil/conductor/elements/start_points") is None diff --git a/tests/test_cli.py b/tests/test_cli.py index 7dac4c3..5b34e15 100644 --- a/tests/test_cli.py +++ b/tests/test_cli.py @@ -21,9 +21,11 @@ def test_list_templates(capsys): assert main(["--list-templates"]) == 0 lines = capsys.readouterr().out.splitlines() assert [line.split()[0] for line in lines] == [ - "core_profiles", "core_sources", "equilibrium", "free_boundary", "nbi", "pf_active" + "core_profiles", "core_sources", "device", "equilibrium", "free_boundary", "nbi", + "pf_active", ] # and the codes whose results they show + assert "device also: autofrd" in lines assert "free_boundary also: freegs, sonamu" in lines assert "nbi also: nubdec" in lines diff --git a/tests/test_extract.py b/tests/test_extract.py index 7c4c074..ed777b7 100644 --- a/tests/test_extract.py +++ b/tests/test_extract.py @@ -9,7 +9,14 @@ Trace, parse_signal_path, ) -from idscope.data.extract import geometry_outline, marker_positions, normalize, plan_path +from idscope.data.extract import ( + annular_outline, + conductor_line, + geometry_outline, + marker_positions, + normalize, + plan_path, +) from idscope.data.sources import DataSource, parse_source from idscope.errors import NoDataError, PathError @@ -454,6 +461,134 @@ def test_geometry_types(): assert np.all(geometry_outline(geometry)[1:16, 1] > 0) +def test_machine_shapes(source_a): + """The vessel between its two contours, the passive loops and the conductor of + the toroidal field coils, as the demo machine has them.""" + from idscope.demo import TF_CONDUCTOR + + vessel = get(source_a, "wall/description_2d[0]/vessel/unit[:]/annular") + assert vessel.shape == "annular" and vessel.names == ["vacuum vessel"] + (ring,), _ = vessel.at(None) + gap = np.flatnonzero(np.isnan(ring[:, 0])) + assert len(gap) == 1 # between the inner and the outer contour, both closed + inner, outer = ring[: gap[0]], ring[gap[0] + 1:] + np.testing.assert_allclose([inner[0], outer[0]], [inner[-1], outer[-1]]) + assert inner[:, 0].min() == pytest.approx(1.0) and outer[:, 0].max() == pytest.approx(2.65) + plates = get(source_a, "pf_passive/loop[:]/element[:]/geometry") + assert plates.shape == "geometry" and plates.names == ["upper plate", "lower plate"] + currents = get(source_a, "pf_passive/loop[:]/current") + assert isinstance(currents, Channels) and currents.names == plates.names + tf = get(source_a, "tf/coil[:]/conductor[:]/elements") + assert tf.shape == "conductor" and tf.names == ["TF"] # the coil, not its conductor + (line,), _ = tf.at(None) + assert not np.isnan(line).any() # the elements meet, all the way round + np.testing.assert_allclose(line[0], line[-1]) + np.testing.assert_allclose([line[:, 0].min(), line[:, 0].max(), line[:, 1].max()], + [0.9, 3.13, 1.93]) + # the arcs are on their circles + for _, kind, start, _, end, centre in TF_CONDUCTOR: + if kind != 2: + continue + radius = np.hypot(*np.subtract(start, centre)) + on_arc = np.isclose(np.hypot(line[:, 0] - centre[0], line[:, 1] - centre[1]), radius) + between = ((line[:, 0] - min(start[0], end[0]) >= -1e-9) + & (line[:, 0] - max(start[0], end[0]) <= 1e-9) + & (line[:, 1] - min(start[1], end[1]) >= -1e-9) + & (line[:, 1] - max(start[1], end[1]) <= 1e-9)) + assert on_arc[between].sum() >= 16 + field = get(source_a, "tf/b_field_phi_vacuum_r/data") + np.testing.assert_allclose(field.values, 3.6) + + +def test_conductor_elements(): + import imas + + tf = imas.IDSFactory().tf() + tf.coil.resize(1) + tf.coil[0].conductor.resize(1) + elements = tf.coil[0].conductor[0].elements + assert conductor_line(elements) is None # nothing filled + nan = np.nan + + def fill(types, start, middle, end, centre, phi=None): + if hasattr(elements, "names"): + elements.names = [f"element {i}" for i in range(len(types))] + elements.types = np.array(types, dtype=np.int32) + for node, points in ((elements.start_points, start), + (elements.intermediate_points, middle), + (elements.end_points, end), (elements.centres, centre)): + node.r = [p[0] for p in points] + node.z = [p[1] for p in points] + node.phi = phi if phi is not None else [0.0] * len(points) + + # a line, then an arc the long way round, through its intermediate point, and a + # full circle about its centre that starts elsewhere + s = np.sqrt(0.5) + fill([1, 2, 3], start=[(1, 0), (2, 0), (3, 3)], middle=[(nan, nan), (1 - s, -s), (nan, nan)], + end=[(2, 0), (1, 1), (nan, nan)], centre=[(nan, nan), (1, 0), (2.5, 3)]) + line = conductor_line(elements) + gaps = np.flatnonzero(np.isnan(line[:, 0])) + assert len(gaps) == 1 # before the circle, which does not start where the arc ends + np.testing.assert_allclose(line[0], [1, 0]) + arc = line[1: gaps[0]] + np.testing.assert_allclose(np.hypot(arc[:, 0] - 1, arc[:, 1]), 1.0) + assert arc[:, 1].min() == pytest.approx(-1.0) # three quarters, past the bottom + np.testing.assert_allclose(arc[-1], [1, 1]) + circle = line[gaps[0] + 1:] + np.testing.assert_allclose(np.hypot(circle[:, 0] - 2.5, circle[:, 1] - 3), 0.5) + np.testing.assert_allclose(circle[0], circle[-1]) + # an arc without its intermediate point goes the shorter way about its centre; + # types that are not known are left out + fill([2, 7], start=[(2, 0), (0, 0)], middle=[(nan, nan)] * 2, end=[(1, 1), (1, 1)], + centre=[(1, 0), (nan, nan)]) + line = conductor_line(elements) + assert not np.isnan(line).any() and line[:, 1].min() >= -1e-12 + # in three dimensions: a ring around the axis is one point in (R, Z); a straight + # line between two toroidal angles comes closer to the axis in between + fill([3], start=[(2, 1)], middle=[(2, 1)], end=[(nan, nan)], centre=[(0, 1)], + phi=[0.0]) + elements.intermediate_points.phi = [np.pi / 2] + np.testing.assert_allclose(conductor_line(elements), [[2, 1]] * 65, atol=1e-12) + fill([1], start=[(2, 0)], middle=[(nan, nan)], end=[(2, 0)], centre=[(nan, nan)]) + elements.end_points.phi = [np.pi / 2] + chord = conductor_line(elements) + assert len(chord) == 17 and chord[:, 0].min() == pytest.approx(np.sqrt(2)) + + +def test_annular_layers(): + import imas + + wall = imas.IDSFactory().wall() + wall.description_2d.resize(1) + wall.description_2d[0].vessel.unit.resize(1) + annular = wall.description_2d[0].vessel.unit[0].annular + assert annular_outline(annular) is None # nothing filled + # a centreline and a thickness per segment: a closed square, 0.2 thick + annular.centreline.r = [0.0, 1.0, 1.0, 0.0] + annular.centreline.z = [0.0, 0.0, 1.0, 1.0] + annular.thickness = [0.2] * 4 + ring = annular_outline(annular) + gap = np.flatnonzero(np.isnan(ring[:, 0]))[0] + np.testing.assert_allclose(ring[:gap], [[0.1, 0.1], [0.9, 0.1], [0.9, 0.9], [0.1, 0.9], + [0.1, 0.1]]) + np.testing.assert_allclose(ring[gap + 1:][:4], [[-0.1, -0.1], [1.1, -0.1], [1.1, 1.1], + [-0.1, 1.1]]) + annular.thickness = [0.2] # one for all + np.testing.assert_allclose(annular_outline(annular), ring) + annular.thickness = [0.2, 0.3] # too few for the segments + assert annular_outline(annular) is None + # open contours: the band between them, one way along one and back the other + z = np.linspace(-1.0, 1.0, 5) + annular.outline_inner.r, annular.outline_inner.z = np.full(5, 1.0), z + annular.outline_outer.r, annular.outline_outer.z = np.full(5, 1.2), z[::-1] + band = annular_outline(annular) + np.testing.assert_allclose(band[:5], np.column_stack([np.full(5, 1.0), z])) + np.testing.assert_allclose(band[5:], [[1.2, 1.0], [1.2, 0.5], [1.2, 0.0], [1.2, -0.5], + [1.2, -1.0], [1.0, -1.0]]) + annular.outline_outer.z = z # the same way along both: still the band + np.testing.assert_allclose(annular_outline(annular), band) + + @pytest.fixture(scope="module") def two_beams(tmp_path_factory): """Two core sources with the same identifier, on grids of different lengths.""" diff --git a/tests/test_templates.py b/tests/test_templates.py index f0fa2ac..37cba74 100644 --- a/tests/test_templates.py +++ b/tests/test_templates.py @@ -8,7 +8,8 @@ def test_template_names(): assert template_names() == [ - "core_profiles", "core_sources", "equilibrium", "free_boundary", "nbi", "pf_active" + "core_profiles", "core_sources", "device", "equilibrium", "free_boundary", "nbi", + "pf_active", ] @@ -17,6 +18,7 @@ def test_template_aliases(): assert set(ALIASES.values()) <= set(template_names()) assert aliases_of("nbi") == ["nubdec"] assert aliases_of("free_boundary") == ["freegs", "sonamu"] + assert aliases_of("device") == ["autofrd"] assert aliases_of("equilibrium") == [] assert load_template("nubdec").title == "nbi" assert load_template("SONAMU").title == "free_boundary" # in any case @@ -445,3 +447,122 @@ def test_free_boundary_template_on_demo_data(template_window, data_dir): assert x_points.data.path.text == "equilibrium/time_slice/contour_tree/node[1:]" np.testing.assert_allclose(np.ravel(x_points.item.getData()), [1.8 + 0.03 * np.sin(5.0) - 0.15, -0.8]) + + +@pytest.fixture(scope="module") +def machine_description(tmp_path_factory): + """What autofrd writes (through tbd-fusion) for a reactor it drew, in Data + Dictionary 3.39.0: the solenoid and the PF coils as rectangles of one turn and + no current, a part of the TF coil as a passive loop, the vessel as an annular + unit and the blanket and a port as limiter units, and the TF coil as the center + line of its conductor. No time base anywhere.""" + import imas + from imas.ids_defs import EMPTY_FLOAT + + path = tmp_path_factory.mktemp("autofrd") / "reactor.nc" + factory = imas.IDSFactory("3.39.0") + + def rectangle(element, r, z, width, height): + element.turns_with_sign = 1.0 + element.geometry.geometry_type = 2 + box = element.geometry.rectangle + box.r, box.z, box.width, box.height = r, z, width, height + + pf = factory.pf_active() + pf.ids_properties.homogeneous_time = 1 + pf.time = [0.0] + boxes = [(f"OH {k + 1}", 0.5, -1.2 + 0.8 * k + 0.4, 0.2, 0.8) for k in range(3)] + boxes += [("PF 1", 1.2, 1.9, 0.3, 0.3), ("PF 2", 2.8, 1.0, 0.3, 0.3), + ("PF 3", 2.8, -1.0, 0.3, 0.3), ("PF 4", 1.2, -1.9, 0.3, 0.3)] + pf.coil.resize(len(boxes)) + for coil, (name, *box) in zip(pf.coil, boxes): + coil.name = name + coil.element.resize(1) + rectangle(coil.element[0], *box) + + passive = factory.pf_passive() + passive.ids_properties.homogeneous_time = 1 + passive.loop.resize(1) + passive.loop[0].name = "TF_2" + passive.loop[0].element.resize(1) + rectangle(passive.loop[0].element[0], 0.7, 0.0, 0.05, 2.4) + + wall = factory.wall() + wall.ids_properties.homogeneous_time = 2 + wall.description_2d.resize(1) + drawing = wall.description_2d[0] + theta = np.linspace(0.0, 2.0 * np.pi, 73)[:-1] + drawing.vessel.unit.resize(1) + drawing.vessel.unit[0].name = "VV" + for contour, (a, b) in ((drawing.vessel.unit[0].annular.outline_inner, (0.9, 1.4)), + (drawing.vessel.unit[0].annular.outline_outer, (1.0, 1.5))): + contour.r, contour.z = 1.8 + a * np.cos(theta), b * np.sin(theta) + drawing.limiter.unit.resize(2) + blanket, port = drawing.limiter.unit + blanket.name, blanket.closed = "Blanket", 1 + blanket.outline.r, blanket.outline.z = 1.8 + 0.8 * np.cos(theta), 1.3 * np.sin(theta) + port.name, port.closed = "Port_E_array", 1 + port.outline.r, port.outline.z = [2.7, 3.3, 3.3, 2.7], [-0.2, -0.2, 0.2, 0.2] + port.phi_extensions = [[0.0, np.pi / 2], [0.3, 0.3]] + + tf = factory.tf() + tf.ids_properties.homogeneous_time = 2 + tf.is_periodic, tf.coils_n = 1, 16 + tf.coil.resize(1) + tf.coil[0].name = tf.coil[0].identifier = "TF" + tf.coil[0].conductor.resize(1) + elements = tf.coil[0].conductor[0].elements + # a D: the inner leg, and an arc of radius 1.4 about (0.75, 0) for the rest, + # through the outboard midplane + elements.names = ["inner leg", "outer arc"] + elements.types = np.array([1, 2], dtype=np.int32) + empty = float(EMPTY_FLOAT) # what autofrd writes where a point does not apply + for node, points in ( + (elements.start_points, [(0.75, -1.4), (0.75, 1.4)]), + (elements.intermediate_points, [(empty, empty), (0.75 + 1.4, 0.0)]), + (elements.end_points, [(0.75, 1.4), (0.75, -1.4)]), + (elements.centres, [(empty, empty), (0.75, 0.0)])): + node.r = [p[0] for p in points] + node.z = [p[1] for p in points] + node.phi = [empty if p[0] == empty else 0.0 for p in points] + with imas.DBEntry(str(path), "w", dd_version="3.39.0") as entry: + for ids in (pf, passive, wall, tf): + entry.put(ids) + return path + + +def test_device_template_on_a_machine_description(template_window, machine_description): + window = template_window("autofrd", [machine_description]) + section, coils, passive, tf_current, field = window.panels + labels = [label.text for _, label in section.legend.items] + assert labels == ["vessel", "wall", "passive", "coils", "TF coil"] + shapes = {c.label: c for c in section.curves} + assert shapes["wall"].data.names == ["Blanket", "Port_E_array"] + assert shapes["coils"].data.names[:4] == ["OH 1", "OH 2", "OH 3", "PF 1"] + assert shapes["passive"].data.names == ["TF_2"] + x, y = shapes["TF coil"].item.getData() + assert shapes["TF coil"].data.names == ["TF"] + np.testing.assert_allclose([np.nanmin(x), np.nanmax(x), np.nanmin(y), np.nanmax(y)], + [0.75, 2.15, -1.4, 1.4]) + # the ring of the vessel is filled + assert shapes["vessel"].item.opts["fillLevel"] == "enclosed" + # The machine description has no currents: noted, not errors + assert not any(panel.errors for panel in window.panels) + assert set(errors_by_panel(window)) == {"Coil currents", "Passive currents", + "TF coil current", "Vacuum toroidal field"} + assert "not filled" in coils.notes[0] + + +def test_device_template_on_demo_data(template_window, data_dir): + window = template_window("device", [data_dir / "a.nc", data_dir / "b.nc"]) + window.set_time(1.0) + section, coils, passive, tf_current, field = window.panels + assert errors_by_panel(window) == {} + assert passive._categories == ["upper plate", "lower plate"] + heights = passive._bars[0].item.opts["height"] + assert heights[0] < -20 and heights[1] < -20 # kA, induced by the ramp-up + # two sources: a line style each, the same color for a kind of shape + tf_lines = [c for c in section.curves if c.label == "TF coil"] + assert len(tf_lines) == 2 + pens = [c.item.opts["pen"] for c in tf_lines] + assert pens[0].color() == pens[1].color() and pens[0].style() != pens[1].style()