diff --git a/cherab/compass/equilibrium/compass_upgrade.py b/cherab/compass/equilibrium/compass_upgrade.py index 1a3ee7f..e11faa0 100644 --- a/cherab/compass/equilibrium/compass_upgrade.py +++ b/cherab/compass/equilibrium/compass_upgrade.py @@ -68,7 +68,6 @@ def from_database(cls, shot_number): 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()) @@ -83,15 +82,11 @@ def from_database(cls, shot_number): 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["boundary_closed_r/fiesta_out"].to_dataset().rename({"boundary_closed_r": "R_lcfs", + "boundary_closed_r_axis1": "lcfs_vertex"}) + tmp = tmp.merge(shot["boundary_closed_z/fiesta_out"].to_dataset().rename({"boundary_closed_z": "Z_lcfs", + "boundary_closed_z_axis1": "lcfs_vertex"})) + equilibrium._data = equilibrium._data.merge(tmp) tmp =shot["Z_limiter/fiesta_out"] tmp = xr.Dataset({"Z_limiter": ("limiter_vertex", tmp.Z_limiter.values)}) @@ -105,16 +100,6 @@ def from_database(cls, shot_number): 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) - - - diff --git a/cherab/compass/math/__init__.py b/cherab/compass/math/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/cherab/compass/math/profile/__init__.py b/cherab/compass/math/profile/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/cherab/compass/math/profile/electron.pxd b/cherab/compass/math/profile/electron.pxd new file mode 100644 index 0000000..b1cee70 --- /dev/null +++ b/cherab/compass/math/profile/electron.pxd @@ -0,0 +1,15 @@ +from raysect.core.math.function cimport Function1D + + +cdef class ModifiedTanhGaussian(Function1D): + + cdef: + double _ped_height, _ped_sol, _ped_width, _ped_pos, _ped_slope, _core_height, _core_width, _core_exp + + cdef double evaluate(self, double x) except? -1e999 + + cdef double pedestal(self, double x) + + cdef double mtanh(self, double x) + + cdef double core(self, double x) \ No newline at end of file diff --git a/cherab/compass/math/profile/electron.pyx b/cherab/compass/math/profile/electron.pyx new file mode 100644 index 0000000..b07b4f4 --- /dev/null +++ b/cherab/compass/math/profile/electron.pyx @@ -0,0 +1,157 @@ +from raysect.core.math.function cimport Function1D + +from libc.math cimport exp + +cdef class ModifiedTanhGaussian(Function1D): + """ + Profile function used at the COMPASS tokamak for electron temperature and density profiles. Pedestal part + is described by a modified hyperbolic tanges. The core part is a Normal-like shape with power as a parameter. + Based on [1]_. + + + :param ped_height: Pedestal height + :param ped_sol: Scrapeoff layer position + :param ped_width: Width of the pedestal + :param ped_pos: Position of the pedestal + :param ped_slope: Slope of the pedestal + :param core_height: Value of the profile at the magnetic axis + :param core_width: Width or shape of the curvature of the core profile part + :param core_exp: Exponential parameter of the bell-shaped core profile + + + [1]Stefanikova, E., Peterka, M., Bohm, P., Bilkova, P., Aftanas, M., Sos, M., Urban, J., Hron, M. and Panek, R., + 2016. Fitting of the Thomson scattering density and temperature profiles on the COMPASS tokamak. + Review of Scientific Instruments, 87(11), p.11E536. + """ + + def __init__(self, ped_height, ped_sol, ped_width, ped_pos, ped_slope, core_height, core_width, core_exp): + super().__init__() + + self.ped_height = ped_height + self.ped_sol = ped_sol + self.ped_width = ped_width + self.ped_pos = ped_pos + self.ped_slope = ped_slope + self.core_height = core_height + self.core_width = core_width + self.core_exp = core_exp + + @property + def ped_height(self): + return self._ped_height + + @ped_height.setter + def ped_height(self, value): + self._ped_height = value + + @property + def ped_sol(self): + return self._ped_sol + + @ped_sol.setter + def ped_sol(self, value): + self._ped_sol = value + + + @property + def ped_width(self): + return self._ped_width + + @ped_width.setter + def ped_width(self, value): + self._ped_width = value + + @property + def ped_pos(self): + return self._ped_pos + + @ped_pos.setter + def ped_pos(self, value): + self._ped_pos = value + + @property + def ped_slope(self): + return self._ped_slope + + @ped_slope.setter + def ped_slope(self, value): + self._ped_slope = value + + @property + def core_height(self): + return self._core_height + + @core_height.setter + def core_height(self, value): + self._core_height = value + + + @property + def core_width(self): + return self._core_width + + @core_width.setter + def core_width(self, value): + self._core_width = value + + @property + def core_exp(self): + return self._core_exp + + @core_exp.setter + def core_exp(self, value): + self._core_exp = value + + cdef double evaluate(self, double x) except? -1e999: + """ + Evaluate the function for x + :param x: Independent variable + :return: Profile value + """ + cdef: + double core, ped + + core = self.core(x) + ped = self.pedestal(x) + + return ped + (self._core_height - ped) * core + + + cdef double pedestal(self, double x): + """ + Pedestal part of the profile. + :param x: Independent variable + :return: Value of pedestal part of profile + """ + cdef: + double mtanh + + mtanh = self.mtanh(x) + + return 0.5 * (self._ped_height - self._ped_sol) * (mtanh + 1) + + cdef double mtanh(self, double x): + """ + Modified hyperbolic tangens used in pedestal shape description + :param x: Independent variable + :return: + """ + cdef: + double xped + + xped = (self._ped_pos - x ) / (2 * self._ped_width) # Independent variable for modified hyp. tangens + return ((1 + self._ped_slope * xped) * exp(xped) - exp(-1 * xped)) / (exp(xped) + exp(-1 * xped)) + + + cdef double core(self, double x): + """ + Core part of the profile described by a Normal-like distribution. + :param x: Independent variable + :return: Core part of the profile + """ + cdef: + double exponent + + exponent = -(x / self._core_width) ** self._core_exp + return exp(exponent) + diff --git a/cherab/compass/metis/cudb.py b/cherab/compass/metis/cudb.py index 56815dc..6e3fbe3 100644 --- a/cherab/compass/metis/cudb.py +++ b/cherab/compass/metis/cudb.py @@ -1,11 +1,14 @@ from cherab.metis.data_source import MetisDataSource_base -import cdb_extras.xarray_support as cdbxr + from pyCDB import client +import cdb_extras.xarray_support as cdbxr + import numpy as np + class METISFromCUDB(MetisDataSource_base): - def __init__(self, shot_number:int = None): + 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/") @@ -15,44 +18,68 @@ def __init__(self, shot_number:int = None): 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): + def shot_number(self, value: int): self._shot_number = value self._shot = cdbxr.Shot(self._shot_number) self.get_data() - def _get_data(self): + @property + def profil0d_signals(self): + return self._profil0d_signals + + @property + def zerod_signals(self): + return self._zerod_signals + + def update_signal_names(self, shot_number=None): + if shot_number is None: + shot_number = self._shot_number + + # get all generic signals belonging to the metis output data source and available for the given shot number 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))) + # get the metis output signal name metis_output = self._cudb.query( "SELECT name FROM data_sources WHERE data_source_id = {0}".format(self._metis_data_source_id)) - zerod = {} - profil0d = {} - + # get dictionariers with profile names as keys and signal names as items + self._profil0d_signals = {} + self._zerod_signals = {} 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 - + self._profil0d_signals[key] = signal_name 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 + self._zerod_signals[key] = signal_name + + def _get_data(self): + # update signal names just in case + self.update_signal_names() + + zerod = {} + profil0d = {} + + for key, item in self._profil0d_signals.items(): + try: + profil0d[key] = self._shot[item].values.T + except Exception: + pass + + for key, item in self._zerod_signals.items(): + try: + zerod[key] = self._shot[item].values + except Exception: + pass profil0d["psin"] = np.divide((profil0d["psi"] - profil0d["psi"][0, :]), profil0d["psi"][-1, :] - profil0d["psi"][0, :]) diff --git a/cherab/compass/utils/nbi/conversions.py b/cherab/compass/utils/nbi/conversions.py new file mode 100644 index 0000000..e69de29