From 43f85bd809fddb2bc499e3e1313d5c1c96c97228 Mon Sep 17 00:00:00 2001 From: je-cook <81617086+je-cook@users.noreply.github.com> Date: Tue, 30 Jul 2024 10:13:38 +0100 Subject: [PATCH] =?UTF-8?q?=F0=9F=9A=9A=20Move=20to=20new=20eqdsk=20packag?= =?UTF-8?q?e=20(#3416)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * 🚚 Move to new eqdsk package * 🐛 17 -> 3 * 🎨 From_cocos_index -> from_cocos etc * 🎨 Add full coil specification flag * ✅ Fix tests based on eqdsk changes * 🎨 Small changes * 👽️ Fix based of eqdsk changes * 📝 Qpsi docs * Update eudemo/eudemo_tests/test_tf_coils.py * 🐛 Fix fixed_boundary writer * 🎨 Set default bluemira cocos as constant * ✅ Fix tests * ⬆️ Update eqdsk package version --- bluemira/equilibria/coils/_grouping.py | 3 +- bluemira/equilibria/coils/_tools.py | 6 +- bluemira/equilibria/constants.py | 3 + bluemira/equilibria/equilibrium.py | 207 +++++-- .../equilibria/fem_fixed_boundary/file.py | 16 +- bluemira/equilibria/file.py | 548 ------------------ bluemira/equilibria/find.py | 3 +- bluemira/equilibria/flux_surfaces.py | 15 +- bluemira/equilibria/grid.py | 50 +- bluemira/equilibria/limiter.py | 9 +- bluemira/equilibria/physics.py | 12 +- bluemira/equilibria/plasma.py | 3 +- bluemira/equilibria/profiles.py | 61 +- bluemira/equilibria/run.py | 3 +- bluemira/equilibria/solve.py | 3 +- eudemo/config/build_config.json | 3 +- eudemo/config/fixed_boundary_eqdsk.json | 10 +- eudemo/eudemo/equilibria/_designer.py | 16 +- eudemo/eudemo/pf_coils/designer.py | 4 +- .../eudemo_tests/equilibria/test_designer.py | 5 +- .../ivc/test_divertor_silhouette.py | 2 +- .../eudemo_tests/ivc/test_wall_silhouette.py | 2 +- eudemo/eudemo_tests/test_tf_coils.py | 2 +- pyproject.toml | 1 + requirements.txt | 1 + tests/equilibria/test_constraint_functions.py | 3 +- tests/equilibria/test_constraints.py | 3 +- tests/equilibria/test_equilibrium.py | 33 +- tests/equilibria/test_file.py | 246 -------- tests/equilibria/test_find.py | 13 +- tests/equilibria/test_flux_surfaces.py | 9 +- tests/equilibria/test_harmonics.py | 25 +- tests/equilibria/test_optimiser.py | 5 +- tests/equilibria/test_st_equilibrium.py | 11 +- .../test_advective_transport.py | 5 +- .../test_flux_surface_maker.py | 5 +- .../test_radiation_profile.py | 3 +- 37 files changed, 384 insertions(+), 965 deletions(-) delete mode 100644 bluemira/equilibria/file.py delete mode 100644 tests/equilibria/test_file.py diff --git a/bluemira/equilibria/coils/_grouping.py b/bluemira/equilibria/coils/_grouping.py index a88c13d06b..b6353aeaae 100644 --- a/bluemira/equilibria/coils/_grouping.py +++ b/bluemira/equilibria/coils/_grouping.py @@ -36,10 +36,9 @@ if TYPE_CHECKING: import numpy.typing as npt + from eqdsk import EQDSKInterface from matplotlib.axes import Axes - from bluemira.equilibria.file import EQDSKInterface - def symmetrise_coilset(coilset: CoilSet) -> CoilSet: """ diff --git a/bluemira/equilibria/coils/_tools.py b/bluemira/equilibria/coils/_tools.py index d33c7958ce..65cbe9da10 100644 --- a/bluemira/equilibria/coils/_tools.py +++ b/bluemira/equilibria/coils/_tools.py @@ -12,13 +12,13 @@ from typing import TYPE_CHECKING -if TYPE_CHECKING: - from bluemira.equilibria.coils import CoilSet - import numpy as np from bluemira.magnetostatics.greens import circular_coil_inductance_elliptic, greens_psi +if TYPE_CHECKING: + from bluemira.equilibria.coils import CoilSet + def make_mutual_inductance_matrix(coilset: CoilSet) -> np.ndarray: """ diff --git a/bluemira/equilibria/constants.py b/bluemira/equilibria/constants.py index e841b3f10e..d97ad3e9ab 100644 --- a/bluemira/equilibria/constants.py +++ b/bluemira/equilibria/constants.py @@ -77,3 +77,6 @@ # Minimum discretisation of a coil [m], limit imposed to avoid massive memory # usage COIL_DISCR = 0.05 + +# bluemira default cocos +BLUEMIRA_DEFAULT_COCOS = 3 diff --git a/bluemira/equilibria/equilibrium.py b/bluemira/equilibria/equilibrium.py index 105c46ac1a..28cfb82487 100644 --- a/bluemira/equilibria/equilibrium.py +++ b/bluemira/equilibria/equilibrium.py @@ -19,16 +19,16 @@ import numpy as np import numpy.typing as npt import tabulate +from eqdsk import EQDSKInterface from scipy.optimize import minimize from bluemira.base.constants import MU_0 from bluemira.base.file import get_bluemira_path -from bluemira.base.look_and_feel import bluemira_print_flush +from bluemira.base.look_and_feel import bluemira_print_flush, bluemira_warn from bluemira.equilibria.boundary import FreeBoundary, apply_boundary from bluemira.equilibria.coils import CoilSet, symmetrise_coilset -from bluemira.equilibria.constants import PSI_NORM_TOL +from bluemira.equilibria.constants import BLUEMIRA_DEFAULT_COCOS, PSI_NORM_TOL from bluemira.equilibria.error import EquilibriaError -from bluemira.equilibria.file import EQDSKInterface from bluemira.equilibria.find import ( Opoint, Xpoint, @@ -62,6 +62,7 @@ from bluemira.utilities.tools import abs_rel_difference if TYPE_CHECKING: + from eqdsk.models import Sign from matplotlib.axes import Axes from bluemira.equilibria.find import Lpoint @@ -118,7 +119,13 @@ def set_grid(self, grid: Grid): def _get_eqdsk( cls, filename: Path | str, - ) -> tuple[EQDSKInterface, npt.NDArray[np.float64], CoilSet, Grid, Limiter | None]: + from_cocos: int | None = 11, + to_cocos: int | None = BLUEMIRA_DEFAULT_COCOS, + qpsi_sign: Sign | int | None = None, + *, + full_coil: bool = False, + **kwargs, + ) -> tuple[EQDSKInterface, Grid]: """ Get eqdsk data from file for read in @@ -126,8 +133,18 @@ def _get_eqdsk( ---------- filename: Filename - force_symmetry: - Whether or not to force symmetrisation in the CoilSet + from_cocos: + The COCOS index of the EQDSK file. Used when the determined + COCOS is ambiguous. Will raise if given and not one of + the determined COCOS indices. + to_cocos: + The COCOS index to convert the EQDSK file to. + qpsi_sign: + The sign of the qpsi, required for identification + when qpsi is not present in the file. + full_coil: + Whether the eqdsk dxc and dzc represents + the full coil width or half coil width Returns ------- @@ -142,24 +159,25 @@ def _get_eqdsk( limiter: Limiter instance if any limiters are in file """ - e = EQDSKInterface.from_file(filename) - if "equilibria" in e.name: - psi = e.psi - elif "SCENE" in e.name and not isinstance(cls, Breakdown): - psi = e.psi - e.dxc /= 2 - e.dzc /= 2 - elif "fiesta" in e.name.lower(): - psi = e.psi - else: # CREATE - psi = e.psi / (2 * np.pi) # V.s as opposed to V.s/rad + e = EQDSKInterface.from_file( + filename, + from_cocos=from_cocos, + to_cocos=to_cocos, + qpsi_sign=qpsi_sign, + **kwargs, + ) + if "SCENE" in e.name and not isinstance(cls, Breakdown): + bluemira_warn( + "eqdsk name is SCENE assuming coil dx and dz in the file" + " represent the full width and height of the coil" + ) + full_coil = True + + if full_coil: e.dxc /= 2 e.dzc /= 2 - e.cplasma = abs(e.cplasma) - - grid = Grid.from_eqdsk(e) - return e, psi, grid + return e, Grid.from_eqdsk(e) def to_eqdsk( self, @@ -232,7 +250,16 @@ def __init__( self.filename = filename @classmethod - def from_eqdsk(cls, filename: Path | str): + def from_eqdsk( + cls, + filename: Path | str, + from_cocos: int | None = 11, + to_cocos: int | None = BLUEMIRA_DEFAULT_COCOS, + qpsi_sign: Sign | int | None = None, + *, + full_coil: bool = False, + **kwargs, + ): """ Initialises a Breakdown Object from an eqdsk file. Note that this will involve recalculation of the magnetic flux. @@ -241,20 +268,43 @@ def from_eqdsk(cls, filename: Path | str): ---------- filename: Filename - force_symmetry: - Whether or not to force symmetrisation in the CoilSet - """ - e, psi, grid = super()._get_eqdsk(filename) + from_cocos: + The COCOS index of the EQDSK file. Used when the determined + COCOS is ambiguous. Will raise if given and not one of + the determined COCOS indices. + to_cocos: + The COCOS index to convert the EQDSK file to. + qpsi_sign: + The sign of the qpsi, required for identification + when qpsi is not present in the file. + full_coil: + Whether the eqdsk dxc and dzc represents + the full coil width or half coil width + """ + e, grid = super()._get_eqdsk( + filename, + from_cocos=from_cocos, + to_cocos=to_cocos, + full_coil=full_coil, + qpsi_sign=qpsi_sign, + **kwargs, + ) psi_ax = e.psimag psi_b = e.psibdry lcfs = Coordinates({"x": e.xbdry, "z": e.zbdry}) lcfs.close() - profiles = CustomProfile.from_eqdsk(filename) + profiles = CustomProfile.from_eqdsk(e) cls._eqdsk = e return cls( - grid, lcfs, profiles, psi=psi, psi_ax=psi_ax, psi_b=psi_b, filename=filename + grid, + lcfs, + profiles, + psi=e.psi, + psi_ax=psi_ax, + psi_b=psi_b, + filename=filename, ) def get_LCFS(self) -> Coordinates: @@ -386,10 +436,15 @@ def __init__(self): def _get_eqdsk( cls, filename: Path | str, + from_cocos: int | None = 11, + to_cocos: int | None = BLUEMIRA_DEFAULT_COCOS, + qpsi_sign: Sign | int | None = None, *, user_coils: CoilSet | None = None, force_symmetry: bool = False, - ) -> tuple[EQDSKInterface, npt.NDArray[np.float64], CoilSet, Grid, Limiter | None]: + full_coil: bool = False, + **kwargs, + ) -> tuple[EQDSKInterface, Grid, CoilSet, Limiter | None]: """ Get eqdsk data from file for read in @@ -397,11 +452,23 @@ def _get_eqdsk( ---------- filename: Filename + from_cocos: + The COCOS index of the EQDSK file. Used when the determined + COCOS is ambiguous. Will raise if given and not one of + the determined COCOS indices. + to_cocos: + The COCOS index to convert the EQDSK file to. + qpsi_sign: + The sign of the qpsi, required for identification + when qpsi is not present in the file. user_coils: Coilset provided by the user. Set current, j_max and b_max to zero in user_coils. force_symmetry: Whether or not to force symmetrisation in the CoilSet + full_coil: + Whether the eqdsk dxc and dzc represents + the full coil width or half coil width Returns ------- @@ -416,19 +483,21 @@ def _get_eqdsk( limiter: Limiter instance if any limiters are in file """ - e, psi, grid = super()._get_eqdsk(filename) + e, grid = super()._get_eqdsk( + filename, + from_cocos=from_cocos, + to_cocos=to_cocos, + qpsi_sign=qpsi_sign, + full_coil=full_coil, + **kwargs, + ) coilset = user_coils if user_coils is not None else CoilSet.from_group_vecs(e) if force_symmetry: coilset = symmetrise_coilset(coilset) - if e.nlim == 0: - limiter = None - elif e.nlim < 5: # noqa: PLR2004 - limiter = Limiter(e.xlim, e.zlim) - else: - limiter = None # CREATE.. + limiter = None if e.nlim > 5 or e.nlim == 0 else Limiter(e.xlim, e.zlim) # noqa: PLR2004 - return e, psi, coilset, grid, limiter + return e, grid, coilset, limiter def _remap_greens(self): """ @@ -557,9 +626,14 @@ def __init__( def from_eqdsk( cls, filename: Path | str, + from_cocos: int | None = 11, + to_cocos: int | None = BLUEMIRA_DEFAULT_COCOS, + qpsi_sign: Sign | int | None = None, *, force_symmetry: bool, user_coils: CoilSet | None = None, + full_coil: bool = False, + **kwargs, ): """ Initialises a Breakdown Object from an eqdsk file. Note that this @@ -569,16 +643,35 @@ def from_eqdsk( ---------- filename: Filename + from_cocos: + The COCOS index of the EQDSK file. Used when the determined + COCOS is ambiguous. Will raise if given and not one of + the determined COCOS indices. + to_cocos: + The COCOS index to convert the EQDSK file to. + qpsi_sign: + The sign of the qpsi, required for identification + when qpsi is not present in the file. force_symmetry: Whether or not to force symmetrisation in the CoilSet user_coils: Coilset provided by the user. Set current, j_max and b_max to zero in user_coils. + full_coil: + Whether the eqdsk dxc and dzc represents + the full coil width or half coil width """ - cls._eqdsk, psi, coilset, grid, limiter = super()._get_eqdsk( - filename, force_symmetry=force_symmetry, user_coils=user_coils + cls._eqdsk, grid, coilset, limiter = super()._get_eqdsk( + filename, + from_cocos, + to_cocos, + qpsi_sign=qpsi_sign, + force_symmetry=force_symmetry, + user_coils=user_coils, + full_coil=full_coil, + **kwargs, ) - return cls(coilset, grid, limiter=limiter, psi=psi, filename=filename) + return cls(coilset, grid, limiter=limiter, psi=cls._eqdsk.psi, filename=filename) def to_dict(self) -> dict[str, Any]: """ @@ -880,9 +973,14 @@ def __init__( def from_eqdsk( cls, filename: Path | str, + from_cocos: int | None = 11, + to_cocos: int | None = BLUEMIRA_DEFAULT_COCOS, + qpsi_sign: Sign | int | None = None, *, force_symmetry: bool = False, user_coils: CoilSet | None = None, + full_coil: bool = False, + **kwargs, ): """ Initialises an Equilibrium Object from an eqdsk file. Note that this @@ -896,24 +994,41 @@ def from_eqdsk( ---------- filename: Filename + from_cocos: + The COCOS index of the EQDSK file. Used when the determined + COCOS is ambiguous. Will raise if given and not one of + the determined COCOS indices. + to_cocos: + The COCOS index to convert the EQDSK file to. + qpsi_sign: + The sign of the qpsi, required for identification + when qpsi is not present in the file. force_symmetry: Whether or not to force symmetrisation in the CoilSet user_coils: Coilset provided by the user. Set current, j_max and b_max to zero in user_coils. + full_coil: + Whether the eqdsk dxc and dzc represents + the full coil width or half coil width """ - e, psi, coilset, grid, limiter = super()._get_eqdsk( + e, grid, coilset, limiter = super()._get_eqdsk( filename, + from_cocos=from_cocos, + to_cocos=to_cocos, + qpsi_sign=qpsi_sign, force_symmetry=force_symmetry, user_coils=user_coils, + full_coil=full_coil, + **kwargs, ) - profiles = CustomProfile.from_eqdsk(filename) + profiles = CustomProfile.from_eqdsk(e) cls._eqdsk = e - o_points, x_points = find_OX_points(grid.x, grid.z, psi, limiter=limiter) - jtor = profiles.jtor(grid.x, grid.z, psi, o_points=o_points, x_points=x_points) + o_points, x_points = find_OX_points(grid.x, grid.z, e.psi, limiter=limiter) + jtor = profiles.jtor(grid.x, grid.z, e.psi, o_points=o_points, x_points=x_points) return cls( coilset, @@ -921,7 +1036,7 @@ def from_eqdsk( profiles=profiles, vcontrol=None, limiter=limiter, - psi=psi, + psi=e.psi, jtor=jtor, filename=filename, ) @@ -1010,7 +1125,7 @@ def to_dict(self, qpsi_calcmode: int = 0) -> dict[str, Any]: def to_eqdsk( self, filename: Path | str, - header: str = "BP_equilibria", + header: str = "bluemira_equilibria", directory: str | None = None, filetype: str = "json", qpsi_calcmode: int = 0, diff --git a/bluemira/equilibria/fem_fixed_boundary/file.py b/bluemira/equilibria/fem_fixed_boundary/file.py index 21a286c928..c759179e47 100644 --- a/bluemira/equilibria/fem_fixed_boundary/file.py +++ b/bluemira/equilibria/fem_fixed_boundary/file.py @@ -9,21 +9,27 @@ File saving for fixed boundary equilibrium """ +from __future__ import annotations + +from typing import TYPE_CHECKING + import numpy as np +from eqdsk import EQDSKInterface from scipy.integrate import quad, quadrature from bluemira.base.look_and_feel import bluemira_warn -from bluemira.equilibria.fem_fixed_boundary.fem_magnetostatic_2D import ( - FixedBoundaryEquilibrium, -) from bluemira.equilibria.fem_fixed_boundary.utilities import ( _interpolate_profile, find_magnetic_axis, get_mesh_boundary, ) -from bluemira.equilibria.file import EQDSKInterface from bluemira.equilibria.grid import Grid +if TYPE_CHECKING: + from bluemira.equilibria.fem_fixed_boundary.fem_magnetostatic_2D import ( + FixedBoundaryEquilibrium, + ) + def _pressure_profile(pprime, psi_norm, psi_mag): pressure = np.zeros(len(psi_norm)) @@ -139,7 +145,7 @@ def save_fixed_boundary_to_file( pressure=pressure, pprime=p_prime, psi=psi, - psibdry=np.zeros(nbdry), + psibdry=0, psimag=psi_mag, xbdry=xbdry, xc=np.array([]), diff --git a/bluemira/equilibria/file.py b/bluemira/equilibria/file.py deleted file mode 100644 index 9b9ba4a42c..0000000000 --- a/bluemira/equilibria/file.py +++ /dev/null @@ -1,548 +0,0 @@ -# SPDX-FileCopyrightText: 2021-present M. Coleman, J. Cook, F. Franza -# SPDX-FileCopyrightText: 2021-present I.A. Maione, S. McIntosh -# SPDX-FileCopyrightText: 2021-present J. Morris, D. Short -# -# SPDX-License-Identifier: LGPL-2.1-or-later -""" -Input and output file interface. EQDSK and json. NOTE: jsons are better :) -""" - -from __future__ import annotations - -import json -import time -from dataclasses import asdict, dataclass -from pathlib import Path -from typing import Any - -import fortranformat as ff -import numpy as np -import numpy.typing as npt - -from bluemira.base.look_and_feel import bluemira_warn -from bluemira.utilities.tools import is_num, json_writer - -EQDSK_EXTENSIONS = [".eqdsk", ".eqdsk_out", ".geqdsk"] - - -@dataclass(repr=False) -class EQDSKInterface: - """ - Container for data from an EQDSK file. - - Inspired by an EQDSK reader originally developed by B. Dudson: - https://github.com/bendudson/pyTokamak/blob/master/tokamak/formats/geqdsk.py - - The G-EQDSK file format is described here: - https://fusion.gat.com/conferences/snowmass/working/mfe/physics/p3/equilibria/g_eqdsk_s.pdf - - Notes - ----- - G-EQDSK is from the 1980's and EQDSK files should generally only be - read and not written. New equilibria should really just be saved as - JSON files. - - Poloidal magnetic flux units not enforced here! - - Plasma current direction is not enforced here! - """ - - bcentre: float - """Magnetic field at the reference radius [T].""" - cplasma: float - """Plasma current [A].""" - dxc: npt.NDArray[np.float64] - """X half-thicknesses of the coils [m].""" - dzc: npt.NDArray[np.float64] - """Z half-thicknesses of the coils [m].""" - ffprime: npt.NDArray[np.float64] - """FF' function on 1-D flux grid [m.T^2/V.s/rad].""" - fpol: npt.NDArray[np.float64] - """Poloidal current function f = R*B on 1-D flux [T.m].""" - Ic: npt.NDArray[np.float64] - """Coil currents [A].""" - name: str - """Name of the equilibrium EQDSK [dimensionless].""" - nbdry: int - """Number of boundary points [dimensionless].""" - ncoil: int - """Number of coils [dimensionless].""" - nlim: int - """Number of limiters [dimensionless].""" - nx: int - """Number of grid points in the radial direction [dimensionless].""" - nz: int - """Number of grid points in the vertical direction [dimensionless].""" - pprime: npt.NDArray[np.float64] - """P' function on 1-D flux grid [N/m^2/V.s/rad].""" - pressure: npt.NDArray[np.float64] - """Plasma pressure function on 1-D flux grid [N/m^2].""" - psi: npt.NDArray[np.float64] - """Poloidal magnetic flux on the 2-D grid [V.s/rad].""" - psibdry: float - """Poloidal flux at the magnetic axis [V.s/rad].""" - psimag: float - """Poloidal flux at the magnetic axis [V.s/rad].""" - xbdry: npt.NDArray[np.float64] - """X coordinates of the plasma boundary [m].""" - xc: npt.NDArray[np.float64] - """X coordinates of the coils [m].""" - xcentre: float - """Radius of the reference toroidal magnetic [m].""" - xdim: float - """Horizontal dimension of the spatial grid [m].""" - xgrid1: float - """Minimum radius of the spatial grid [m].""" - xlim: npt.NDArray[np.float64] - """X coordinates of the limiters [m].""" - xmag: float - """Radius of the magnetic axis [m].""" - zbdry: npt.NDArray[np.float64] - """Z coordinates of the plasma boundary [m].""" - zc: npt.NDArray[np.float64] - """Z coordinates of the coils [m].""" - zdim: float - """Vertical dimension of the spatial grid [m].""" - zlim: npt.NDArray[np.float64] - """Z coordinates of the limiters [m].""" - zmag: float - """Z coordinate of the magnetic axis [m].""" - zmid: float - """Z coordinate of the middle of the spatial grid [m].""" - x: npt.NDArray[np.float64] | None = None - """X 1-D vector [m] (calculated if not given).""" - z: npt.NDArray[np.float64] | None = None - """Z 1-D vector [m] (calculated if not given).""" - psinorm: npt.NDArray[np.float64] | None = None - """Normalised psi vector [A] (calculated if not given).""" - qpsi: npt.NDArray[np.float64] | None = None - """Safety factor values on the 1-D flux grid [dimensionless].""" - file_name: str | None = None - """The EQDSK file the data originates from.""" - coil_names: list[str] | None = None - """Name of the coils""" - coil_types: list[str] | None = None - """Type of the coils""" - - def __post_init__(self): - """Calculate derived parameters if they're not given.""" - if self.x is None: - self.x = _derive_x(self.xgrid1, self.xdim, self.nx) - if self.z is None: - self.z = _derive_z(self.zmid, self.zdim, self.nz) - if self.psinorm is None: - self.psinorm = _derive_psinorm(self.fpol) - - @classmethod - def from_file(cls, file_path: Path | str) -> EQDSKInterface: - """ - Create an EQDSKInterface object from a file. - - Parameters - ---------- - file_path: - Path to a file of one of the following formats: - - * JSON - * eqdsk - * eqdsk_out - * geqdsk - - Returns - ------- - An instance of this class containing the EQDSK file's data. - """ - file_extension = Path(file_path).suffix - file_name = Path(file_path).name - if file_extension.lower() in EQDSK_EXTENSIONS: - return cls(file_name=file_name, **_read_eqdsk(file_path)) - if file_extension.lower() == ".json": - return cls(file_name=file_name, **_read_json(file_path)) - raise ValueError(f"Unrecognised file format '{file_extension}'.") - - def to_dict(self) -> dict: - """Return a dictionary of the EQDSK data.""" - d = asdict(self) - # Remove the file name as this is metadata, not EQDSK data - del d["file_name"] - return d - - def write( - self, - file_path: str, - file_format: str = "json", - json_kwargs: dict | None = None, - ): - """ - Write the EQDSK data to file in the given format. - - Parameters - ---------- - file_path: - Path to where the file should be written. - file_format: - The format to save the file in. One of 'json', 'eqdsk', or - 'geqdsk'. - json_kwargs: - Key word arguments to pass to the ``json.dump`` call. Only - used if ``format`` is 'json'. - """ - if file_format == "json": - json_kwargs = {} if json_kwargs is None else json_kwargs - json_writer(self.to_dict(), file_path, **json_kwargs) - elif file_format in {"eqdsk", "geqdsk"}: - bluemira_warn( - "You are in the 21st century. Are you sure you want to be making an" - " EDQSK in this day and age?" - ) - _write_eqdsk(Path(file_path).as_posix(), self.to_dict()) - - def update(self, eqdsk_data: dict[str, Any]): - """ - Update this object's data with values from a dictionary. - - Parameters - ---------- - eqdsk_data: - A dict containing the new eqdsk data. - - Raises - ------ - ValueError - If a key in ``eqdsk_data`` does not correspond to an - attribute of this class. - """ - for key, value in eqdsk_data.items(): - if hasattr(self, key): - setattr(self, key, value) - else: - raise ValueError( - f"Cannot update EQDSKInterface from dict. Unrecognised key '{key}'." - ) - - -def _read_json(file) -> dict[str, Any]: - if isinstance(file, Path | str): - with open(file) as f_h: - return _read_json(f_h) - - data = json.load(file) - data_has_pnorm = False - data_has_psinorm = False - for k, value in data.items(): - if isinstance(value, list) and k not in {"coil_type", "coil_names"}: - data[k] = np.asarray(value) - data_has_pnorm |= k == "pnorm" - data_has_psinorm |= k == "psinorm" - - # For backward compatibility where 'psinorm' was sometimes 'pnorm' - if data_has_pnorm: - if data_has_psinorm: - del data["pnorm"] - else: - data["psinorm"] = data.pop("pnorm") - - return data - - -def _read_array(tokens, n, name="Unknown"): - data = np.zeros([n]) - try: - for i in np.arange(n): - data[i] = float(next(tokens)) - except StopIteration: - raise ValueError(f"Failed reading array {name} of size {n}") from None - return data - - -def _read_2d_array(tokens, n_x, n_y, name="Unknown"): - data = np.zeros([n_y, n_x]) - for i in np.arange(n_y): - data[i, :] = _read_array(tokens, n_x, name + "[" + str(i) + "]") - return np.transpose(data) - - -def _eqdsk_generator(file): - """ - Transforms a file object into a generator, following G-EQDSK number - conventions - - Parameters - ---------- - file: - The file to read - - Returns - ------- - The generator of the file handle being read - """ - while True: - line = file.readline() - if not line: - break - - # Distinguish negative/positive numbers from negative/positive exponent - if "E" in line or "e" in line: - line = line.replace("E-", "*") - line = line.replace("e-", "*") - line = line.replace("-", " -") - line = line.replace("*", "e-") - line = line.replace("E+", "*") - line = line.replace("e+", "*") - line = line.replace("+", " ") - line = line.replace("*", "e+") - generator_list = line.split() - yield from generator_list - - -def _read_eqdsk(file) -> dict: - if isinstance(file, Path | str): - with open(file) as f_handle: - return _read_eqdsk(f_handle) - - description = file.readline() - if not description: - raise OSError(f"Could not read the file '{file}'.") - description = description.split() - - ints = [value for value in description if is_num(value)] - - if len(ints) < 3: # noqa: PLR2004 - raise OSError( - f"Should be at least 3 numbers in the first line of the EQDSK {file}." - ) - - data = {} - n_x = int(ints[-2]) - n_z = int(ints[-1]) - data["name"] = description[0] - data["nx"] = n_x - data["nz"] = n_z - - tokens = _eqdsk_generator(file) - for name in [ - "xdim", - "zdim", - "xcentre", - "xgrid1", - "zmid", - "xmag", - "zmag", - "psimag", - "psibdry", - "bcentre", - "cplasma", - "psimag", - None, - "xmag", - None, - "zmag", - None, - "psibdry", - None, - None, - ]: - if name is not None: # Lots of dummies and duplication - data[name] = float(next(tokens)) - else: - next(tokens) # Dummy - - for name in ["fpol", "pressure", "ffprime", "pprime"]: - data[name] = _read_array(tokens, n_x, name) - - data["psi"] = _read_2d_array(tokens, n_x, n_z, "psi") - data["qpsi"] = _read_array(tokens, n_x, "qpsi") - nbdry = int(next(tokens)) - nlim = int(next(tokens)) - data["nbdry"] = nbdry - data["nlim"] = nlim - - xbdry = np.zeros(nbdry) - zbdry = np.zeros(nbdry) - for i in range(nbdry): - xbdry[i] = float(next(tokens)) - zbdry[i] = float(next(tokens)) - data["xbdry"] = xbdry - data["zbdry"] = zbdry - - xlim = np.zeros(nlim) - zlim = np.zeros(nlim) - for i in range(nlim): - xlim[i] = float(next(tokens)) - zlim[i] = float(next(tokens)) - data["xlim"] = xlim - data["zlim"] = zlim - - try: - ncoil = int(next(tokens)) - except StopIteration: # No coils in file - ncoil = 0 - - x_c = np.zeros(ncoil) - z_c = np.zeros(ncoil) - dxc = np.zeros(ncoil) - dzc = np.zeros(ncoil) - i_c = np.zeros(ncoil) - for i in range(ncoil): - x_c[i] = float(next(tokens)) - z_c[i] = float(next(tokens)) - dxc[i] = float(next(tokens)) - dzc[i] = float(next(tokens)) - i_c[i] = float(next(tokens)) - data["ncoil"] = ncoil - data["xc"] = x_c - data["zc"] = z_c - data["dxc"] = dxc - data["dzc"] = dzc - data["Ic"] = i_c - - # Additional utility data - data["x"] = _derive_x(data["xgrid1"], data["xdim"], data["nx"]) - data["z"] = _derive_z(data["zmid"], data["zdim"], data["nz"]) - data["psinorm"] = _derive_psinorm(data["fpol"]) - return data - - -def _derive_x(xgrid1, xdim, nx): - return np.linspace(xgrid1, xgrid1 + xdim, nx) - - -def _derive_z(zmid, zdim, nz): - return np.linspace(zmid - zdim / 2, zmid + zdim / 2, nz) - - -def _derive_psinorm(fpol): - return np.linspace(0, 1, len(fpol)) - - -def _write_eqdsk(file: Path | str, data: dict): - """ - Writes data out to a text file in G-EQDSK format. - - Parameters - ---------- - file: - The full path string of the file to be created - data: - Dictionary of EQDSK data. - """ - if isinstance(file, str | Path): - file = Path(file) - if file.suffix not in EQDSK_EXTENSIONS: - file = Path(file).with_suffix("").with_suffix(".eqdsk") - with open(file, "w") as f_handle: - return _write_eqdsk(f_handle, data) - - def write_header( - fortran_format: ff.FortranRecordWriter, id_string: str, var_list: list[str] - ): - """ - Writes G-EQDSK header out to file. - - Parameters - ---------- - fortran_format: - FortranRecordWriter object for Fortran format edit descriptor - to be used for header output. - id_string: - String containing name of file to be used as identification - string. Will be trimmed if length exceeds 39 characters, - so it will fit within the permitted header length of the - GEQDSK specification when a timestamp is added. - var_list: - List of names of keys in EQDSKInterface.data identifying - variables to add to the header following the id_string. - Empty strings will be recorded as 0. - """ - line = [id_string] - line += [data[v] if v else 0 for v in var_list] - file.write(fortran_format.write(line)) - file.write("\n") - - def write_line(fortran_format: ff.FortranRecordWriter, var_list: list[str]): - """ - Writes a line of variable values out to a G-EQDSK file. - - Parameters - ---------- - fortran_format: - FortranRecordWriter object for Fortran format edit descriptor - to be used for the format of the line output. - var_list: - List of names of keys in EQDSKInterface.data identifying - variables to added to the current line. - Empty strings will be recorded as 0. - """ - line = [data[v] if v else 0 for v in var_list] - file.write(fortran_format.write(line)) - file.write("\n") - - def write_array( - fortran_format: ff.FortranRecordWriter, array: npt.NDArray[np.float64] - ): - """ - Writes a numpy array out to a G-EQDSK file. - - Parameters - ---------- - fortran_format: - FortranRecordWriter object for Fortran format edit descriptor - to be used for the format of the line output. - array: - Numpy array of variables to be written to file. - Array will be flattened in column-major (Fortran) - order if is more than one-dimensional. - """ - if array.ndim > 1: - flat_array = array.flatten(order="F") - file.write(fortran_format.write(flat_array)) - else: - file.write(fortran_format.write(array)) - file.write("\n") - - # Create id string for file comprising of timestamp and trimmed filename - # that fits the 48 character limit of strings in EQDSK headers. - timestamp = time.strftime("%d%m%Y") - trimmed_name = data["name"][0 : 48 - len(timestamp) - 1] - file_id_string = f"{trimmed_name}_{timestamp}" - - # Define dummy data for qpsi if it has not been previously defined. - qpsi = np.zeros(data["nx"]) if data["qpsi"] is None else data["qpsi"] - - # Create array containing coilset information. - coil = np.zeros(5 * data["ncoil"]) - for i, value in enumerate(["xc", "zc", "dxc", "dzc", "Ic"]): - coil[i::5] = data[value] - - # Create FortranRecordWriter objects with the Fortran format - # edit descriptors to be used in the G-EQDSK output. - f2000 = ff.FortranRecordWriter("a48,3i4") - f2020 = ff.FortranRecordWriter("5e16.9") - f2022 = ff.FortranRecordWriter("2i5") - fCSTM = ff.FortranRecordWriter("i5") - - # Write header in f2000 (6a8,3i4) format. - write_header(f2000, file_id_string, ["", "nx", "nz"]) - # Write out lines containing floats in f2020 (5e16.9) format. - write_line(f2020, ["xdim", "zdim", "xcentre", "xgrid1", "zmid"]) - write_line(f2020, ["xmag", "zmag", "psimag", "psibdry", "bcentre"]) - write_line(f2020, ["cplasma", "psimag", "", "xmag", ""]) - write_line(f2020, ["zmag", "", "psibdry", "", ""]) - # Write out arrays in in f2020 (5e16.9) format. - write_array(f2020, data["fpol"]) - write_array(f2020, data["pressure"]) - write_array(f2020, data["ffprime"]) - write_array(f2020, data["pprime"]) - write_array(f2020, data["psi"]) - write_array(f2020, qpsi) - # Write out number of boundary points and limiters f2022 (2i5) format. - write_line(f2022, ["nbdry", "nlim"]) - # Write out boundary point and limiter data as array of ordered pairs. - write_array(f2020, np.array([data["xbdry"], data["zbdry"]])) - write_array(f2020, np.array([data["xlim"], data["zlim"]])) - - # Output of coilset information. This is an extension to the - # regular eqdsk format. - write_line(fCSTM, ["ncoil"]) - write_array(f2020, coil) - return None diff --git a/bluemira/equilibria/find.py b/bluemira/equilibria/find.py index 10842a2106..730ed1a160 100644 --- a/bluemira/equilibria/find.py +++ b/bluemira/equilibria/find.py @@ -14,7 +14,6 @@ import numba as nb import numpy as np -import numpy.typing as npt from contourpy import LineType, contour_generator from scipy.interpolate import RectBivariateSpline @@ -31,6 +30,8 @@ if TYPE_CHECKING: from collections.abc import Iterable + import numpy.typing as npt + from bluemira.equilibria.limiter import Limiter __all__ = [ diff --git a/bluemira/equilibria/flux_surfaces.py b/bluemira/equilibria/flux_surfaces.py index 60d6e122d1..14a5a90bee 100644 --- a/bluemira/equilibria/flux_surfaces.py +++ b/bluemira/equilibria/flux_surfaces.py @@ -15,14 +15,6 @@ from functools import lru_cache from typing import TYPE_CHECKING -from bluemira.utilities.tools import floatify - -if TYPE_CHECKING: - from collections.abc import Iterable - - from bluemira.equilibria.equilibrium import Equilibrium - from bluemira.equilibria.find import PsiPoint - import matplotlib.pyplot as plt import numba as nb import numpy as np @@ -44,6 +36,13 @@ ) from bluemira.geometry.plane import BluemiraPlane from bluemira.geometry.tools import _signed_distance_2D +from bluemira.utilities.tools import floatify + +if TYPE_CHECKING: + from collections.abc import Iterable + + from bluemira.equilibria.equilibrium import Equilibrium + from bluemira.equilibria.find import PsiPoint @nb.jit(nopython=True, cache=True) diff --git a/bluemira/equilibria/grid.py b/bluemira/equilibria/grid.py index 9912dbdc17..8abb9b226a 100644 --- a/bluemira/equilibria/grid.py +++ b/bluemira/equilibria/grid.py @@ -21,7 +21,11 @@ from bluemira.geometry.coordinates import get_area_2d, get_centroid_2d if TYPE_CHECKING: - from bluemira.equilibria.file import EQDSKInterface + from collections.abc import Iterable + + import numpy.typing as npt + from eqdsk import EQDSKInterface + __all__ = ["Grid", "integrate_dx_dz", "revolved_volume", "volume_integral"] @@ -33,17 +37,17 @@ class Grid: Parameters ---------- - x_min: float > 0 - Minimum x grid coordinate [m] - x_max: float + x_min: + Minimum x grid coordinate (>0) [m] + x_max: Maximum x grid coordinate [m] - z_min: float + z_min: Minimum z grid coordinate [m] - z_max: float + z_max: Maximum z grid coordinate [m] - nx: int + nx: Number of x grid points - nz: int + nz: Number of z grid points """ @@ -68,12 +72,13 @@ class Grid: "z_size", ) - def __init__(self, x_min, x_max, z_min, z_max, nx, nz): + def __init__( + self, x_min: float, x_max: float, z_min: float, z_max: float, nx: int, nz: int + ): if x_min == x_max or z_min == z_max: raise EquilibriaError("Invalid Grid dimensions specified.") if x_min > x_max: - print() # stdout flusher # noqa: T201 bluemira_warn( f"x_min should be < x_max {x_min:.2f} > {x_max:.2f}. Switching x_min and" " x_max." @@ -81,7 +86,6 @@ def __init__(self, x_min, x_max, z_min, z_max, nx, nz): x_min, x_max = x_max, x_min if z_min > z_max: - print() # stdout flusher # noqa: T201 bluemira_warn( f"z_min should be < z_max {z_min:.2f} > {z_max:.2f}. Switching z_min and" " z_max." @@ -92,14 +96,12 @@ def __init__(self, x_min, x_max, z_min, z_max, nx, nz): x_min = X_AXIS_MIN if nx < MIN_N_DISCR: - print() # stdout flusher # noqa: T201 bluemira_warn( f"Insufficient nx discretisation: {nx}, setting to {MIN_N_DISCR}." ) nx = MIN_N_DISCR if nz < MIN_N_DISCR: - print() # stdout flusher # noqa: T201 bluemira_warn( f"Insufficient nx discretisation: {nz}, setting to {MIN_N_DISCR}." ) @@ -131,7 +133,7 @@ def __init__(self, x_min, x_max, z_min, z_max, nx, nz): ]) @classmethod - def from_eqdict(cls, e): + def from_eqdict(cls, e) -> Grid: """ Initialise a Grid object from an EQDSK dictionary. @@ -150,13 +152,13 @@ def from_eqdict(cls, e): ) @classmethod - def from_eqdsk(cls, e: EQDSKInterface): + def from_eqdsk(cls, e: EQDSKInterface) -> Grid: """ Initialise a Grid object from an EQDSKInterface. Parameters ---------- - e: EQDSKInterface + e: """ return cls( @@ -168,15 +170,15 @@ def from_eqdsk(cls, e: EQDSKInterface): e.nz, ) - def point_inside(self, x, z=None): + def point_inside(self, x: float | Iterable[float], z: float | None = None) -> bool: """ Determine if a point is inside the rectangular grid (includes edges). Parameters ---------- - x: Union[float, Iterable] + x: The x coordinate of the point. Or the 2-D point. - z: Optional[float] + z: The z coordinate of the point Returns @@ -193,20 +195,22 @@ def point_inside(self, x, z=None): and (z <= self.z_max) ) - def distance_to(self, x, z=None): + def distance_to( + self, x: float | Iterable[float], z: float | None = None + ) -> npt.NDArray[np.float64]: """ Get the distances of a point to the edges of the Grid. Parameters ---------- - x: Union[float, Iterable] + x: The x coordinate of the point. Or the 2-D point. - z: Optional[float] + z: The z coordinate of the point Returns ------- - distances: np.ndarray + distances: Distances to the edges of the Grid. """ if z is None: diff --git a/bluemira/equilibria/limiter.py b/bluemira/equilibria/limiter.py index 104ac0544c..9b76feac9e 100644 --- a/bluemira/equilibria/limiter.py +++ b/bluemira/equilibria/limiter.py @@ -8,14 +8,19 @@ Limiter object class """ +from __future__ import annotations + from itertools import cycle +from typing import TYPE_CHECKING import numpy as np -import numpy.typing as npt -from matplotlib.pyplot import Axes from bluemira.equilibria.plotting import LimiterPlotter +if TYPE_CHECKING: + import numpy.typing as npt + from matplotlib.pyplot import Axes + __all__ = ["Limiter"] diff --git a/bluemira/equilibria/physics.py b/bluemira/equilibria/physics.py index 75fd7014d9..99aaca32c7 100644 --- a/bluemira/equilibria/physics.py +++ b/bluemira/equilibria/physics.py @@ -12,12 +12,6 @@ from typing import TYPE_CHECKING -if TYPE_CHECKING: - from collections.abc import Iterable - - from bluemira.equilibria.equilibrium import Equilibrium - from bluemira.equilibria.find import Opoint, Xpoint - import numpy as np import numpy.typing as npt from scipy.interpolate import RectBivariateSpline @@ -26,6 +20,12 @@ from bluemira.equilibria.find import in_plasma from bluemira.equilibria.grid import revolved_volume, volume_integral +if TYPE_CHECKING: + from collections.abc import Iterable + + from bluemira.equilibria.equilibrium import Equilibrium + from bluemira.equilibria.find import Opoint, Xpoint + def calc_psi_norm( psi: npt.ArrayLike, opsi: float, xpsi: float diff --git a/bluemira/equilibria/plasma.py b/bluemira/equilibria/plasma.py index 478216a63c..b545bc2dcf 100644 --- a/bluemira/equilibria/plasma.py +++ b/bluemira/equilibria/plasma.py @@ -13,7 +13,6 @@ from typing import TYPE_CHECKING import numpy as np -import numpy.typing as npt from scipy.interpolate import RectBivariateSpline from bluemira.equilibria.constants import J_TOR_MIN @@ -23,6 +22,8 @@ from bluemira.utilities.tools import floatify if TYPE_CHECKING: + import numpy.typing as npt + from bluemira.equilibria.grid import Grid diff --git a/bluemira/equilibria/profiles.py b/bluemira/equilibria/profiles.py index 2264418d65..a25b627cfb 100644 --- a/bluemira/equilibria/profiles.py +++ b/bluemira/equilibria/profiles.py @@ -16,13 +16,14 @@ import numba as nb import numpy as np import numpy.typing as npt +from eqdsk import EQDSKInterface from scipy.integrate import quad from scipy.interpolate import interp1d from scipy.optimize import curve_fit from bluemira.base.constants import MU_0 +from bluemira.equilibria.constants import BLUEMIRA_DEFAULT_COCOS from bluemira.equilibria.error import EquilibriaError -from bluemira.equilibria.file import EQDSKInterface from bluemira.equilibria.find import find_LCFS_separatrix, in_plasma, in_zone from bluemira.equilibria.grid import integrate_dx_dz, revolved_volume, volume_integral from bluemira.equilibria.plotting import ProfilePlotter @@ -31,7 +32,9 @@ from collections.abc import Callable from pathlib import Path - from bluemira.equilibria.find import Opoint, Optional, Xpoint + from eqdsk.models import Sign + + from bluemira.equilibria.find import Opoint, Xpoint __all__ = [ "BetaIpProfile", @@ -50,7 +53,7 @@ def fitfunc( func: Callable[[float], float], data: npt.NDArray[np.float64], - order: Optional[int] = None, + order: int | None = None, ) -> npt.NDArray[np.float64]: """ Uses scipy's curve_fit to fit 1-D data to a custom function @@ -217,7 +220,7 @@ class ShapeFunction: """ @classmethod - def from_datafit(cls, data: npt.NDArray[np.float64], order: Optional[int] = None): + def from_datafit(cls, data: npt.NDArray[np.float64], order: int | None = None): """ Defines function from a dataset, fit using scipy curve_fit """ @@ -417,7 +420,7 @@ def _jtor( psi: npt.NDArray[np.float64], o_points: list[Opoint], x_points: list[Xpoint], - lcfs: Optional[np.ndarray] = None, + lcfs: np.ndarray | None = None, ) -> tuple[float, float, npt.NDArray[np.float64]]: """ Do-not-repeat-yourself utility @@ -534,7 +537,7 @@ def __init__( I_p: float, R_0: float, B_0: float, - shape: Optional[ShapeFunction] = None, + shape: ShapeFunction | None = None, ): self.betap = betap self.I_p = I_p @@ -672,7 +675,7 @@ def __init__( I_p: float, R_0: float, B_0: float, - shape: Optional[ShapeFunction] = None, + shape: ShapeFunction | None = None, li_rel_tol: float = 0.015, li_min_iter: int = 5, ): @@ -708,9 +711,9 @@ def __init__( ffprime_func: npt.NDArray[np.float64] | Callable[[float]] | float, R_0: float, B_0: float, - p_func: Optional[npt.NDArray[np.float64] | Callable[[float]] | float] = None, - f_func: Optional[npt.NDArray[np.float64] | Callable[[float]] | float] = None, - I_p: Optional[float] = None, + p_func: npt.NDArray[np.float64] | Callable[[float]] | float | None = None, + f_func: npt.NDArray[np.float64] | Callable[[float]] | float | None = None, + I_p: float | None = None, ): self._pprime_in = self.parse_to_callable(pprime_func) self._ffprime_in = self.parse_to_callable(ffprime_func) @@ -799,17 +802,37 @@ def fRBpol(self, psinorm: npt.ArrayLike) -> float | npt.NDArray[np.float64]: return super().fRBpol(psinorm) @classmethod - def from_eqdsk(cls, filename: Path | str) -> CustomProfile: + def from_eqdsk_file( + cls, + filename: Path | str, + from_cocos: int | None = 11, + to_cocos: int | None = BLUEMIRA_DEFAULT_COCOS, + qpsi_sign: Sign | int | None = None, + **kwargs, + ) -> CustomProfile: """ Initialises a CustomProfile object from an eqdsk file """ - e = EQDSKInterface.from_file(filename) + e = EQDSKInterface.from_file( + filename, + from_cocos=from_cocos, + to_cocos=to_cocos, + qpsi_sign=qpsi_sign, + **kwargs, + ) + return cls.from_eqdsk(e) + + @classmethod + def from_eqdsk(cls, eq: EQDSKInterface) -> CustomProfile: + """ + Initialises a CustomProfile object from an eqdsk object + """ return cls( - e.pprime, - e.ffprime, - R_0=e.xcentre, - B_0=abs(e.bcentre), - p_func=e.pressure, - f_func=e.fpol, - I_p=abs(e.cplasma), + eq.pprime, + eq.ffprime, + R_0=eq.xcentre, + B_0=eq.bcentre, + p_func=eq.pressure, + f_func=eq.fpol, + I_p=eq.cplasma, ) diff --git a/bluemira/equilibria/run.py b/bluemira/equilibria/run.py index 9125f2f952..9d45d252e1 100644 --- a/bluemira/equilibria/run.py +++ b/bluemira/equilibria/run.py @@ -17,7 +17,6 @@ import matplotlib.pyplot as plt import numpy as np -import numpy.typing as npt from bluemira.base.look_and_feel import bluemira_print, bluemira_warn from bluemira.base.parameter_frame import Parameter, ParameterFrame @@ -51,6 +50,8 @@ if TYPE_CHECKING: from collections.abc import Iterable + import numpy.typing as npt + from bluemira.equilibria.coils import CoilSet from bluemira.equilibria.grid import Grid from bluemira.equilibria.limiter import Limiter diff --git a/bluemira/equilibria/solve.py b/bluemira/equilibria/solve.py index 7e90999c52..08aca99bcd 100644 --- a/bluemira/equilibria/solve.py +++ b/bluemira/equilibria/solve.py @@ -15,7 +15,6 @@ import matplotlib.pyplot as plt import numpy as np -import numpy.typing as npt from bluemira.base.file import try_get_bluemira_path from bluemira.base.look_and_feel import ( @@ -30,6 +29,8 @@ if TYPE_CHECKING: from collections.abc import Iterator + import numpy.typing as npt + from bluemira.equilibria.equilibrium import Equilibrium from bluemira.equilibria.optimisation.problem import CoilsetOptimisationProblem from bluemira.equilibria.optimisation.problem.base import CoilsetOptimiserResult diff --git a/eudemo/config/build_config.json b/eudemo/config/build_config.json index 040a62b46c..92e47caf58 100644 --- a/eudemo/config/build_config.json +++ b/eudemo/config/build_config.json @@ -77,7 +77,8 @@ "nz": 127, "grid_scale_x": 1.6, "grid_scale_z": 1.7 - } + }, + "cocos": 3 }, "IVC": { "Wall silhouette": { diff --git a/eudemo/config/fixed_boundary_eqdsk.json b/eudemo/config/fixed_boundary_eqdsk.json index 691ae26726..15dca4fd86 100644 --- a/eudemo/config/fixed_boundary_eqdsk.json +++ b/eudemo/config/fixed_boundary_eqdsk.json @@ -2913,15 +2913,7 @@ -134.9306040434294 ] ], - "psibdry": [ - 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 - ], + "psibdry": 0.0, "psimag": 213.13029781275134, "xbdry": [ 6.0207, 6.022742199299811, 6.028883552470619, 6.039168849839246, diff --git a/eudemo/eudemo/equilibria/_designer.py b/eudemo/eudemo/equilibria/_designer.py index 2119d77617..5c157654b6 100644 --- a/eudemo/eudemo/equilibria/_designer.py +++ b/eudemo/eudemo/equilibria/_designer.py @@ -14,6 +14,8 @@ import matplotlib.pyplot as plt import numpy as np +from eqdsk import EQDSKInterface +from eqdsk.models import Sign from bluemira.base.designer import Designer from bluemira.base.file import get_bluemira_path, get_bluemira_root @@ -28,7 +30,6 @@ ) from bluemira.equilibria.fem_fixed_boundary.file import save_fixed_boundary_to_file from bluemira.equilibria.fem_fixed_boundary.utilities import get_mesh_boundary -from bluemira.equilibria.file import EQDSKInterface from bluemira.equilibria.optimisation.problem import ( UnconstrainedTikhonovCurrentGradientCOP, ) @@ -141,7 +142,11 @@ def run(self) -> Equilibrium: def read(self) -> Equilibrium: """Load an equilibrium from a file.""" - eq = Equilibrium.from_eqdsk(self.file_path) + eq = Equilibrium.from_eqdsk( + self.file_path, + from_cocos=self.build_config.get("cocos"), + qpsi_sign=self.build_config.get("qpsi", Sign.NEGATIVE), + ) self._update_params_from_eq(eq) return eq @@ -367,7 +372,12 @@ def read(self) -> tuple[Coordinates, CustomProfile]: """ Read in a fixed boundary equilibrium """ - data = EQDSKInterface.from_file(self.file_path) + data = EQDSKInterface.from_file( + self.file_path, + clockwise_phi=False, + volt_seconds_per_radian=True, + qpsi_sign=Sign.NEGATIVE, + ) lcfs_coords = Coordinates({"x": data.xbdry, "y": 0, "z": data.zbdry}) lcfs_coords.close() diff --git a/eudemo/eudemo/pf_coils/designer.py b/eudemo/eudemo/pf_coils/designer.py index 389b2f0b3a..c1470ef637 100644 --- a/eudemo/eudemo/pf_coils/designer.py +++ b/eudemo/eudemo/pf_coils/designer.py @@ -15,12 +15,12 @@ import matplotlib.pyplot as plt import numpy as np +from eqdsk import EQDSKInterface from bluemira.base.designer import Designer from bluemira.base.look_and_feel import bluemira_print from bluemira.base.parameter_frame import Parameter, ParameterFrame from bluemira.equilibria.coils import CoilSet -from bluemira.equilibria.file import EQDSKInterface from bluemira.equilibria.optimisation.constraints import ( CoilFieldConstraints, CoilForceConstraints, @@ -133,7 +133,7 @@ def read(self) -> CoilSet: # TODO: Load up equilibria from files and add states to manager - eqdsk = EQDSKInterface(**data[next(iter(data))]) + eqdsk = EQDSKInterface(**data[next(iter(data))], from_cocos=3) return CoilSet.from_group_vecs(eqdsk) def run(self) -> CoilSet: diff --git a/eudemo/eudemo_tests/equilibria/test_designer.py b/eudemo/eudemo_tests/equilibria/test_designer.py index e744572453..def11a88d2 100644 --- a/eudemo/eudemo_tests/equilibria/test_designer.py +++ b/eudemo/eudemo_tests/equilibria/test_designer.py @@ -8,6 +8,7 @@ from pathlib import Path import pytest +from eqdsk.models import Sign from bluemira.base.file import get_bluemira_path from bluemira.equilibria import Equilibrium @@ -43,12 +44,12 @@ def test_designer_converges_on_run(self): def test_designer_reads_file_in_read_mode(self): eqdsk = self.EQDSK_FILE designer = EquilibriumDesigner( - self.param_dict, {"run_mode": "read", "file_path": eqdsk} + self.param_dict, {"run_mode": "read", "file_path": eqdsk, "cocos": 3} ) eq = designer.execute() - ref_eq = Equilibrium.from_eqdsk(eqdsk) + ref_eq = Equilibrium.from_eqdsk(eqdsk, from_cocos=3, qpsi_sign=Sign.NEGATIVE) assert eq.analyse_plasma() == ref_eq.analyse_plasma() def test_ValueError_on_init_given_read_mode_and_no_file_path(self): diff --git a/eudemo/eudemo_tests/ivc/test_divertor_silhouette.py b/eudemo/eudemo_tests/ivc/test_divertor_silhouette.py index 14fa0ae3a6..ff7a6356f7 100644 --- a/eudemo/eudemo_tests/ivc/test_divertor_silhouette.py +++ b/eudemo/eudemo_tests/ivc/test_divertor_silhouette.py @@ -42,7 +42,7 @@ class TestDivertorSilhouetteDesigner: @classmethod def setup_class(cls): - cls.eq = Equilibrium.from_eqdsk(Path(DATA, "eqref_OOB.json")) + cls.eq = Equilibrium.from_eqdsk(Path(DATA, "eqref_OOB.json"), from_cocos=7) cls.separatrix = make_polygon(cls.eq.get_separatrix().xyz.T) _, cls.x_points = find_OX_points(cls.eq.x, cls.eq.z, cls.eq.psi()) diff --git a/eudemo/eudemo_tests/ivc/test_wall_silhouette.py b/eudemo/eudemo_tests/ivc/test_wall_silhouette.py index 501130eaf6..191ac033f5 100644 --- a/eudemo/eudemo_tests/ivc/test_wall_silhouette.py +++ b/eudemo/eudemo_tests/ivc/test_wall_silhouette.py @@ -55,7 +55,7 @@ class TestWallSilhouetteDesigner: @classmethod def setup_class(cls): - cls.eq = Equilibrium.from_eqdsk(Path(EQDATA, "eqref_OOB.json")) + cls.eq = Equilibrium.from_eqdsk(Path(EQDATA, "eqref_OOB.json"), from_cocos=7) _, cls.x_points = find_OX_points(cls.eq.x, cls.eq.z, cls.eq.psi()) def test_parameterisation_read(self): diff --git a/eudemo/eudemo_tests/test_tf_coils.py b/eudemo/eudemo_tests/test_tf_coils.py index ba7f3e38a4..8a87f4c001 100644 --- a/eudemo/eudemo_tests/test_tf_coils.py +++ b/eudemo/eudemo_tests/test_tf_coils.py @@ -56,7 +56,7 @@ class TestTFCoilDesigner: @classmethod def setup_class(cls): - cls.eq = Equilibrium.from_eqdsk(Path(EQDATA, "eqref_OOB.json")) + cls.eq = Equilibrium.from_eqdsk(Path(EQDATA, "eqref_OOB.json"), from_cocos=7) _, cls.x_points = find_OX_points(cls.eq.x, cls.eq.z, cls.eq.psi()) cls.lcfs = make_polygon(cls.eq.get_LCFS().xyz, closed=True) cls.vvts_koz = make_circle(10, center=(15, 0, 0), axis=(0.0, 1.0, 0.0)) diff --git a/pyproject.toml b/pyproject.toml index cf90acaae7..ce3f882596 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -21,6 +21,7 @@ dependencies = [ "click", "contourpy", "CoolProp", + "eqdsk>=0.4.0", "fortranformat", "gmsh", "imageio", diff --git a/requirements.txt b/requirements.txt index 6ab6f2c21d..e00fec2ebf 100644 --- a/requirements.txt +++ b/requirements.txt @@ -6,6 +6,7 @@ click==8.1.7 contourpy==1.2.1 CoolProp==6.6.0 cycler==0.12.1 +eqdsk==0.4.0 fonttools==4.53.0 fortranformat==2.0.0 gmsh==4.13.1 diff --git a/tests/equilibria/test_constraint_functions.py b/tests/equilibria/test_constraint_functions.py index 4dbc4c8fe4..73daf85070 100644 --- a/tests/equilibria/test_constraint_functions.py +++ b/tests/equilibria/test_constraint_functions.py @@ -68,8 +68,7 @@ def test_L2NormConstraint(self): class TestEquilibriumInput: @classmethod def setup_class(cls): - eq_name = "eqref_OOB.json" - cls.eq = Equilibrium.from_eqdsk(Path(TEST_PATH, eq_name)) + cls.eq = Equilibrium.from_eqdsk(Path(TEST_PATH, "eqref_OOB.json"), from_cocos=7) cls.coilset = cls.eq.coilset cls.vector = cls.coilset.current cls.scale = 1.0 diff --git a/tests/equilibria/test_constraints.py b/tests/equilibria/test_constraints.py index f741db7197..7272f4d162 100644 --- a/tests/equilibria/test_constraints.py +++ b/tests/equilibria/test_constraints.py @@ -6,6 +6,7 @@ import os +from eqdsk.models import Sign import numpy as np from bluemira.base.file import get_bluemira_path @@ -28,7 +29,7 @@ def test_constraint_weights(self): # Generate a test equilibrium path = get_bluemira_path("equilibria/test_data", subfolder="tests") fn = os.sep.join([path, "DN-DEMO_eqref.json"]) - eq = Equilibrium.from_eqdsk(fn) + eq = Equilibrium.from_eqdsk(fn, from_cocos=3, qpsi_sign=Sign.NEGATIVE) # Test that both default weights and custom weights can be applied for apply_weights in (True, False): diff --git a/tests/equilibria/test_equilibrium.py b/tests/equilibria/test_equilibrium.py index b52a301c00..d04a6cf117 100644 --- a/tests/equilibria/test_equilibrium.py +++ b/tests/equilibria/test_equilibrium.py @@ -10,12 +10,13 @@ import numpy as np import pytest +from eqdsk import EQDSKInterface +from eqdsk.models import Sign from matplotlib import pyplot as plt from bluemira.base.file import get_bluemira_path, try_get_bluemira_private_data_root from bluemira.equilibria.coils import CoilGroup, CoilSet from bluemira.equilibria.equilibrium import Equilibrium, FixedPlasmaEquilibrium -from bluemira.equilibria.file import EQDSKInterface from bluemira.equilibria.grid import Grid from bluemira.equilibria.optimisation.constraints import ( FieldNullConstraint, @@ -330,14 +331,18 @@ def test_betapliip_profile(self, shape): class TestEquilibrium: def test_double_null(self): path = get_bluemira_path("equilibria/test_data", subfolder="tests") - dn = Equilibrium.from_eqdsk(Path(path, "DN-DEMO_eqref.json")) + dn = Equilibrium.from_eqdsk( + Path(path, "DN-DEMO_eqref.json"), from_cocos=3, qpsi_sign=Sign.NEGATIVE + ) assert dn.is_double_null - sn = Equilibrium.from_eqdsk(Path(path, "eqref_OOB.json")) + sn = Equilibrium.from_eqdsk(Path(path, "eqref_OOB.json"), from_cocos=7) assert not sn.is_double_null def test_qpsi_calculation_modes(self): path = get_bluemira_path("equilibria/test_data", subfolder="tests") - dn = Equilibrium.from_eqdsk(Path(path, "DN-DEMO_eqref.json")) + dn = Equilibrium.from_eqdsk( + Path(path, "DN-DEMO_eqref.json"), from_cocos=3, qpsi_sign=Sign.NEGATIVE + ) with patch.object(dn, "q") as eq_q: res = dn.to_dict(qpsi_calcmode=0) assert eq_q.call_count == 0 @@ -355,7 +360,7 @@ def test_qpsi_calculation_modes(self): @pytest.mark.parametrize("grouping", [CoilSet, CoilGroup]) def test_woops_no_coils(self, grouping): testfile = Path(get_bluemira_path("eqdsk", subfolder="data"), "jetto.eqdsk_out") - e = EQDSKInterface.from_file(testfile) + e = EQDSKInterface.from_file(testfile, from_cocos=11) coil = grouping.from_group_vecs(e) assert isinstance(coil, grouping), "Check classmethod is making the right class" assert coil.current.any() == 0 @@ -370,7 +375,7 @@ def test_eq_coilnames(self): get_bluemira_path("equilibria/test_data", subfolder="tests"), "DN-DEMO_eqref_withCoilNames.json", ) - e = Equilibrium.from_eqdsk(testfile) + e = Equilibrium.from_eqdsk(testfile, from_cocos=3, qpsi_sign=Sign.NEGATIVE) assert e.coilset.name == [ *("PF_1", "PF_2", "PF_3", "PF_4", "PF_5", "PF_6"), *("CS_1", "CS_2", "CS_3", "CS_4", "CS_5"), @@ -388,7 +393,10 @@ def test_read_write(self, qpsi_calcmode, file_format): new_file_name = f"eqref_OOB_temp1.{file_format}" new_file_path = Path(data_path, new_file_name) - eq = Equilibrium.from_eqdsk(Path(data_path, file_name)) + eq = Equilibrium.from_eqdsk( + Path(data_path, file_name), from_cocos=7, to_cocos=None + ) + # Note we have recalculated the qpsi data here eq.to_eqdsk( directory=data_path, filename=new_file_name, @@ -397,7 +405,12 @@ def test_read_write(self, qpsi_calcmode, file_format): ) d1 = eq.to_dict(qpsi_calcmode=qpsi_calcmode) - eq2 = Equilibrium.from_eqdsk(new_file_path) + eq2 = Equilibrium.from_eqdsk( + new_file_path, + from_cocos=7 if qpsi_calcmode else 3, + to_cocos=None, + qpsi_sign=None if qpsi_calcmode else Sign.NEGATIVE, + ) d2 = eq2.to_dict(qpsi_calcmode=qpsi_calcmode) new_file_path.unlink() if file_format == "eqdsk": @@ -413,10 +426,10 @@ def setup_class(cls): root = try_get_bluemira_private_data_root() path = Path(root, "equilibria", "STEP_SPR_08") jetto_file = "SPR-008_3_Inputs_jetto.eqdsk_out" - jetto = EQDSKInterface.from_file(Path(path, jetto_file)) + jetto = EQDSKInterface.from_file(Path(path, jetto_file), from_cocos=11) cls.q_ref = jetto.qpsi eq_file = "SPR-008_3_Outputs_STEP_eqref.eqdsk" - cls.eq = Equilibrium.from_eqdsk(Path(path, eq_file)) + cls.eq = Equilibrium.from_eqdsk(Path(path, eq_file), from_cocos=7) def test_q_benchmark(self): n = len(self.q_ref) diff --git a/tests/equilibria/test_file.py b/tests/equilibria/test_file.py deleted file mode 100644 index 398d4effca..0000000000 --- a/tests/equilibria/test_file.py +++ /dev/null @@ -1,246 +0,0 @@ -# SPDX-FileCopyrightText: 2021-present M. Coleman, J. Cook, F. Franza -# SPDX-FileCopyrightText: 2021-present I.A. Maione, S. McIntosh -# SPDX-FileCopyrightText: 2021-present J. Morris, D. Short -# -# SPDX-License-Identifier: LGPL-2.1-or-later - -import copy -import json -from pathlib import Path -from typing import get_type_hints -from unittest import mock - -import fortranformat as ff -import numpy as np -import pytest -from typeguard import TypeCheckError, check_type - -from bluemira.base.file import get_bluemira_path -from bluemira.equilibria.file import EQDSKInterface -from bluemira.utilities.tools import compare_dicts -from tests._helpers import combine_text_mock_write_calls - -OPEN = "builtins.open" - - -class TestEQDSKInterface: - path = get_bluemira_path("equilibria/test_data", subfolder="tests") - testfiles = ( - Path(get_bluemira_path("eqdsk", subfolder="data"), "jetto.eqdsk_out"), - Path(path, "DN-DEMO_eqref.json"), - Path(path, "eqref_OOB.json"), - Path(path, "DN-DEMO_eqref_withCoilNames.json"), - ) - - @classmethod - def setup_class(cls): - data_dir = get_bluemira_path("equilibria", subfolder="data") - data_file = Path(data_dir, "DN-DEMO_eqref.json") - with open(data_file) as f: - cls.eudemo_sof_data = json.load(f) - - def read_strict_geqdsk(self, file_path): - """ - Reads an input EQDSK file in, assuming strict adherence to the - GEQDSK format. Used to check bluemira outputs can be read by - external readers. - - Note: The main bluemira GEQDSK reader is more forgiving to - format variations than this! - - Parameters - ---------- - file_path: str - Full path string of the file - - """ - - # Create FortranRecordReader objects with the Fortran format - # edit descriptors to be used to parse the G-EQDSK input. - f2000 = ff.FortranRecordReader("a48,3i4") - f2020 = ff.FortranRecordReader("5e16.9") - f2022 = ff.FortranRecordReader("2i5") - fCSTM = ff.FortranRecordReader("i5") - - # Define helper function to read in flattened arrays. - def read_flat_array(fortran_format, array_size): - """ - Reads in a flat (1D) numpy array from a G-EQDSK file. - - Parameters - ---------- - fortran_format: ff.FortranRecordReader - FortranRecordReader object for Fortran format edit descriptor - to be used to parse the format of each line of the output. - array_size: int - Number of elements in array to be read in. - - Returns - ------- - array: np.array - 1D Numpy array of length array_size populated by elements from - the GEQDSK file input. - """ - - # Initialise numpy array and read in first line. - array = np.zeros((array_size,)) - line = fortran_format.read(file.readline()) - # Define a counter to track which column in a line - # is currently being saved into the array. - col = 0 - # Populate array. If the column index moves past the - # end of a line, it is reset to zero and the next line is read. - for i in range(array_size): - if col == len(line): - line = fortran_format.read(file.readline()) - col = 0 - array[i] = line[col] - col += 1 - return array - - # Open file. - with open(file_path) as file: - # Read in data. Variable names are for readability; - # strict format is as defined at - # https://w3.pppl.gov/ntcc/TORAY/G_EQDSK.pdf - _id, _, nx, nz = f2000.read(file.readline()) - xdim, zdim, xcentre, xgrid1, zmid = f2020.read(file.readline()) - xmag, zmag, psimag, psibdry, bcentre = f2020.read(file.readline()) - cplasma, psimag, _, xmag, _ = f2020.read(file.readline()) - zmag, _, psibdry, _, _ = f2020.read(file.readline()) - fpol = read_flat_array(f2020, nx) - pressure = read_flat_array(f2020, nx) - ffprime = read_flat_array(f2020, nx) - pprime = read_flat_array(f2020, nx) - psi = read_flat_array(f2020, nx * nz) - qpsi = read_flat_array(f2020, nx) - nbdry, nlim = f2022.read(file.readline()) - xbdry_zbdry = read_flat_array(f2020, 2 * nbdry) - xlim_zlim = read_flat_array(f2020, 2 * nlim) - - # Read in coil information, as found in the GEQDSK extension - # used by bluemira. - (ncoil,) = fCSTM.read(file.readline()) - coil = read_flat_array(f2020, 5 * ncoil) - - @pytest.mark.parametrize("file", testfiles) - def test_read(self, file): - # Create EQDSK file interface and read data to a dict - eqdsk = EQDSKInterface.from_file(file) - d1 = eqdsk.to_dict() - - # Write data read in from test file into a new EQDSK - # file, with the suffix "_temp" - name = Path(file).stem + "_temp" - fname = Path(self.path, f"{name}.eqdsk") - eqdsk.write(fname, file_format="eqdsk") - d2 = eqdsk.to_dict() - - # Check eqdsk is readable by Fortran readers. - # This demands stricter adherence to the G-EQDSK - # format than bluemira's main reader. - self.read_strict_geqdsk(fname) - - # Write data read in from test file into a new JSON - # file, with the suffix "_temp" - jname = fname.with_suffix("").with_suffix(".json") - eqdsk.write(jname, file_format="json") - d3 = EQDSKInterface.from_file(jname).to_dict() - - # Clean up temporary files - fname.unlink() - jname.unlink() - - # Compare dictionaries to check data hasn't - # been changed. - assert compare_dicts(d1, d2, verbose=True) - assert compare_dicts(d1, d3, verbose=True) - assert compare_dicts(d2, d3, verbose=True) - - def test_read_matches_values_in_file(self): - eq = EQDSKInterface.from_file(self.testfiles[0]) - - assert eq.nz == 151 - assert eq.nx == 151 - assert eq.xdim == pytest.approx(3.14981545) - assert eq.ncoil == 0 - assert eq.xc.size == 0 - assert eq.nbdry == 72 - np.testing.assert_allclose( - eq.xbdry[:3], [0.399993127e01, 0.399150254e01, 0.396906908e01] - ) - np.testing.assert_allclose( - eq.zbdry[-3:], [-0.507187454e00, -0.240712636e00, 0.263892047e-01] - ) - - def test_values_match_annotated_types(self): - eq = EQDSKInterface.from_file(self.testfiles[0]) - - mismatched = [] - for key, value_type in get_type_hints(EQDSKInterface).items(): - value = getattr(eq, key) - try: - check_type(value, value_type) - except TypeCheckError: - mismatched.append((key, type(value), value_type)) - assert not mismatched - - def test_write_then_read_in_json_format(self): - # Test with jetto_eqdsk.out (no coil information) - eq = EQDSKInterface.from_file(self.testfiles[0]) - - with mock.patch(OPEN, new_callable=mock.mock_open) as open_mock: - eq.write("some/path.json", file_format="json") - written = combine_text_mock_write_calls(open_mock) - - with mock.patch( - OPEN, new_callable=mock.mock_open, read_data=written - ) as open_mock: - eq2 = EQDSKInterface.from_file("/some/path.json") - - assert eq2.nz == 151 - assert eq2.nbdry == 72 - assert eq2.coil_names is None - assert eq2.coil_types is None - - # Test with DN-DEMO_eqref_withCoilNames.json (with coil information) - eq3 = EQDSKInterface.from_file(self.testfiles[3]) - - with open(str(self.testfiles[3])) as data_file: - eq3_json_dict = json.load(data_file) - - with mock.patch(OPEN, new_callable=mock.mock_open) as open_mock: - eq3.write("some/path.json", file_format="json") - written = combine_text_mock_write_calls(open_mock) - - # check the contents of the written mock variable - compare_dicts(eq3_json_dict, json.loads(written)) - - with mock.patch( - OPEN, new_callable=mock.mock_open, read_data=written - ) as open_mock: - eq4 = EQDSKInterface.from_file("/some/path.json") - - assert eq4.ncoil == 11 - assert eq4.coil_names[2] == "PF_3" - assert eq4.coil_names[10] == "CS_5" - assert eq4.coil_types[1] == "PF" - assert eq4.coil_types[7] == "CS" - - def test_derived_field_is_calculated_if_not_given(self): - data = copy.deepcopy(self.eudemo_sof_data) - for field in ["x", "z", "psinorm"]: - del data[field] - - open_mock = mock.mock_open(read_data=json.dumps(data)) - with mock.patch(OPEN, new=open_mock): - eqdsk = EQDSKInterface.from_file("/some/file.json") - - np.testing.assert_allclose(eqdsk.x, self.eudemo_sof_data["x"]) - np.testing.assert_allclose(eqdsk.z, self.eudemo_sof_data["z"]) - # The calculation used for psinorm has changed since the - # eudemo_sof_data was created - so we can't compare to that in - # this case. - np.testing.assert_allclose( - eqdsk.psinorm, np.linspace(0, 1, len(self.eudemo_sof_data["fpol"])) - ) diff --git a/tests/equilibria/test_find.py b/tests/equilibria/test_find.py index 91b4806cd4..566aec850e 100644 --- a/tests/equilibria/test_find.py +++ b/tests/equilibria/test_find.py @@ -10,6 +10,7 @@ import numpy as np import pytest +from eqdsk.models import Sign from bluemira.base.file import get_bluemira_path from bluemira.equilibria.equilibrium import Equilibrium @@ -61,7 +62,7 @@ def test_inv_2x2_jacobian(): class TestFindLCFSSeparatrix: def test_other_grid(self): - sof = Equilibrium.from_eqdsk(Path(DATA, "eqref_OOB.json")) + sof = Equilibrium.from_eqdsk(Path(DATA, "eqref_OOB.json"), from_cocos=7) psi = sof.psi() o_points, x_points = sof.get_OX_points(psi) grid_tol = np.hypot(sof.grid.dx, sof.grid.dz) @@ -83,7 +84,9 @@ def test_other_grid(self): assert np.amin(distances) <= grid_tol def test_double_null(self): - sof = Equilibrium.from_eqdsk(Path(DATA, "DN-DEMO_eqref.json")) + sof = Equilibrium.from_eqdsk( + Path(DATA, "DN-DEMO_eqref.json"), from_cocos=3, qpsi_sign=Sign.NEGATIVE + ) psi = sof.psi() o_points, x_points = sof.get_OX_points(psi) grid_tol = np.hypot(sof.grid.dx, sof.grid.dz) @@ -126,8 +129,10 @@ def test_recursion(self): class TestGetLegs: @classmethod def setup_class(cls): - cls.sn_eq = Equilibrium.from_eqdsk(Path(DATA, "eqref_OOB.json")) - cls.dn_eq = Equilibrium.from_eqdsk(Path(DATA, "DN-DEMO_eqref.json")) + cls.sn_eq = Equilibrium.from_eqdsk(Path(DATA, "eqref_OOB.json"), from_cocos=7) + cls.dn_eq = Equilibrium.from_eqdsk( + Path(DATA, "DN-DEMO_eqref.json"), from_cocos=3, qpsi_sign=Sign.NEGATIVE + ) cls.falsified_dn_eq = deepcopy(cls.sn_eq) def test_legflux(self): diff --git a/tests/equilibria/test_flux_surfaces.py b/tests/equilibria/test_flux_surfaces.py index 65b4b7a6c1..b2817cb93d 100644 --- a/tests/equilibria/test_flux_surfaces.py +++ b/tests/equilibria/test_flux_surfaces.py @@ -9,6 +9,7 @@ import numpy as np import pytest +from eqdsk.models import Sign from bluemira.base.constants import EPS from bluemira.base.file import get_bluemira_path @@ -38,7 +39,7 @@ class TestOpenFluxSurfaceStuff: @classmethod def setup_class(cls): eq_name = "eqref_OOB.json" - cls.eq = Equilibrium.from_eqdsk(Path(TEST_PATH, eq_name)) + cls.eq = Equilibrium.from_eqdsk(Path(TEST_PATH, eq_name), from_cocos=7) def test_bad_geometry(self): closed_coords = Coordinates({"x": [0, 4, 5, 8, 0], "z": [1, 2, 3, 4, 1]}) @@ -145,7 +146,7 @@ class TestFieldLine: @classmethod def setup_class(cls): eq_name = "eqref_OOB.json" - cls.eq = Equilibrium.from_eqdsk(Path(TEST_PATH, eq_name)) + cls.eq = Equilibrium.from_eqdsk(Path(TEST_PATH, eq_name), from_cocos=7) cls.flt = FieldLineTracer(cls.eq) cls.field_line = cls.flt.trace_field_line(13, 0, n_points=1000) @@ -219,7 +220,9 @@ def _extract_endpoint(self, field_line): def test_poloidal_angle(): eq_name = "DN-DEMO_eqref.json" - eq = Equilibrium.from_eqdsk(Path(TEST_PATH, eq_name)) + eq = Equilibrium.from_eqdsk( + Path(TEST_PATH, eq_name), from_cocos=3, qpsi_sign=Sign.NEGATIVE + ) # Building inputs x_strike = 10.0 z_strike = -7.5 diff --git a/tests/equilibria/test_harmonics.py b/tests/equilibria/test_harmonics.py index 4ef4f95e15..38365c847b 100644 --- a/tests/equilibria/test_harmonics.py +++ b/tests/equilibria/test_harmonics.py @@ -9,6 +9,7 @@ import numpy as np import pytest +from eqdsk.models import Sign from bluemira.base.constants import EPS from bluemira.base.file import get_bluemira_path @@ -144,7 +145,11 @@ def test_collocation_points(): def test_coils_outside_sphere_vacuum_psi(): - eq = Equilibrium.from_eqdsk(Path(TEST_PATH, "SH_test_file.json").as_posix()) + eq = Equilibrium.from_eqdsk( + Path(TEST_PATH, "SH_test_file.json").as_posix(), + from_cocos=3, + qpsi_sign=Sign.NEGATIVE, + ) sh_coil_names, bdry_r = coils_outside_lcfs_sphere(eq) assert len(sh_coil_names) == 16 @@ -165,7 +170,11 @@ def test_coils_outside_sphere_vacuum_psi(): def test_get_psi_harmonic_amplitudes(): - eq = Equilibrium.from_eqdsk(Path(TEST_PATH, "SH_test_file.json").as_posix()) + eq = Equilibrium.from_eqdsk( + Path(TEST_PATH, "SH_test_file.json").as_posix(), + from_cocos=3, + qpsi_sign=Sign.NEGATIVE, + ) test_colocation = collocation_points( plasma_boundary=eq.get_LCFS(), @@ -206,7 +215,11 @@ def test_get_psi_harmonic_amplitudes(): def test_spherical_harmonic_approximation(): - eq = Equilibrium.from_eqdsk(Path(TEST_PATH, "SH_test_file.json").as_posix()) + eq = Equilibrium.from_eqdsk( + Path(TEST_PATH, "SH_test_file.json").as_posix(), + from_cocos=3, + qpsi_sign=Sign.NEGATIVE, + ) ( _, @@ -314,7 +327,11 @@ def test_SphericalHarmonicConstraintFunction(): def test_SphericalHarmonicConstraint(): - eq = Equilibrium.from_eqdsk(Path(TEST_PATH, "SH_test_file.json").as_posix()) + eq = Equilibrium.from_eqdsk( + Path(TEST_PATH, "SH_test_file.json").as_posix(), + from_cocos=3, + qpsi_sign=Sign.NEGATIVE, + ) sh_coil_names, _ = coils_outside_lcfs_sphere(eq) ref_harmonics = np.array([ diff --git a/tests/equilibria/test_optimiser.py b/tests/equilibria/test_optimiser.py index 2fa4739f2d..d7c314ca80 100644 --- a/tests/equilibria/test_optimiser.py +++ b/tests/equilibria/test_optimiser.py @@ -8,6 +8,7 @@ from pathlib import Path import numpy as np +from eqdsk.models import Sign from bluemira.base.file import get_bluemira_path from bluemira.equilibria.coils import ( @@ -144,7 +145,9 @@ def setup_class(cls): Path( get_bluemira_path("equilibria/test_data", subfolder="tests"), "DN-DEMO_eqref.json", - ).as_posix() + ).as_posix(), + from_cocos=3, + qpsi_sign=Sign.NEGATIVE, ) _dummy_eq_partial = deepcopy(_dummy_eq_none) _dummy_eq_sym = deepcopy(_dummy_eq_partial) diff --git a/tests/equilibria/test_st_equilibrium.py b/tests/equilibria/test_st_equilibrium.py index d0db7c4e87..d7d0da68c0 100644 --- a/tests/equilibria/test_st_equilibrium.py +++ b/tests/equilibria/test_st_equilibrium.py @@ -14,6 +14,8 @@ import numpy as np import pytest +from eqdsk import EQDSKInterface +from eqdsk.models import Sign from bluemira.base.file import get_bluemira_root from bluemira.equilibria import ( @@ -25,7 +27,6 @@ PicardIterator, SymmetricCircuit, ) -from bluemira.equilibria.file import EQDSKInterface from bluemira.equilibria.optimisation.constraints import ( IsofluxConstraint, MagneticConstraintSet, @@ -46,11 +47,13 @@ def setup_class(cls): private = os.path.split(root)[0] private = Path(private, "bluemira-private-data/equilibria/STEP_SPR_08") eq_name = "STEP_SPR08_BLUEPRINT.json" - cls.eq_blueprint = Equilibrium.from_eqdsk(Path(private, eq_name)) + cls.eq_blueprint = Equilibrium.from_eqdsk( + Path(private, eq_name), from_cocos=3, to_cocos=11, qpsi_sign=Sign.NEGATIVE + ) jeq_name = "jetto.eqdsk_out" filename = Path(private, jeq_name) - cls.profiles = CustomProfile.from_eqdsk(filename) - cls.jeq_dict = EQDSKInterface.from_file(filename) + cls.jeq_dict = EQDSKInterface.from_file(filename, from_cocos=11) + cls.profiles = CustomProfile.from_eqdsk(cls.jeq_dict) def test_equilibrium(self): build_tweaks = { diff --git a/tests/radiation_transport/test_advective_transport.py b/tests/radiation_transport/test_advective_transport.py index 0286e1aca9..96f61a4e0f 100644 --- a/tests/radiation_transport/test_advective_transport.py +++ b/tests/radiation_transport/test_advective_transport.py @@ -9,6 +9,7 @@ import numpy as np import pytest +from eqdsk.models import Sign import bluemira.radiation_transport.flux_surfaces_maker as fsm from bluemira.base.file import get_bluemira_path @@ -73,7 +74,7 @@ class TestChargedParticleRecursionSN: def setup_class(cls): eq_name = "EU-DEMO_EOF.json" filename = Path(EQ_PATH, eq_name) - eq = Equilibrium.from_eqdsk(filename) + eq = Equilibrium.from_eqdsk(filename, from_cocos=3, qpsi_sign=Sign.NEGATIVE) fw_name = "first_wall.json" filename = Path(TEST_PATH, fw_name) fw = Coordinates.from_json(filename) @@ -199,7 +200,7 @@ class TestChargedParticleRecursionDN: def setup_class(cls): eq_name = "DN-DEMO_eqref.json" filename = Path(EQ_PATH, eq_name) - eq = Equilibrium.from_eqdsk(filename) + eq = Equilibrium.from_eqdsk(filename, from_cocos=3, qpsi_sign=Sign.NEGATIVE) fw_name = "DN_fw_shape.json" filename = Path(TEST_PATH, fw_name) fw = Coordinates.from_json(filename) diff --git a/tests/radiation_transport/test_flux_surface_maker.py b/tests/radiation_transport/test_flux_surface_maker.py index b9e899a8cb..827d219c3d 100644 --- a/tests/radiation_transport/test_flux_surface_maker.py +++ b/tests/radiation_transport/test_flux_surface_maker.py @@ -7,6 +7,7 @@ from pathlib import Path import pytest +from eqdsk.models import Sign from bluemira.base.file import get_bluemira_path from bluemira.equilibria.equilibrium import Equilibrium @@ -67,7 +68,9 @@ class TestMakeFS: @classmethod def setup_class(cls): eq_name = "DN-DEMO_eqref.json" - cls.eq = Equilibrium.from_eqdsk(Path(TEST_PATH, eq_name)) + cls.eq = Equilibrium.from_eqdsk( + Path(TEST_PATH, eq_name), from_cocos=3, qpsi_sign=Sign.NEGATIVE + ) cls.o_point = cls.eq.get_OX_points()[0][0] # 1st o_point cls.yz_plane = BluemiraPlane.from_3_points( [0, 0, cls.o_point.z], [1, 0, cls.o_point.z], [1, 1, cls.o_point.z] diff --git a/tests/radiation_transport/test_radiation_profile.py b/tests/radiation_transport/test_radiation_profile.py index d81f0bd69e..515f6e3871 100644 --- a/tests/radiation_transport/test_radiation_profile.py +++ b/tests/radiation_transport/test_radiation_profile.py @@ -8,6 +8,7 @@ import numpy as np import pytest +from eqdsk.models import Sign from bluemira.base import constants from bluemira.base.constants import raw_uc @@ -45,7 +46,7 @@ class TestCoreRadiation: def setup_class(cls): eq_name = "DN-DEMO_eqref.json" filename = Path(EQ_PATH, eq_name) - eq = Equilibrium.from_eqdsk(filename) + eq = Equilibrium.from_eqdsk(filename, from_cocos=3, qpsi_sign=Sign.NEGATIVE) fw_name = "DN_fw_shape.json" filename = Path(TEST_PATH, fw_name) fw_shape = Coordinates.from_json(filename)