From dd4c1313e87b59e49f8b016a9a5d99b46c27cfd6 Mon Sep 17 00:00:00 2001 From: MatejTomes Date: Mon, 7 Jan 2019 16:31:11 +0100 Subject: [PATCH 1/6] files cleanup --- cherab/compass/__init__.py | 4 +- cherab/compass/plasma.py | 96 -------------------------------------- 2 files changed, 2 insertions(+), 98 deletions(-) delete mode 100644 cherab/compass/plasma.py diff --git a/cherab/compass/__init__.py b/cherab/compass/__init__.py index 63b69e8..1449d52 100644 --- a/cherab/compass/__init__.py +++ b/cherab/compass/__init__.py @@ -1,2 +1,2 @@ - -from .plasma import load_core_plasma_from_files +from .equilibrium import * +#from .plasma import load_core_plasma_from_files diff --git a/cherab/compass/plasma.py b/cherab/compass/plasma.py deleted file mode 100644 index 68fb30f..0000000 --- a/cherab/compass/plasma.py +++ /dev/null @@ -1,96 +0,0 @@ - -# Core and external imports -import numpy as np -from scipy.constants import electron_mass, atomic_mass -from raysect.core import Vector3D, translate -from raysect.primitive import Cylinder - -# Cherab and raysect imports -from cherab.core import Species, Plasma, Maxwellian -from cherab.core.atomic import carbon -from cherab.core.math import Interpolate1DCubic, Interpolate2DCubic, IsoMapper3D, AxisymmetricMapper, Constant3D, Blend3D - - -def load_core_plasma_from_files(world, atomic_data, psirz_file, psin_file, electron_file, carbon_file): - - # Load equilibrium psi function - loaded = np.loadtxt(psirz_file, delimiter="\t", skiprows=1) - psi_R = loaded[:, 0] - psi_Z = loaded[:, 1] - psi_map = np.loadtxt(psin_file, delimiter="\t", skiprows=1) - psin_2d = Interpolate2DCubic(psi_R, psi_Z, psi_map, extrapolate=True) - psin_3d = AxisymmetricMapper(psin_2d) - - # Load electron species - loaded = np.loadtxt(electron_file, delimiter="\t", skiprows=1) - psi_coord = loaded[:, 0] - electrons_temperature_axis = loaded[:, 1] - electrons_densiy_axis = loaded[:, 2] - - flow_velocity = lambda x, y, z: Vector3D(0, 0, 0) # Ignoring velocity, necessary but not used for this model. - - electorn_temperature_vs_psi = Interpolate1DCubic(psi_coord, electrons_temperature_axis, extrapolate=True) - electron_temperature = Blend3D(Constant3D(0.0), IsoMapper3D(psin_3d, electorn_temperature_vs_psi), _inside_lcfs) - electron_density_vs_psi = Interpolate1DCubic(psi_coord, electrons_densiy_axis, extrapolate=True) - electron_density = Blend3D(Constant3D(0.0), IsoMapper3D(psin_3d, electron_density_vs_psi), _inside_lcfs) - electron_distribution = Maxwellian(electron_density, electron_temperature, flow_velocity, electron_mass) - - # Load carbon species - c_abundance_data = np.loadtxt(carbon_file, delimiter="\t", skiprows=1) - - c0_density_psi = Interpolate1DCubic(psi_coord, c_abundance_data[:, 0], extrapolate=True) - c0_density = Blend3D(Constant3D(0.0), IsoMapper3D(psin_3d, c0_density_psi), _inside_lcfs) - c0_distribution = Maxwellian(c0_density, electron_temperature, flow_velocity, carbon.atomic_weight * atomic_mass) - c0_species = Species(carbon, 0, c0_distribution) - - c1_density_psi = Interpolate1DCubic(psi_coord, c_abundance_data[:, 1], extrapolate=True) - c1_density = Blend3D(Constant3D(0.0), IsoMapper3D(psin_3d, c1_density_psi), _inside_lcfs) - c1_distribution = Maxwellian(c1_density, electron_temperature, flow_velocity, carbon.atomic_weight * atomic_mass) - c1_species = Species(carbon, 1, c1_distribution) - - c2_density_psi = Interpolate1DCubic(psi_coord, c_abundance_data[:, 2], extrapolate=True) - c2_density = Blend3D(Constant3D(0.0), IsoMapper3D(psin_3d, c2_density_psi), _inside_lcfs) - c2_distribution = Maxwellian(c2_density, electron_temperature, flow_velocity, carbon.atomic_weight * atomic_mass) - c2_species = Species(carbon, 2, c2_distribution) - - c3_density_psi = Interpolate1DCubic(psi_coord, c_abundance_data[:, 3], extrapolate=True) - c3_density = Blend3D(Constant3D(0.0), IsoMapper3D(psin_3d, c3_density_psi), _inside_lcfs) - c3_distribution = Maxwellian(c3_density, electron_temperature, flow_velocity, carbon.atomic_weight * atomic_mass) - c3_species = Species(carbon, 3, c3_distribution) - - c4_density_psi = Interpolate1DCubic(psi_coord, c_abundance_data[:, 4], extrapolate=True) - c4_density = Blend3D(Constant3D(0.0), IsoMapper3D(psin_3d, c4_density_psi), _inside_lcfs) - c4_distribution = Maxwellian(c4_density, electron_temperature, flow_velocity, carbon.atomic_weight * atomic_mass) - c4_species = Species(carbon, 4, c4_distribution) - - c5_density_psi = Interpolate1DCubic(psi_coord, c_abundance_data[:, 5], extrapolate=True) - c5_density = Blend3D(Constant3D(0.0), IsoMapper3D(psin_3d, c5_density_psi), _inside_lcfs) - c5_distribution = Maxwellian(c5_density, electron_temperature, flow_velocity, carbon.atomic_weight * atomic_mass) - c5_species = Species(carbon, 5, c5_distribution) - - c6_density_psi = Interpolate1DCubic(psi_coord, c_abundance_data[:, 6], extrapolate=True) - c6_density = Blend3D(Constant3D(0.0), IsoMapper3D(psin_3d, c6_density_psi), _inside_lcfs) - c6_distribution = Maxwellian(c6_density, electron_temperature, flow_velocity, carbon.atomic_weight * atomic_mass) - c6_species = Species(carbon, 6, c6_distribution) - - # Assemble the plasma object - plasma = Plasma(world) - plasma.atomic_data = atomic_data - plasma.electron_distribution = electron_distribution - plasma.composition = [c0_species, c1_species, c2_species, c3_species, c4_species, c5_species, c6_species] - plasma.b_field = lambda x, y, z: Vector3D(0, 0, 0.) # ignoring magnetic field, necessary but not used for this model. - - # Give the plasma some information about COMPASS' cylindrical geometry - plasma.geometry = Cylinder(1.0, 1.0) - plasma.geometry_transform = translate(0, 0, -0.5) - - return plasma - - -# Criteria for being inside the plasma -def _inside_lcfs(x, y, z): - if np.sqrt(x**2 + y**2) < 0.35 or np.sqrt(x**2 + y**2) > 0.75: - return 0 - else: - return 1.0 - From 1c277e12c9dffb54a22e8c854c6d7cd0209a9019 Mon Sep 17 00:00:00 2001 From: MatejTomes Date: Mon, 3 Jun 2019 17:35:46 +0200 Subject: [PATCH 2/6] add q_profile as reaction to cherab evolution loading of compass efit file should be working now. It was needed to react to cherab changes. --- cherab/compass/equilibrium/equilibrium.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/cherab/compass/equilibrium/equilibrium.py b/cherab/compass/equilibrium/equilibrium.py index 9ef3419..aaa31b6 100644 --- a/cherab/compass/equilibrium/equilibrium.py +++ b/cherab/compass/equilibrium/equilibrium.py @@ -75,7 +75,8 @@ def _read_efitfile(self, path): self._xpoint_z = f["output/separatrixGeometry/xpointCoordsZ"].value - self._rBphi = f["output/fluxFunctionProfiles/rBphi"].value + self._f = f["output/fluxFunctionProfiles/rBphi"].value + self._q = f["output/fluxFunctionProfiles/q"].value self._psin = f["output/fluxFunctionProfiles/normalizedPoloidalFlux"].value @@ -114,12 +115,14 @@ def time(self, time): xpoints = self._process_efit_points(self._xpoint_r[time_index, :], self._xpoint_z[time_index, :]) - f_profile = np.array([self._psin, self._rBphi[time_index, :]]) + f_profile = np.array([self._psin, self._f[time_index, :]]) + + q_profile = np.array([self._psin, self._q[time_index, :]]) equilibrium = EFITEquilibrium(r=r, z=z, psi_grid=psi_grid, psi_axis=psi_axis, psi_lcfs=psi_lcfs, magnetic_axis=magnetic_axis, x_points = xpoints, strike_points=strikepoints, - f_profile=f_profile, b_vacuum_radius=b_vacuum_radius, + f_profile=f_profile, q_profile=q_profile, b_vacuum_radius=b_vacuum_radius, b_vacuum_magnitude=b_vacuum_magnitude, lcfs_polygon=lcfs_polygon, limiter_polygon=limiter_polygon, time=time_slice) From ce00829bab3d326b211c8cccbad8f0ba3a6a6733 Mon Sep 17 00:00:00 2001 From: MatejTomes Date: Fri, 25 Oct 2019 18:42:57 +0200 Subject: [PATCH 3/6] rewrite of COMPASSEquilibrium --- cherab/__init__.py | 1 - cherab/compass/equilibrium/__init__.py | 1 + cherab/compass/equilibrium/compass.py | 54 ++++++ cherab/compass/equilibrium/compass_upgrade.py | 13 ++ cherab/compass/equilibrium/equilibrium.py | 174 ------------------ .../compass/equilibrium/equilibrium_base.py | 88 +++++++++ 6 files changed, 156 insertions(+), 175 deletions(-) create mode 100644 cherab/compass/equilibrium/compass.py create mode 100644 cherab/compass/equilibrium/compass_upgrade.py delete mode 100644 cherab/compass/equilibrium/equilibrium.py create mode 100644 cherab/compass/equilibrium/equilibrium_base.py diff --git a/cherab/__init__.py b/cherab/__init__.py index 7d29d19..de40ea7 100644 --- a/cherab/__init__.py +++ b/cherab/__init__.py @@ -1,2 +1 @@ - __import__('pkg_resources').declare_namespace(__name__) diff --git a/cherab/compass/equilibrium/__init__.py b/cherab/compass/equilibrium/__init__.py index e69de29..e0d5f0d 100644 --- a/cherab/compass/equilibrium/__init__.py +++ b/cherab/compass/equilibrium/__init__.py @@ -0,0 +1 @@ +from .compass import COMPASSEquilibrium \ No newline at end of file diff --git a/cherab/compass/equilibrium/compass.py b/cherab/compass/equilibrium/compass.py new file mode 100644 index 0000000..0fcc495 --- /dev/null +++ b/cherab/compass/equilibrium/compass.py @@ -0,0 +1,54 @@ +from cherab.compass.equilibrium.equilibrium_base import COMPASSEquilibrium_base + +class COMPASSEquilibrium(COMPASSEquilibrium_base): + + @classmethod + def from_cdb(cls, shot_number): + """ + Obtains EFIT equilibrium data from CDB database for the specified shot. + :param shot_number: COMPASS shot number + :return: COMPASSEquilibrium class + """ + import cdb_extras.xarray_support as cdbxr #if we want to read equilibrium data from a file without cdb installed.... + + shot = cdbxr.Shot(shot_number) + + equilibrium = cls() + equilibrium.shot_number = shot_number + + + equilibrium._data = shot["psi_RZ/EFIT"].copy() + equilibrium._data.name = "psi_grid" + equilibrium._data = equilibrium._data.to_dataset() + + equilibrium._data = equilibrium._data.merge(shot["psi_mag_axis/EFIT"].to_dataset()) + equilibrium._data = equilibrium._data.rename({"psi_mag_axis":"psi_axis"}) + + equilibrium._data = equilibrium._data.merge(shot["psi_lcfs/EFIT"].to_dataset()) + + equilibrium._data = equilibrium._data.merge(shot["RZ_mag_axis/EFIT"]) + equilibrium._data = equilibrium._data.rename({"R_mag_axis": "R_magnetic_axis", "Z_mag_axis": "Z_magnetic_axis"}) + + equilibrium._data = equilibrium._data.merge(shot["RZ_xpoint/EFIT"]) + + + equilibrium._data = equilibrium._data.merge(shot["RZ_strike_points/EFIT"]) + + equilibrium._data = equilibrium._data.merge(shot["RBphi/EFIT"].to_dataset()) + equilibrium._data = equilibrium._data.rename({"RBphi": "f_profile"}) + + equilibrium._data = equilibrium._data.merge(shot["q/EFIT"].to_dataset()) + equilibrium._data = equilibrium._data.rename({"q": "q_profile"}) + + equilibrium._data = equilibrium._data.merge(shot["B_vac_R_geom/EFIT"].to_dataset()) + equilibrium._data = equilibrium._data.rename({"B_vac_R_geom": "b_vacuum_magnitude"}) + equilibrium._data.attrs["b_vacuum_radius"] = 0.56 + + equilibrium._data = equilibrium._data.merge(shot["RZ_bound/EFIT"]) + equilibrium._data = equilibrium._data.rename({"R_bound": "R_lcfs", "Z_bound": "Z_lcfs" }) + + equilibrium._data = equilibrium._data.merge(shot["R_limiter_input/EFIT"].to_dataset()) + equilibrium._data = equilibrium._data.merge(shot["Z_limiter_input/EFIT"].to_dataset()) + equilibrium._data = equilibrium._data.rename({"R_limiter_input": "R_limiter", "Z_limiter_input": "Z_limiter" }) + + return equilibrium \ No newline at end of file diff --git a/cherab/compass/equilibrium/compass_upgrade.py b/cherab/compass/equilibrium/compass_upgrade.py new file mode 100644 index 0000000..341047f --- /dev/null +++ b/cherab/compass/equilibrium/compass_upgrade.py @@ -0,0 +1,13 @@ +from cherab.compass.equilibrium.equilibrium_base import COMPASSEquilibrium_base +import cdb_extras.xarray_support as cdbxr +from pyCDB import client +cdb = client.CDBClient(host="cudb.tok.ipp.cas.cz",data_root="/compass/ +CC19_COMPASS-U_data/") +#class COMPASSUpgradeEquilibrium(COMPASSEquilibrium_base): + +# def from_cdb(self, shot_number): + +shot_number = 6400 + +shot = cudbxr.Shot(shot_number) +shot["Bphi/Fiesta_OUT"] \ No newline at end of file diff --git a/cherab/compass/equilibrium/equilibrium.py b/cherab/compass/equilibrium/equilibrium.py deleted file mode 100644 index 9ef3419..0000000 --- a/cherab/compass/equilibrium/equilibrium.py +++ /dev/null @@ -1,174 +0,0 @@ -import h5py -import numpy as np -from raysect.core import Point2D -from cherab.tools.equilibrium import EFITEquilibrium - -""" -COMPASS equilibrium data reading routines -""" - - -class COMPASSEquilibrium: - - """ - Loads COMPASS efit reconstruction. Either shot number or path has to be passed. - If shot number is specified, path to the efit file is obtained via cdb and data - is loaded. Otherwise path to the efit file has to be specified. - - :param shot_number: COMPASS shot number - :param path: Path to the efit file - :param revision: - """ - def __init__(self, shot_number=None, path=None, revision=1): - self.shot_number = shot_number - - if shot_number is None and path is None: - raise ValueError("Equilibrium has to be specified by either shot number or a path to the efit file") - - if shot_number is not None: - import pyCDB.client - cdb = pyCDB.client.CDBClient() - sig_ref = cdb.get_signal_references(record_number=shot_number, - generic_signal_id=2860, - revision=revision)[0] - path = cdb.get_data_file_reference(**sig_ref) - - self.file_path = path - - self._read_efitfile(self.file_path) - - def _read_efitfile(self, path): - with h5py.File(path, "r") as f: - # read time coordinates - self._time_slice = f["time"].value - - # read 2D poloidal flux map - self._r = f["output/profiles2D/r"].value - self._z = f["output/profiles2D/z"].value - self._psi_grid = f["output/profiles2D/poloidalFlux"].value - - # read poloidal flux on axis and lcfs - self._psi_axis = f["output/globalParameters/psiAxis"].value - self._psi_lcfs = f["output/globalParameters/psiBoundary"].value - - # read positions of magnetic axis - self._magnetic_axis_r = f["output/globalParameters/magneticAxisR"].value - self._magnetic_axis_z = f["output/globalParameters/magneticAxisZ"].value - - # get vacuum toroidal field on axis - self._b_vacuum_radius = 0.56 - self._b_vacuum_magnitude = f["output/globalParameters/bvacRgeom"].value - - # read last close flux surface contour points - self._lcfs_r = f["output/separatrixGeometry/boundaryCoordsR"].value - self._lcfs_z = f["output/separatrixGeometry/boundaryCoordsZ"].value - - # read limiter contour points - self._limiter_r = f["input/limiter/rValues"].value - self._limiter_z = f["input/limiter/zValues"].value - - #get equilibrium important points - self._strikepoint_r = f["output/separatrixGeometry/strikepointCoordsR"].value - self._strikepoint_z = f["output/separatrixGeometry/strikepointCoordsZ"].value - - self._xpoint_r = f["output/separatrixGeometry/xpointCoordsR"].value - self._xpoint_z = f["output/separatrixGeometry/xpointCoordsZ"].value - - - self._rBphi = f["output/fluxFunctionProfiles/rBphi"].value - - self._psin = f["output/fluxFunctionProfiles/normalizedPoloidalFlux"].value - - - def time(self, time): - """ - Returns an equilibrium object for the time-slice closest to the requested time. - The specific time-slice returned is held in the time attribute of the returned object. - :param time: The equilibrium time point. - :returns: An EFITEquilibrium object. - """ - time_index, time_slice = self._find_nearest(self._time_slice, time) - - r = self._r[time_index, :] # get the R coords for psi_grid - z = self._z[time_index, :] # get the Z coords for psi_grid - - psi_grid = self._psi_grid[time_index, :, :] - psi_axis = self._psi_axis[time_index] - psi_lcfs = self._psi_lcfs[time_index] - - magnetic_axis = Point2D(self._magnetic_axis_r[time_index], - self._magnetic_axis_z[time_index]) - - b_vacuum_radius = self._b_vacuum_radius - b_vacuum_magnitude = self._b_vacuum_magnitude[time_index] - - lcfs_polygon = self._process_efit_polygon(self._lcfs_r[time_index, :], - self._lcfs_z[time_index, :]) - - limiter_polygon = self._process_efit_polygon(self._limiter_r[time_index, :], - self._limiter_z[time_index, :]) - - strikepoints = self._process_efit_points(self._strikepoint_r[time_index, :], - self._strikepoint_z[time_index, :]) - - xpoints = self._process_efit_points(self._xpoint_r[time_index, :], - self._xpoint_z[time_index, :]) - - f_profile = np.array([self._psin, self._rBphi[time_index, :]]) - - equilibrium = EFITEquilibrium(r=r, z=z, psi_grid=psi_grid, psi_axis=psi_axis, - psi_lcfs=psi_lcfs, magnetic_axis=magnetic_axis, - x_points = xpoints, strike_points=strikepoints, - f_profile=f_profile, b_vacuum_radius=b_vacuum_radius, - b_vacuum_magnitude=b_vacuum_magnitude, - lcfs_polygon=lcfs_polygon, limiter_polygon=limiter_polygon, - time=time_slice) - - return equilibrium - - def _process_efit_points(self, point_r, point_z): - #avoid using the point if any coordinate is nan - point_list = [] - for i in range(point_r.shape[0]): - if not np.isnan(point_r[i]) and not np.isnan(point_z[i]): - point_list.append(Point2D(point_r[i], point_z[i])) - - return point_list - - def _find_nearest(self, array, value): - if array.min() > value or array.max() < value: - raise IndexError("Requested value is outside the range of the data.") - - idx = np.argmin(np.abs(array - value)) - return idx, array[idx] - - def _process_efit_polygon(self, poly_r, poly_z): - - if poly_r.shape != poly_z.shape: - raise ValueError("EFIT polygon coordinate arrays are inconsistent in length.") - - n = poly_r.shape[0] - if n < 2: - raise ValueError("EFIT polygon coordinate contain less than 2 points.") - - # boundary polygon contains redundant points that must be removed - unique = (poly_r != poly_r[0]) | (poly_z != poly_z[0]) - unique[0] = True # first point must be included! - poly_r = poly_r[unique] - poly_z = poly_z[unique] - - # generate single array containing coordinates - poly_coords = np.zeros((2, len(poly_r))) - poly_coords[0, :] = poly_r - poly_coords[1, :] = poly_z - - return poly_coords - -if __name__ == "__main__": - from os.path import expanduser - from cherab.tools.equilibrium import plot_equilibrium - path = expanduser("~/EFIT/17636.1.h5") - - equilibrium = COMPASSEquilibrium(path=path) - eq_slice = equilibrium.time(1.135) - plot_equilibrium(equilibrium=eq_slice, detail=True, resolution=0.01) \ No newline at end of file diff --git a/cherab/compass/equilibrium/equilibrium_base.py b/cherab/compass/equilibrium/equilibrium_base.py new file mode 100644 index 0000000..f0df165 --- /dev/null +++ b/cherab/compass/equilibrium/equilibrium_base.py @@ -0,0 +1,88 @@ +from raysect.core import Point2D + +from cherab.tools.equilibrium import EFITEquilibrium +from cherab.tools.equilibrium.example import example_equilibrium +from cherab.tools.equilibrium.plot import plot_equilibrium +import cdb_extras.xarray_support as cdbxr +import xarray as xr +import numpy as np + + +class COMPASSEquilibrium_base: + + @property + def data(self): + return self._data + + def _process_efit_polygon(self, poly_r, poly_z): + """ + Processes polygon data to match EFITEquilibrium requirements + :param poly_r: Array with r polygon coordinated. + :param poly_z: Array with z polygon coordinated. + :return: polygon + """ + + if poly_r.shape != poly_z.shape: + raise ValueError("EFIT polygon coordinate arrays are inconsistent in length.") + + n = poly_r.shape[0] + if n < 2: + raise ValueError("EFIT polygon coordinate contain less than 2 points.") + + # boundary polygon contains redundant points that must be removed + unique = (poly_r != poly_r[0]) | (poly_z != poly_z[0]) + unique[0] = True # first point must be included! + poly_r = poly_r[unique] + poly_z = poly_z[unique] + + # generate single array containing coordinates + poly_coords = np.zeros((2, len(poly_r))) + poly_coords[0, :] = poly_r + poly_coords[1, :] = poly_z + + return poly_coords + + def time_sliece(self, time): + """ + Constructs the EFITEquilibrium object for the closes available time to the passed time + :param time: Time to generate EFITEquilibrium object. The nearest available time slice is chosen + :return: EFITEquilibrium object + """ + + time_slice = self._data.sel(time = time, method="nearest") + + r = time_slice.R.values + z = time_slice.Z.values + psi_grid = time_slice.psi_grid.values + psi_axis = time_slice.psi_axis.values + psi_lcfs = time_slice.psi_lcfs.values + magnetic_axis = Point2D(time_slice.R_magnetic_axis.values, time_slice.Z_magnetic_axis.values) + + x_points = [] + for i, v in enumerate(zip(time_slice.R_xpoint.values, time_slice.Z_xpoint.values)): + if not np.isnan(v[0]) and not np.isnan(v[1]): + x_points.append(Point2D(v[0], v[1])) + + + strike_points = [] + for i, v in enumerate(zip(time_slice.R_strike_points.values, time_slice.Z_strike_points.values)): + if not np.isnan(v[0]) and not np.isnan(v[1]): + strike_points.append(Point2D(v[0], v[1])) + + psi_n = time_slice.psi_n.values + + f_profile = np.concatenate((psi_n[np.newaxis, :], time_slice.f_profile.values[np.newaxis, :])) + q_profile = np.concatenate((psi_n[np.newaxis, :], time_slice.q_profile.values[np.newaxis, :])) + + b_vacuum_radius = time_slice.b_vacuum_radius + b_vacuum_magnitude = time_slice.b_vacuum_magnitude.values + + lcfs_polygon = self._process_efit_polygon(time_slice.R_lcfs.values, time_slice.Z_lcfs.values) + limiter_polygon = self._process_efit_polygon(time_slice.R_limiter.values, time_slice.Z_limiter.values) + + + return EFITEquilibrium(r=r, z=z, psi_grid=psi_grid, psi_axis=psi_axis,psi_lcfs=psi_lcfs, magnetic_axis=magnetic_axis, + x_points=x_points, strike_points=strike_points, f_profile=f_profile, q_profile=q_profile, + b_vacuum_radius=b_vacuum_radius, b_vacuum_magnitude=b_vacuum_magnitude, + lcfs_polygon=lcfs_polygon, limiter_polygon=limiter_polygon, + time=time) From 202f8832b8ae004dc6a355d0924f08b4a5d977c9 Mon Sep 17 00:00:00 2001 From: MatejTomes Date: Wed, 30 Oct 2019 20:38:24 +0100 Subject: [PATCH 4/6] Added equilibrium class for compass Upgrade, minor changes in compass eq both classes can be now initialized only from CDB and CUDB databases. --- cherab/compass/equilibrium/compass.py | 56 ++++++-- cherab/compass/equilibrium/compass_upgrade.py | 125 ++++++++++++++++-- .../compass/equilibrium/equilibrium_base.py | 8 +- 3 files changed, 166 insertions(+), 23 deletions(-) diff --git a/cherab/compass/equilibrium/compass.py b/cherab/compass/equilibrium/compass.py index 0fcc495..c91c00d 100644 --- a/cherab/compass/equilibrium/compass.py +++ b/cherab/compass/equilibrium/compass.py @@ -1,9 +1,14 @@ from cherab.compass.equilibrium.equilibrium_base import COMPASSEquilibrium_base +import numpy as np +import xarray as xr class COMPASSEquilibrium(COMPASSEquilibrium_base): + def __init__(self): + super().__init__() + @classmethod - def from_cdb(cls, shot_number): + def from_database(cls, shot_number): """ Obtains EFIT equilibrium data from CDB database for the specified shot. :param shot_number: COMPASS shot number @@ -22,7 +27,7 @@ def from_cdb(cls, shot_number): equilibrium._data = equilibrium._data.to_dataset() equilibrium._data = equilibrium._data.merge(shot["psi_mag_axis/EFIT"].to_dataset()) - equilibrium._data = equilibrium._data.rename({"psi_mag_axis":"psi_axis"}) + equilibrium._data = equilibrium._data.rename({"psi_mag_axis": "psi_axis"}) equilibrium._data = equilibrium._data.merge(shot["psi_lcfs/EFIT"].to_dataset()) @@ -30,9 +35,12 @@ def from_cdb(cls, shot_number): equilibrium._data = equilibrium._data.rename({"R_mag_axis": "R_magnetic_axis", "Z_mag_axis": "Z_magnetic_axis"}) equilibrium._data = equilibrium._data.merge(shot["RZ_xpoint/EFIT"]) - + equilibrium._data = equilibrium._data.rename({"RZ_xpoint_axis1": "xpoint"}) equilibrium._data = equilibrium._data.merge(shot["RZ_strike_points/EFIT"]) + equilibrium._data = equilibrium._data.rename({"RZ_strike_points_axis1": "strike_point", + "R_strike_points": "R_strike_point", + "Z_strike_points": "Z_strike_point"}) equilibrium._data = equilibrium._data.merge(shot["RBphi/EFIT"].to_dataset()) equilibrium._data = equilibrium._data.rename({"RBphi": "f_profile"}) @@ -40,15 +48,43 @@ def from_cdb(cls, shot_number): equilibrium._data = equilibrium._data.merge(shot["q/EFIT"].to_dataset()) equilibrium._data = equilibrium._data.rename({"q": "q_profile"}) - equilibrium._data = equilibrium._data.merge(shot["B_vac_R_geom/EFIT"].to_dataset()) - equilibrium._data = equilibrium._data.rename({"B_vac_R_geom": "b_vacuum_magnitude"}) - equilibrium._data.attrs["b_vacuum_radius"] = 0.56 + equilibrium._data = equilibrium._data.merge(shot["B_vac_R_mag/EFIT"].to_dataset()) + equilibrium._data = equilibrium._data.rename({"B_vac_R_mag": "Btor_vacuum_magnitude"}) + + equilibrium._data = equilibrium._data.merge(shot["R_mag_axis/EFIT"].to_dataset()) + equilibrium._data = equilibrium._data.rename({"R_mag_axis": "Btor_vacuum_radius"}) equilibrium._data = equilibrium._data.merge(shot["RZ_bound/EFIT"]) equilibrium._data = equilibrium._data.rename({"R_bound": "R_lcfs", "Z_bound": "Z_lcfs" }) - equilibrium._data = equilibrium._data.merge(shot["R_limiter_input/EFIT"].to_dataset()) - equilibrium._data = equilibrium._data.merge(shot["Z_limiter_input/EFIT"].to_dataset()) - equilibrium._data = equilibrium._data.rename({"R_limiter_input": "R_limiter", "Z_limiter_input": "Z_limiter" }) + r_limiter = shot["R_limiter_input/EFIT"] + r_limiter.name = "R_limiter" + r_limiter = r_limiter.rename({"R_limiter_input_axis1":"limiter_vertex"}) + + + z_limiter = shot["Z_limiter_input/EFIT"] + z_limiter.name = "Z_limiter" + z_limiter = z_limiter.rename({"Z_limiter_input_axis1": "limiter_vertex"}) + + limiter = xr.merge((r_limiter, z_limiter)) + equilibrium._data = equilibrium._data.merge(limiter) + + equilibrium._data = equilibrium._data.rename({"RZ_bound_axis1": "lcfs_vertex", + "R_axis1": "psi_grid_R_vertex", + "Z_axis2":"psi_grid_Z_vertex"}) + + return equilibrium + + + +if __name__ == "__main__": + eq_shot = COMPASSEquilibrium.from_database(15001) + time = 1100 + + eq = eq_shot.time_slice(time) + + import matplotlib.pyplot as plt + time_slice = eq_shot.data.sel(time = 1100, method="nearest") - return equilibrium \ No newline at end of file + sadf, ax = plt.subplots() + ax.plot(time_slice.R_lcfs.values, time_slice.Z_lcfs.values) \ No newline at end of file diff --git a/cherab/compass/equilibrium/compass_upgrade.py b/cherab/compass/equilibrium/compass_upgrade.py index 341047f..1a3ee7f 100644 --- a/cherab/compass/equilibrium/compass_upgrade.py +++ b/cherab/compass/equilibrium/compass_upgrade.py @@ -1,13 +1,120 @@ from cherab.compass.equilibrium.equilibrium_base import COMPASSEquilibrium_base -import cdb_extras.xarray_support as cdbxr -from pyCDB import client -cdb = client.CDBClient(host="cudb.tok.ipp.cas.cz",data_root="/compass/ -CC19_COMPASS-U_data/") -#class COMPASSUpgradeEquilibrium(COMPASSEquilibrium_base): +import xarray as xr +import numpy as np + + +class COMPASSUpgradeEquilibrium(COMPASSEquilibrium_base): + + def __init__(self): + super().__init__() + + @classmethod + def from_database(cls, shot_number): + """ + Obtains EFIT equilibrium data from CDB database for the specified shot. + :param shot_number: COMPASS shot number + :return: COMPASSEquilibrium class + """ + #set of imports of packages local to COMPASS-site done inside this function to maintain functionality + #off-site + import cdb_extras.xarray_support as cdbxr + from pyCDB import client + + cudb = client.CDBClient(host="cudb.tok.ipp.cas.cz", data_root="/compass/CC19_COMPASS-U_data/") + cdbxr.access_wrappers._global_CDBClient = cudb + + shot = cdbxr.Shot(shot_number) + equilibrium = cls() + + equilibrium._data = shot["psi/fiesta_out"].copy() + equilibrium._data.name = "psi_grid" + equilibrium._data = equilibrium._data.to_dataset() + equilibrium._data = equilibrium._data.rename({"fiesta_r": "R", + "fiesta_z": "Z"}) + + equilibrium._data = equilibrium._data.merge(shot["psi_0/fiesta_out"].to_dataset()) + equilibrium._data = equilibrium._data.rename({"psi_0": "psi_axis"}) + + equilibrium._data = equilibrium._data.merge(shot["psi_boundary/fiesta_out"].to_dataset()) + equilibrium._data = equilibrium._data.rename({"psi_boundary": "psi_lcfs"}) + + equilibrium._data = equilibrium._data.merge(shot["r_mag/fiesta_out"].to_dataset()) + equilibrium._data = equilibrium._data.rename({"r_mag": "R_magnetic_axis"}) + + equilibrium._data = equilibrium._data.merge(shot["z_mag/fiesta_out"].to_dataset()) + equilibrium._data = equilibrium._data.rename({"z_mag": "Z_magnetic_axis"}) + + #add xpoint data + tmp = shot["xp_lower_r/fiesta_out"].expand_dims({"xpoint": 1}) + tmp = xr.concat((tmp, shot["xp_upper_r/fiesta_out"].expand_dims({"xpoint": 1})), dim="xpoint") + tmp.name = "R_xpoint" + tmp = tmp.to_dataset() + equilibrium._data = xr.merge((equilibrium._data, tmp)) + + tmp = shot["xp_lower_z/fiesta_out"].expand_dims({"xpoint": 1}) + tmp = xr.concat((tmp, shot["xp_upper_z/fiesta_out"].expand_dims({"xpoint": 1})), dim="xpoint") + tmp.name = "Z_xpoint" + tmp = tmp.to_dataset() + equilibrium._data = xr.merge((equilibrium._data, tmp)) + + tmp = shot["sp_hfs_r/fiesta_out"].expand_dims({"strike_point": 1}) + tmp = xr.concat((tmp, shot["sp_hfs_r/fiesta_out"].expand_dims({"strike_point": 1})), dim="strike_point") + tmp.name = "R_strike_point" + tmp = tmp.to_dataset() + equilibrium._data = xr.merge((equilibrium._data, tmp)) + + tmp = shot["sp_hfs_z/fiesta_out"].expand_dims({"strike_point": 1}) + tmp = xr.concat((tmp, shot["sp_hfs_z/fiesta_out"].expand_dims({"strike_point": 1})), dim="strike_point") + tmp.name = "Z_strike_point" + tmp = tmp.to_dataset() + + + equilibrium._data = xr.merge((equilibrium._data, tmp)) + + equilibrium._data = equilibrium._data.merge(shot["f/fiesta_out"].to_dataset()) + equilibrium._data = equilibrium._data.rename({"f": "f_profile"}) + + equilibrium._data = equilibrium._data.merge(shot["q/fiesta_out"].to_dataset()) + equilibrium._data = equilibrium._data.rename({"q": "q_profile"}) + + equilibrium._data = equilibrium._data.merge(shot["Bt_vac_mag_axis/fiesta_out"].to_dataset()) + equilibrium._data = equilibrium._data.rename({"Bt_vac_mag_axis": "Btor_vacuum_magnitude"}) + + equilibrium._data = equilibrium._data.merge(shot["Bt_vac_mag_axis/fiesta_out"].to_dataset()) + equilibrium._data = equilibrium._data.rename({"Bt_vac_mag_axis": "Btor_vacuum_radius"}) + + + r_boundary_signal = cudb.get_signal("boundary_r/fiesta_out:{0}".format(shot_number)) + z_boundary_signal = cudb.get_signal("boundary_z/fiesta_out:{0}".format(shot_number)) + + boundary = xr.Dataset({"R_lcfs": (("time", "lcfs_vertex"), r_boundary_signal.data.T), + "Z_lcfs": (("time", "lcfs_vertex"), z_boundary_signal.data.T)}, + coords={"time": r_boundary_signal.time_axis.data}) + + equilibrium._data = xr.merge((equilibrium._data, boundary)) + + tmp =shot["Z_limiter/fiesta_out"] + tmp = xr.Dataset({"Z_limiter": ("limiter_vertex", tmp.Z_limiter.values)}) + equilibrium._data = equilibrium._data.merge(tmp) + + tmp =shot["R_limiter/fiesta_out"] + tmp = xr.Dataset({"R_limiter": ("limiter_vertex", tmp.R_limiter.values)}) + equilibrium._data = equilibrium._data.merge(tmp) + + equilibrium._data = equilibrium._data.rename({"fiesta_psi_norm": "psi_n"}) + + return equilibrium + +if __name__ == "__main__": + from cherab.tools.equilibrium.plot import plot_equilibrium + eq_shot = COMPASSUpgradeEquilibrium.from_database(7400) + time = 2.100 + + time_slice = eq_shot.time_slice(time) + plot_equilibrium(time_slice, detail=True) + + + -# def from_cdb(self, shot_number): -shot_number = 6400 -shot = cudbxr.Shot(shot_number) -shot["Bphi/Fiesta_OUT"] \ No newline at end of file diff --git a/cherab/compass/equilibrium/equilibrium_base.py b/cherab/compass/equilibrium/equilibrium_base.py index f0df165..4c5bc70 100644 --- a/cherab/compass/equilibrium/equilibrium_base.py +++ b/cherab/compass/equilibrium/equilibrium_base.py @@ -42,7 +42,7 @@ def _process_efit_polygon(self, poly_r, poly_z): return poly_coords - def time_sliece(self, time): + def time_slice(self, time): """ Constructs the EFITEquilibrium object for the closes available time to the passed time :param time: Time to generate EFITEquilibrium object. The nearest available time slice is chosen @@ -65,7 +65,7 @@ def time_sliece(self, time): strike_points = [] - for i, v in enumerate(zip(time_slice.R_strike_points.values, time_slice.Z_strike_points.values)): + for i, v in enumerate(zip(time_slice.R_strike_point.values, time_slice.Z_strike_point.values)): if not np.isnan(v[0]) and not np.isnan(v[1]): strike_points.append(Point2D(v[0], v[1])) @@ -74,8 +74,8 @@ def time_sliece(self, time): f_profile = np.concatenate((psi_n[np.newaxis, :], time_slice.f_profile.values[np.newaxis, :])) q_profile = np.concatenate((psi_n[np.newaxis, :], time_slice.q_profile.values[np.newaxis, :])) - b_vacuum_radius = time_slice.b_vacuum_radius - b_vacuum_magnitude = time_slice.b_vacuum_magnitude.values + b_vacuum_radius = time_slice.Btor_vacuum_radius + b_vacuum_magnitude = time_slice.Btor_vacuum_magnitude.values lcfs_polygon = self._process_efit_polygon(time_slice.R_lcfs.values, time_slice.Z_lcfs.values) limiter_polygon = self._process_efit_polygon(time_slice.R_limiter.values, time_slice.Z_limiter.values) From 8decc5ea71eed8b4239ec09eb58b2e968c32c4fd Mon Sep 17 00:00:00 2001 From: MatejTomes Date: Wed, 30 Oct 2019 20:50:36 +0100 Subject: [PATCH 5/6] adding equilibrium imports into __init__.py --- cherab/compass/equilibrium/__init__.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/cherab/compass/equilibrium/__init__.py b/cherab/compass/equilibrium/__init__.py index e0d5f0d..47d8996 100644 --- a/cherab/compass/equilibrium/__init__.py +++ b/cherab/compass/equilibrium/__init__.py @@ -1 +1,2 @@ -from .compass import COMPASSEquilibrium \ No newline at end of file +from .compass import COMPASSEquilibrium +from .compass_upgrade import COMPASSUpgradeEquilibrium \ No newline at end of file From ffaf80f602ffe6b4c3261735414f6452a14d8332 Mon Sep 17 00:00:00 2001 From: MatejTomes Date: Fri, 1 Nov 2019 19:19:39 +0100 Subject: [PATCH 6/6] add metis cudb data source --- cherab/compass/metis/__init__.py | 0 cherab/compass/metis/cudb.py | 61 ++++++++++++++++++++++++++++++++ 2 files changed, 61 insertions(+) create mode 100644 cherab/compass/metis/__init__.py create mode 100644 cherab/compass/metis/cudb.py diff --git a/cherab/compass/metis/__init__.py b/cherab/compass/metis/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/cherab/compass/metis/cudb.py b/cherab/compass/metis/cudb.py new file mode 100644 index 0000000..56815dc --- /dev/null +++ b/cherab/compass/metis/cudb.py @@ -0,0 +1,61 @@ +from cherab.metis.data_source import MetisDataSource_base +import cdb_extras.xarray_support as cdbxr +from pyCDB import client +import numpy as np + +class METISFromCUDB(MetisDataSource_base): + + def __init__(self, shot_number:int = None): + super().__init__() + + self._cudb = client.CDBClient(host="cudb.tok.ipp.cas.cz", data_root="/compass/CC19_COMPASS-U_data/") + cdbxr.access_wrappers._global_CDBClient = self._cudb + self._metis_data_source_id = 2 + + if shot_number is not None: + self.shot_number = shot_number + + + @property + def shot_number(self): + return self._shot_number + + @shot_number.setter + def shot_number(self, value:int): + self._shot_number = value + self._shot = cdbxr.Shot(self._shot_number) + + self.get_data() + + def _get_data(self): + metis_signals = self._cudb.query( + "SELECT gs.generic_signal_name FROM generic_signals gs INNER JOIN data_signals ds USING(generic_signal_id) WHERE gs.data_source_id = {1} AND ds.record_number = {0}".format( + int(self._shot_number), int(self._metis_data_source_id))) + metis_output = self._cudb.query( + "SELECT name FROM data_sources WHERE data_source_id = {0}".format(self._metis_data_source_id)) + + zerod = {} + profil0d = {} + + for i in metis_signals: + if "profil0d" in i["generic_signal_name"]: + key = i["generic_signal_name"].replace("profil0d_", "") + signal_name = i["generic_signal_name"] + "/" + metis_output[0]["name"] + try: + profil0d[key] = self._shot[signal_name].values.T + except Exception: + pass + + elif "zerod" in i["generic_signal_name"]: + key = i["generic_signal_name"].replace("zerod_", "") + signal_name = i["generic_signal_name"] + "/" + metis_output[0]["name"] + try: + zerod[key] = self._shot[signal_name].values + except Exception: + pass + + profil0d["psin"] = np.divide((profil0d["psi"] - profil0d["psi"][0, :]), + profil0d["psi"][-1, :] - profil0d["psi"][0, :]) + + self._zerod_data = zerod + self._profile0d_data = profil0d