From d3c890e7ab9221c235ec8ef0386818bbf9ea9461 Mon Sep 17 00:00:00 2001 From: marcosertoli <74655352+marcosertoli@users.noreply.github.com> Date: Mon, 18 Sep 2023 17:39:11 +0100 Subject: [PATCH 1/8] Adding capability of recomputing rho-poloidal if fbnd or faxs are corrupted. (#287) Co-authored-by: Marco Sertoli --- indica/equilibrium.py | 24 ++++++++++++++++++------ 1 file changed, 18 insertions(+), 6 deletions(-) diff --git a/indica/equilibrium.py b/indica/equilibrium.py index 239505b2..04020285 100644 --- a/indica/equilibrium.py +++ b/indica/equilibrium.py @@ -65,12 +65,29 @@ def __init__( self.faxs = equilibrium_data["faxs"] self.fbnd = equilibrium_data["fbnd"] self.ftor = equilibrium_data["ftor"] + self.rmag = equilibrium_data["rmag"] + self.rbnd = equilibrium_data["rbnd"] + self.zmag = equilibrium_data["zmag"] + self.zbnd = equilibrium_data["zbnd"] + self.zx = self.zbnd.min("arbitrary_index") self.rhotor = np.sqrt( (self.ftor - self.ftor.sel(rho_poloidal=0.0)) / (self.ftor.sel(rho_poloidal=1.0) - self.ftor.sel(rho_poloidal=0.0)) ) self.psi = equilibrium_data["psi"] - self.rho = np.sqrt((self.psi - self.faxs) / (self.fbnd - self.faxs)) + + # Including workaround in case faxs or fbnd had messy data + rho: DataArray = np.sqrt((self.psi - self.faxs) / (self.fbnd - self.faxs)) + if np.any(np.isnan(rho.interp(R=self.rmag, z=self.zmag))): + self.faxs = self.psi.interp(R=self.rmag, z=self.zmag).drop(["R", "z"]) + self.fbnd = self.psi.interp(R=self.rbnd, z=self.zbnd).mean( + "arbitrary_index" + ) + rho = np.sqrt((self.psi - self.faxs) / (self.fbnd - self.faxs)) + if np.any(np.isnan(rho)): + rho = xr.where(rho > 0, rho, 0.0) + + self.rho = rho self.t = self.rho.t if "vjac" in equilibrium_data and "ajac" in equilibrium_data: self.psin = equilibrium_data["psin"] @@ -85,11 +102,6 @@ def __init__( if "rmji" and "rmjo" in equilibrium_data: self.rmji = equilibrium_data["rmji"] self.rmjo = equilibrium_data["rmjo"] - self.rmag = equilibrium_data["rmag"] - self.rbnd = equilibrium_data["rbnd"] - self.zmag = equilibrium_data["zmag"] - self.zbnd = equilibrium_data["zbnd"] - self.zx = self.zbnd.min("arbitrary_index") self.R_offset = R_shift self.z_offset = z_shift From dd13acdddf172d952777d0fb3613108f392b4fed Mon Sep 17 00:00:00 2001 From: marcosertoli <74655352+marcosertoli@users.noreply.github.com> Date: Tue, 19 Sep 2023 13:33:14 +0100 Subject: [PATCH 2/8] Added BLOM_XY1 reader. TODO: channel mapping still unclear and untouched. SXRC channel mapping also internally reorganised and to be sorted in MDS+. (#286) Co-authored-by: Marco Sertoli --- indica/readers/read_st40.py | 2 ++ indica/readers/st40reader.py | 4 ++++ 2 files changed, 6 insertions(+) diff --git a/indica/readers/read_st40.py b/indica/readers/read_st40.py index 62cdaf62..0927b61e 100644 --- a/indica/readers/read_st40.py +++ b/indica/readers/read_st40.py @@ -24,6 +24,7 @@ "sxr_camera_4": 0, "sxrc_xy1": 0, "sxrc_xy2": 0, + "blom_xy1": 0, "sxr_diode_1": 0, "xrcs": 0, "ts": 0, @@ -42,6 +43,7 @@ "sxr_camera_4": {"brightness": (0, np.inf)}, "sxrc_xy1": {"brightness": (0, np.inf)}, "sxrc_xy2": {"brightness": (0, np.inf)}, + "blom_xy1": {"brightness": (0, np.inf)}, "ts": {"te": (0, 10.0e3), "ne": (0, 1.0e21)}, "pi": {"spectra": (0, np.inf)}, "tws_c": {"spectra": (0, np.inf)}, diff --git a/indica/readers/st40reader.py b/indica/readers/st40reader.py index af54c9ee..de6d3ba4 100644 --- a/indica/readers/st40reader.py +++ b/indica/readers/st40reader.py @@ -87,6 +87,7 @@ class ST40Reader(DataReader): "sxr_camera_4": "get_radiation", "sxrc_xy1": "get_radiation", "sxrc_xy2": "get_radiation", + "blom_xy1": "get_radiation", "cxff_pi": "get_charge_exchange", "cxff_tws_c": "get_charge_exchange", "cxqf_tws_c": "get_charge_exchange", @@ -179,6 +180,9 @@ class ST40Reader(DataReader): "sxrc_xy2": { "brightness": ".profiles:emission", }, + "blom_xy1": { + "brightness": ".profiles:emission", + }, "sxr_diode_1": { "brightness": ".filter_001:signal", }, From 677dcd189985fe025220ca9d8b8e1be38e32d176 Mon Sep 17 00:00:00 2001 From: marcosertoli <74655352+marcosertoli@users.noreply.github.com> Date: Tue, 19 Sep 2023 13:56:58 +0100 Subject: [PATCH 3/8] Marcosertoli/read spd (#275) * Fixed run_tomo_1d minor inconsistencies. Workflows should be refactored to make the code more modular * Refactored readers for new LINES and SXR_SPD diagnostics * Refactored abstractreader test for new diode_filters method * Separated PI and SXRC examples * Deleted commented-out line, generalised test call. * Corrected bug for plotting TS modelled results vs. rho * Simplified general call or read_st40.py wrapper & added example read and plot SXR_SPD * Empty commit --------- Co-authored-by: Marco Sertoli --- indica/models/thomson_scattering.py | 6 +- indica/readers/abstractreader.py | 3 +- indica/readers/read_st40.py | 76 +++---- indica/readers/st40reader.py | 90 ++++----- indica/workflows/run_tomo_1d.py | 186 ++++++++++-------- .../readers/test_abstract_reader_pytest.py | 5 +- 6 files changed, 194 insertions(+), 172 deletions(-) diff --git a/indica/models/thomson_scattering.py b/indica/models/thomson_scattering.py index d74498ee..5396fe3e 100644 --- a/indica/models/thomson_scattering.py +++ b/indica/models/thomson_scattering.py @@ -163,7 +163,8 @@ def example_run( alpha=0.7, ) Ne = bckc["ne"].sel(t=t, method="nearest") - plt.scatter(Ne.rho_poloidal, Ne, color=cols_time[i], marker="o", alpha=0.7) + rho = Ne.transform.rho.sel(t=t, method="nearest") + plt.scatter(rho, Ne, color=cols_time[i], marker="o", alpha=0.7) plt.xlabel("Channel") plt.ylabel("Measured electron density (m^-3)") plt.legend() @@ -176,7 +177,8 @@ def example_run( alpha=0.7, ) Te = bckc["te"].sel(t=t, method="nearest") - plt.scatter(Te.rho_poloidal, Te, color=cols_time[i], marker="o", alpha=0.7) + rho = Te.transform.rho.sel(t=t, method="nearest") + plt.scatter(rho, Te, color=cols_time[i], marker="o", alpha=0.7) plt.xlabel("Channel") plt.ylabel("Measured electron temperature (eV)") plt.legend() diff --git a/indica/readers/abstractreader.py b/indica/readers/abstractreader.py index e34f253e..1bfeec52 100644 --- a/indica/readers/abstractreader.py +++ b/indica/readers/abstractreader.py @@ -862,6 +862,7 @@ def get_diode_filters( t = database_results["times"] t = DataArray(t, coords=[("t", t)], attrs={"long_name": "t", "units": "s"}) + label = database_results["labels"] coords = [("t", t)] if database_results["length"] > 1: coords.append(("channel", np.arange(database_results["length"]))) @@ -876,7 +877,7 @@ def get_diode_filters( coords, transform=transform, ) - data[quantity] = quant_data + data[quantity] = quant_data.assign_coords(label=("channel", label)) return data diff --git a/indica/readers/read_st40.py b/indica/readers/read_st40.py index 0927b61e..34502099 100644 --- a/indica/readers/read_st40.py +++ b/indica/readers/read_st40.py @@ -12,25 +12,25 @@ from indica.readers import ST40Reader from indica.utilities import print_like -REVISIONS = { - "efit": 0, - "brems": -1, - "halpha": -1, - "nirh1": 0, - "smmh1": 0, - "cxff_pi": 0, - "cxff_tws_c": 0, - "cxqf_tws_c": 0, - "sxr_camera_4": 0, - "sxrc_xy1": 0, - "sxrc_xy2": 0, - "blom_xy1": 0, - "sxr_diode_1": 0, - "xrcs": 0, - "ts": 0, - "pi": 0, - "tws_c": 0, -} +INSTRUMENTS: list = [ + "efit", + "lines", + "nirh1", + "smmh1", + "smmh", + "cxff_pi", + "cxff_tws_c", + "cxqf_tws_c", + "sxr_spd", + "sxr_camera_4", + "sxrc_xy1", + "sxrc_xy2", + "sxr_diode_1", + "xrcs", + "ts", + "pi", + "tws_c", +] FILTER_LIMITS = { "cxff_pi": {"ti": (0, np.inf), "vtor": (0, np.inf)}, @@ -39,6 +39,7 @@ "xrcs": {"ti_w": (0, np.inf), "te_kw": (0, np.inf), "te_n3w": (0, np.inf)}, "brems": {"brightness": (0, np.inf)}, "halpha": {"brightness": (0, np.inf)}, + "sxr_spd": {"brightness": (0, np.inf)}, "sxr_diode_1": {"brightness": (0, np.inf)}, "sxr_camera_4": {"brightness": (0, np.inf)}, "sxrc_xy1": {"brightness": (0, np.inf)}, @@ -177,6 +178,7 @@ def map_diagnostics(self, instruments: list, map_raw: bool = False): break def filter_data(self, instruments: list): + if not hasattr(self, "binned_data"): raise ValueError( "Bin data before filtering. No action permitted on raw data structure!" @@ -298,8 +300,8 @@ def plot_profile( def __call__( self, - instruments: list = None, - revisions: dict = None, + instruments: list = [], + revisions: list = None, map_raw: bool = False, tstart: float = None, tend: float = None, @@ -311,10 +313,11 @@ def __call__( debug: bool = False, ): self.debug = debug - if instruments is None: - instruments = list(REVISIONS.keys()) + + if len(instruments) == 0: + instruments = INSTRUMENTS if revisions is None: - revisions = REVISIONS + revisions = [0] * len(instruments) if tstart is None: tstart = self.tstart if tend is None: @@ -324,13 +327,13 @@ def __call__( self.reset_data() self.get_equilibrium(R_shift=R_shift) - for instrument in instruments: + for i, instrument in enumerate(instruments): print(f"Reading {instrument}") if debug: - self.get_raw_data("", instrument, revisions[instrument]) + self.get_raw_data("", instrument, revisions[i]) else: try: - self.get_raw_data("", instrument, revisions[instrument]) + self.get_raw_data("", instrument, revisions[i]) except Exception as e: print(f"Error reading {instrument}: {e}") @@ -363,12 +366,17 @@ def astra_equilibrium(pulse: int, revision: RevisionLike): """Assign ASTRA to equilibrium class""" -def read_cxff_pi(): +def sxr_spd(): import indica.readers.read_st40 as read_st40 + from indica.utilities import set_axis_sci + + st40 = read_st40.ReadST40(11215) + st40(["sxr_spd"]) + + st40.raw_data["sxr_spd"]["brightness"].transform.plot() + + plt.figure() + st40.raw_data["sxr_spd"]["brightness"].sel(channel=0).plot() + set_axis_sci() - st40 = read_st40.ReadST40(10607) - st40(["cxff_pi"]) - st40.raw_data["cxff_pi"]["ti"].los_transform.set_equilibrium( - st40.raw_data["cxff_pi"]["ti"].transform.equilibrium - ) - st40.raw_data["cxff_pi"]["ti"].los_transform.plot() + return st40.raw_data["sxr_spd"]["brightness"] diff --git a/indica/readers/st40reader.py b/indica/readers/st40reader.py index de6d3ba4..1776a0eb 100644 --- a/indica/readers/st40reader.py +++ b/indica/readers/st40reader.py @@ -74,12 +74,13 @@ class ST40Reader(DataReader): "efit": "get_equilibrium", "xrcs": "get_helike_spectroscopy", "princeton": "get_charge_exchange", - "brems": "get_diode_filters", - "halpha": "get_diode_filters", + "lines": "get_diode_filters", "nirh1": "get_interferometry", "nirh1_bin": "get_interferometry", "smmh1": "get_interferometry", + "smmh": "get_interferometry", "astra": "get_astra", + "sxr_spd": "get_diode_filters", "sxr_diode_1": "get_diode_filters", "sxr_diode_2": "get_diode_filters", "sxr_diode_3": "get_diode_filters", @@ -98,12 +99,9 @@ class ST40Reader(DataReader): UIDS_MDS = { "xrcs": "sxr", "princeton": "spectrom", - "brems": "spectrom", - "halpha": "spectrom", "nirh1": "interferom", "nirh1_bin": "interferom", "smmh1": "interferom", - "astra": "", "sxr_diode_1": "sxr", "sxr_diode_2": "sxr", "sxr_diode_3": "sxr", @@ -158,11 +156,14 @@ class ST40Reader(DataReader): "smmh1": { "ne": ".line_int:ne", }, - "brems": { - "brightness": ".brem_mp1:intensity", + "smmh": { + "ne": ".global:ne_int", }, - "halpha": { - "brightness": ".h_alpha_mp1:intensity", + "lines": { + "brightness": ":emission", + }, + "sxr_spd": { + "brightness": ":emission", }, "sxr_camera_4": { "brightness": ".middle_head.filter_4:", @@ -913,47 +914,23 @@ def _get_diode_filters( revision: RevisionLike, quantities: Set[str], ) -> Dict[str, Any]: - + """ + TODO: labels are np.bytes_ type...is this correct? + """ if len(uid) == 0 and instrument in self.UIDS_MDS: uid = self.UIDS_MDS[instrument] - # TODO: change once new MDS+ standardisation has been completed - try: - location, location_path = self._get_signal( - uid, instrument, ".geometry:location", revision - ) - direction, position_path = self._get_signal( - uid, instrument, ".geometry:direction", revision - ) - except TreeNNF: - location = np.array( - [ - [1.0, 0, 0], - ] - ) - direction = ( - np.array( - [ - [0.17, 0, 0], - ] - ) - - location - ) + location, location_path = self._get_signal( + uid, instrument, ".geometry:location", revision + ) + direction, position_path = self._get_signal( + uid, instrument, ".geometry:direction", revision + ) if len(np.shape(location)) == 1: location = np.array([location]) direction = np.array([direction]) length = location[:, 0].size - if instrument == "brems": - _instrument = "lines" - revision = -1 - elif instrument == "halpha": - _instrument = "lines" - revision = -1 - elif "sxr_diode" in instrument: - _instrument = "diode_detr" - else: - raise ValueError(f"{instrument} not yet supported") results: Dict[str, Any] = { "length": length, @@ -961,23 +938,27 @@ def _get_diode_filters( } results["location"] = location results["direction"] = direction - results["revision"] = self._get_revision(uid, _instrument, revision) + results["revision"] = self._get_revision(uid, instrument, revision) results["revision"] = revision revision = results["revision"] quantity = "brightness" qval, q_path = self._get_signal( - uid, _instrument, self.QUANTITIES_MDS[instrument][quantity], revision + uid, instrument, self.QUANTITIES_MDS[instrument][quantity], revision ) - times, _ = self._get_signal_dims(q_path, len(qval.shape)) - times = times[0] + times, _ = self._get_signal(uid, instrument, ":time", revision) + _labels, _ = self._get_signal(uid, instrument, ":label", revision) + if type(_labels[0]) == np.bytes_: + labels = np.array([label.decode("UTF-8") for label in _labels]) + results["times"] = times + results["labels"] = labels results[quantity + "_records"] = q_path results[quantity] = qval try: qval_err, q_path_err = self._get_signal( uid, - _instrument, + instrument, self.QUANTITIES_MDS[instrument][quantity] + "_ERR", revision, ) @@ -996,6 +977,10 @@ def _get_interferometry( revision: RevisionLike, quantities: Set[str], ) -> Dict[str, Any]: + """ + TODO: SMMH 2023 launcher/receiver cross plasma on different poloidal paths! + Currently setting location and direction as average of the two!!!! + """ if len(uid) == 0 and instrument in self.UIDS_MDS: uid = self.UIDS_MDS[instrument] @@ -1014,6 +999,17 @@ def _get_interferometry( direction, direction_path = self._get_signal( uid, instrument, ".geometry:direction", revision ) + + if instrument == "smmh": + location_r, _ = self._get_signal( + uid, instrument, ".geometry:location_r", revision + ) + direction_r, _ = self._get_signal( + uid, instrument, ".geometry:direction_r", revision + ) + location = (location + location_r) / 2.0 + direction = (direction + direction_r) / 2.0 + if len(np.shape(location)) == 1: location = np.array([location]) direction = np.array([direction]) diff --git a/indica/workflows/run_tomo_1d.py b/indica/workflows/run_tomo_1d.py index c4a73fa3..cc5f02b9 100644 --- a/indica/workflows/run_tomo_1d.py +++ b/indica/workflows/run_tomo_1d.py @@ -95,13 +95,15 @@ def phantom_examples( return inverted_emissivity, data_tomo, bckc_tomo -def experimental_examples( +def sxrc_example( + pulse: int = None, instrument: str = "sxrc_xy2", reg_level_guess: float = 0.5, phantom_data: bool = True, plot: bool = True, ): - pulse = PULSES[instrument] + if pulse is None: + pulse = PULSES[instrument] tstart = 0.02 tend = 0.1 @@ -112,17 +114,101 @@ def experimental_examples( quantity = list(st40.binned_data[instrument])[0] los_transform = st40.binned_data[instrument][quantity].transform los_transform.set_equilibrium(equilibrium, force=True) - if instrument == "pi": - BREMS_MODEL.set_los_transform(los_transform) - attrs = st40.binned_data[instrument]["spectra"].attrs - background, brightness = BREMS_MODEL.integrate_spectra( - st40.binned_data[instrument]["spectra"] + model = SXR_MODEL + model.set_los_transform(los_transform) + + if phantom_data: + plasma = example_plasma( + pulse, + tstart=tstart, + tend=tend, + dt=dt, ) - background.attrs = attrs - brightness.attrs = attrs - st40.binned_data[instrument]["background"] = background - st40.binned_data[instrument]["brightness"] = brightness - model = SXR_MODEL + plasma.build_atomic_data() + plasma.set_equilibrium(equilibrium) + model.set_plasma(plasma) + bckc = model() + emissivity = model.emissivity + brightness = bckc["brightness"] + else: + emissivity = None + brightness = st40.binned_data[instrument]["brightness"] + + z = los_transform.z + R = los_transform.R + dl = los_transform.dl + + data_t0 = brightness.isel(t=0).data + has_data = np.logical_not(np.isnan(brightness.isel(t=0).data)) & (data_t0 > 0) + rho_equil = equilibrium.rho.interp(t=brightness.t) + input_dict = dict( + brightness=brightness.data, + dl=dl, + t=brightness.t.data, + R=R, + z=z, + rho_equil=dict( + R=rho_equil.R.data, + z=rho_equil.z.data, + t=rho_equil.t.data, + rho=rho_equil.data, + ), + has_data=has_data, + debug=False, + ) + if emissivity is not None: + input_dict["emissivity"] = emissivity + + tomo = tomo_1D.SXR_tomography(input_dict, reg_level_guess=reg_level_guess) + tomo() + + if plot: + model.los_transform.plot() + tomo.show_reconstruction() + + inverted_emissivity = DataArray( + tomo.emiss, coords=[("t", tomo.tvec), ("rho_poloidal", tomo.rho_grid_centers)] + ) + inverted_error = DataArray( + tomo.emiss_err, + coords=[("t", tomo.tvec), ("rho_poloidal", tomo.rho_grid_centers)], + ) + inverted_emissivity.attrs["error"] = inverted_error + data_tomo = brightness + bckc_tomo = DataArray(tomo.backprojection, coords=data_tomo.coords) + + return inverted_emissivity, data_tomo, bckc_tomo + + +def pi_bremsstrahlung_example( + pulse: int = None, + instrument: str = "pi", + reg_level_guess: float = 0.5, + phantom_data: bool = True, + plot: bool = True, +): + if pulse is None: + pulse = PULSES[instrument] + + tstart = 0.02 + tend = 0.1 + dt = 0.01 + st40 = ReadST40(pulse, tstart, tend, dt=dt) + st40(instruments=[instrument, "efit"]) + equilibrium = Equilibrium(st40.raw_data["efit"]) + quantity = list(st40.binned_data[instrument])[0] + los_transform = st40.binned_data[instrument][quantity].transform + los_transform.set_equilibrium(equilibrium, force=True) + model = BREMS_MODEL + model.set_los_transform(los_transform) + attrs = st40.binned_data[instrument]["spectra"].attrs + background, brightness = BREMS_MODEL.integrate_spectra( + st40.binned_data[instrument]["spectra"] + ) + background.attrs = attrs + brightness.attrs = attrs + st40.binned_data[instrument]["background"] = background + st40.binned_data[instrument]["brightness"] = brightness if phantom_data: plasma = example_plasma( @@ -198,7 +284,7 @@ def sxrc_xy( reg_level_guess: float = 0.5, input_dict: dict = None, channels=slice(0, 15), - save_fig: bool = True, + save_fig: bool = False, instrument="sxrc_xy2", ): if input_dict is None: @@ -212,7 +298,7 @@ def sxrc_xy( equil = st40.equilibrium rho, theta = data.transform.convert_to_rho_theta(t=data.t) - rho = rho.sel(channel=channels) + # rho = rho.sel(channel=channels) brightness = data @@ -341,75 +427,3 @@ def old_camera( tomo.show_reconstruction() return input_dict - - -def fake_data( - pulse: int = 9229, - plasma=None, - model=None, - bckc=None, - debug=False, - exclude_bad_points=True, - plot=True, - reg_level_guess: float = 0.5, - nchannels=12, - input_dict: dict = None, -): - from indica.models.sxr_camera import example_run as sxr_example - - if input_dict is None: - if plasma is None or model is None or bckc is None: - plasma, model, bckc = sxr_example(pulse=pulse, nchannels=nchannels) - - tstart = plasma.t.min().values - tend = plasma.t.max().values - st40 = ReadST40(pulse, tstart, tend) - st40(instruments=["sxrc_xy2", "sxr_camera_4"]) - - try: - data = st40.binned_data["sxrc_xy2"]["brightness"] - except KeyError: - data = st40.binned_data["sxr_camera_4"]["brightness"] - model.set_los_transform(data.transform) - model.los_transform.set_equilibrium(plasma.equilibrium, force=True) - model() - - z = model.los_transform.z - R = model.los_transform.R - dl = model.los_transform.dl - - brightness = model.bckc["brightness"] - data_t0 = brightness.isel(t=0).data - if exclude_bad_points: - has_data = np.logical_not(np.isnan(data_t0)) & (data_t0 >= 1.0e3) - else: - has_data = np.logical_not(np.isnan(data_t0)) - - rho_equil = plasma.equilibrium.rho.interp(t=brightness.t) - input_dict = dict( - emissivity=model.emissivity, - brightness=brightness.data, - dl=dl, - t=brightness.t.data, - R=R, - z=z, - rho_equil=dict( - R=rho_equil.R.data, - z=rho_equil.z.data, - t=rho_equil.t.data, - rho=rho_equil.data, - ), - debug=debug, - has_data=has_data, - ) - - # return input_dict - - tomo = tomo_1D.SXR_tomography(input_dict, reg_level_guess=reg_level_guess) - - tomo() - - if plot: - tomo.show_reconstruction() - - return plasma diff --git a/tests/unit/readers/test_abstract_reader_pytest.py b/tests/unit/readers/test_abstract_reader_pytest.py index 894e2c92..cccd6aba 100644 --- a/tests/unit/readers/test_abstract_reader_pytest.py +++ b/tests/unit/readers/test_abstract_reader_pytest.py @@ -36,7 +36,7 @@ class Reader(DataReader): "radiation": "get_radiation", "helike_spectroscopy": "get_helike_spectroscopy", "interferometry": "get_interferometry", - "filters": "get_diode_filters", + "diode_filters": "get_diode_filters", } def __init__( @@ -368,6 +368,7 @@ def _get_diode_filters( results["location"] = np.array([[1.0, 2.0, 3.0]] * results["length"]) results["direction"] = np.array([[1.0, 2.0, 3.0]] * results["length"]) + results["labels"] = np.array(["label"] * results["length"]) for quantity in quantities: results[quantity] = np.random.uniform(0, 1.0e6, (nt, results["length"])) @@ -494,7 +495,7 @@ def test_get_helike_spectroscopy(): def test_get_diode_filters(): _test_get_methods( - instrument="filters", + instrument="diode_filters", nsamples=10, ) From 0ef4332a47041f1d7eebe24950a31747a94eee0d Mon Sep 17 00:00:00 2001 From: marcosertoli <74655352+marcosertoli@users.noreply.github.com> Date: Thu, 21 Sep 2023 13:09:29 +0100 Subject: [PATCH 4/8] Marcosertoli/xrcs wavelenght hack (#288) * Adding capability of recomputing rho-poloidal if fbnd or faxs are corrupted. * Adding capability of recomputing rho-poloidal if fbnd or faxs are corrupted. * Merging latest reader mods from st40 * Fixed linting issue --------- Co-authored-by: Marco Sertoli --- indica/models/diode_filters.py | 2 +- indica/readers/read_st40.py | 4 ++-- indica/readers/st40reader.py | 11 +++++++---- 3 files changed, 10 insertions(+), 7 deletions(-) diff --git a/indica/models/diode_filters.py b/indica/models/diode_filters.py index 860a610a..f42f5253 100644 --- a/indica/models/diode_filters.py +++ b/indica/models/diode_filters.py @@ -31,7 +31,7 @@ def __init__( filter_fwhm: float = 1, # 1 filter_type: str = "boxcar", etendue: float = 1.0, - calibration: float = 2.0e-5, + calibration: float = 1, instrument_method="get_diode_filters", channel_mask: slice = None, # =slice(18, 28), ): diff --git a/indica/readers/read_st40.py b/indica/readers/read_st40.py index 34502099..8b6539b2 100644 --- a/indica/readers/read_st40.py +++ b/indica/readers/read_st40.py @@ -301,7 +301,7 @@ def plot_profile( def __call__( self, instruments: list = [], - revisions: list = None, + revisions: dict = None, map_raw: bool = False, tstart: float = None, tend: float = None, @@ -317,7 +317,7 @@ def __call__( if len(instruments) == 0: instruments = INSTRUMENTS if revisions is None: - revisions = [0] * len(instruments) + revisions = {instrument: 0 for instrument in instruments} if tstart is None: tstart = self.tstart if tend is None: diff --git a/indica/readers/st40reader.py b/indica/readers/st40reader.py index 1776a0eb..dc92d5ad 100644 --- a/indica/readers/st40reader.py +++ b/indica/readers/st40reader.py @@ -702,11 +702,14 @@ def _get_helike_spectroscopy( direction = np.array([direction]) results["times"], _ = self._get_signal(uid, instrument, ":time", revision) - results["wavelength"], _ = self._get_signal( - uid, instrument, ":wavelength", revision - ) + wavelength, _ = self._get_signal(uid, instrument, ":wavelength", revision) # TODO: change once wavelength in MDS+ has been fixed to nanometers! - results["wavelength"] /= 10.0 + wavelength /= 10.0 + if self.pulse >= 10307: + dlambda = float(np.abs(wavelength[1] - wavelength[0])) * 4 + wavelength += dlambda + results["wavelength"] = wavelength + for q in quantities: qval, q_path = self._get_signal( uid, instrument, self.QUANTITIES_MDS[instrument][q], revision From fb3529414ab30a51eb070f573dc3d8e56ba50751 Mon Sep 17 00:00:00 2001 From: marcosertoli <74655352+marcosertoli@users.noreply.github.com> Date: Tue, 3 Oct 2023 14:04:32 +0100 Subject: [PATCH 5/8] Fixed bug to read desired revision in ReadST40 wrapper that was affecting the assigned equilibrium. (#290) Co-authored-by: Marco Sertoli --- indica/readers/read_st40.py | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/indica/readers/read_st40.py b/indica/readers/read_st40.py index 8b6539b2..bcd2c577 100644 --- a/indica/readers/read_st40.py +++ b/indica/readers/read_st40.py @@ -318,6 +318,12 @@ def __call__( instruments = INSTRUMENTS if revisions is None: revisions = {instrument: 0 for instrument in instruments} + for instr in instruments: + if instr not in revisions.keys(): + revisions[instr] = 0 + if "efit" not in revisions: + revisions["efit"] = 0 + if tstart is None: tstart = self.tstart if tend is None: @@ -326,14 +332,14 @@ def __call__( dt = self.dt self.reset_data() - self.get_equilibrium(R_shift=R_shift) - for i, instrument in enumerate(instruments): + self.get_equilibrium(R_shift=R_shift, revision=revisions["efit"]) + for instrument in instruments: print(f"Reading {instrument}") if debug: - self.get_raw_data("", instrument, revisions[i]) + self.get_raw_data("", instrument, revisions[instrument]) else: try: - self.get_raw_data("", instrument, revisions[i]) + self.get_raw_data("", instrument, revisions[instrument]) except Exception as e: print(f"Error reading {instrument}: {e}") From 683004d97d03217827d7b79ee9cfbee36a80e77e Mon Sep 17 00:00:00 2001 From: marcosertoli <74655352+marcosertoli@users.noreply.github.com> Date: Mon, 9 Oct 2023 09:18:42 +0100 Subject: [PATCH 6/8] Updated st40 reader to use "smaug" server. (#291) Co-authored-by: Marco Sertoli --- indica/plotters/plot_time_evolution.py | 216 +++++++++++++++---------- indica/readers/st40reader.py | 2 +- indica/writers/bda_tree.py | 2 +- 3 files changed, 134 insertions(+), 86 deletions(-) diff --git a/indica/plotters/plot_time_evolution.py b/indica/plotters/plot_time_evolution.py index bc4e897c..2c41cb2c 100644 --- a/indica/plotters/plot_time_evolution.py +++ b/indica/plotters/plot_time_evolution.py @@ -1,11 +1,12 @@ +from copy import deepcopy import getpass import matplotlib.pyplot as plt import numpy as np import xarray as xr -from indica.numpy_typing import ArrayLike from indica.readers.read_st40 import ReadST40 +from indica.utilities import save_figure from indica.utilities import set_axis_sci from indica.utilities import set_plot_colors from indica.utilities import set_plot_rcparams @@ -14,99 +15,146 @@ CMAP, COLORS = set_plot_colors() LINESTYLES = ["solid", "dashed", "dotted"] set_plot_rcparams("profiles") +QUANTITIES: list = [ + "smmh:ne", + "xrcs:ti_w", + "xrcs:te_n3w", + "xrcs:te_kw", + "xrcs:spectra", + "sxrc_xy2:brightness", + "efit:ipla", + "efit:rmag", + "efit:zmag", + "sxr_spd:brightness", + # "cxff_pi:ti", + # "cxff_pi:vtor", + # "cxff_tws_c:ti", + # "cxff_tws_c:vtor", + # "ts:te", + # "ts:ne", + # "lines", +] +Y0 = {} +Y0["nirh1"] = True +Y0["smmh1"] = True +Y0["xrcs"] = True +Y0["sxr_diode_1"] = True +Y0["brems"] = True +Y0["efit"] = True +Y0["cxff_tws_c"] = True +Y0["cxff_pi"] = True -def plot_sxrc( - pulse: int = 10607, # 10605 - what_to_read: list = ["sxrc_xy2:brightness"], - tplot: ArrayLike = [0.067, 0.069], # [0.07, 0.09] - save_fig=False, - plot_raw: bool = False, - xvar: str = "t", - yvar: str = "channel", - tstart=0.02, - tend=0.1, - dt: float = 0.001, - data_key="binned_data", + +def plot_st40_data( + pulses: list, + tstart: float = 0.02, + tend: float = 0.1, + dt: float = 0.005, + quantities: list = [], + tplot: float = None, + save_fig: bool = False, + fig_path: str = None, + fig_style: str = "profiles", + plot_binned: bool = True, ): - instruments = [] - quantities = [] - linestyles = {} - for i, identifier in enumerate(what_to_read): - instr, quant = identifier.split(":") - instruments.append(instr) - quantities.append(quant) - linestyles[instr] = LINESTYLES[i] - - st40 = ReadST40(pulse, tstart, tend, dt=dt) - st40(instruments=instruments) - - data: dict = {} - for instr, quant in zip(instruments, quantities): - if instr not in data.keys(): - data[instr] = {} - data[instr][quant] = getattr(st40, data_key)[instr][quant].transpose(yvar, xvar) - _rho, _theta = data[instr][quant].transform.convert_to_rho_theta( - t=data[instr][quant].t - ) + if tplot is None: + tplot = np.mean([tstart, tend]) + if len(quantities) == 0: + quantities = QUANTITIES + instruments = list(np.unique([quant.split(":")[0] for quant in quantities])) + if fig_path is None: + fig_path = f"{FIG_PATH}" - for instr, quant in zip(instruments, quantities): - plt.figure() - plot = data[instr][quant].plot(label=f"{instr}:{quant}") - set_axis_sci(plot_object=plot) - plt.title(f"{instr.upper()} for pulse {pulse}") + set_plot_rcparams(fig_style) + xr.set_options(keep_attrs=True) + colors = CMAP(np.linspace(0.75, 0.1, len(pulses), dtype=float)) - if tplot is not None: - cols = CMAP(np.linspace(0.75, 0.1, np.size(tplot), dtype=float)) + raw: dict = {quant: {} for quant in quantities} + binned = deepcopy(raw) + data = {"raw": raw, "binned": binned} + for i, pulse in enumerate(pulses): + fig_path += f"_{pulse}" + st40 = ReadST40(pulse, tstart=tstart, tend=tend, dt=dt) + st40(instruments) + for quantity in quantities: + instr, quant = quantity.split(":") + data["binned"][quantity][pulse] = st40.binned_data[instr][quant] + data["raw"][quantity][pulse] = st40.raw_data[instr][quant] + for quantity in quantities: + print(quantity) + instr, _ = quantity.split(":") plt.figure() - for icol, t in enumerate(tplot): - for instr, quant in zip(instruments, quantities): - _data = data[instr][quant].sel(t=t, method="nearest") - if "error" in _data.attrs: - _err = (_data.error + _data.stdev).sel(t=t, method="nearest") - else: - _err = xr.full_like(_data, 0.0) - _R_diff = ( - _data.transform.impact_parameter.value - - _data.transform.equilibrium.rmag - ).sel(t=t, method="nearest") - plt.fill_between( - _R_diff.values, - _data.values + _err.values, - _data.values - _err.values, - color=cols[icol], - alpha=0.5, - ) - plt.plot( - _R_diff.values, - _data.values, - label=f"t={t:.3f}", - linestyle=linestyles[instr], - color=cols[icol], + for i, pulse in enumerate(pulses): + color = colors[i] + if len(data["raw"][quantity].keys()) == 0: + continue + plot_data(data, quantity, pulse, tplot, key="raw", color=color) + if plot_binned: + plot_data(data, quantity, pulse, tplot, key="binned", color=color) + + set_axis_sci() + if instr in Y0.keys(): + plt.ylim( + 0, ) - plt.xlabel("Impact R - R$_{mag}$") - plt.ylabel(f"{_data.long_name} [{_data.units}]") - plt.title(f"{instr.upper()} for pulse {pulse}") - plt.legend() - set_axis_sci() - - return st40 - - -# def plot_time_surface( -# pulse: int = 10605, -# instruments: list = ["sxrc_xy1"], -# quantity: str = "brightness", -# tplot: ArrayLike = None, -# save_fig=False, -# plot_raw: bool = False, -# xvar: str = "time", -# yvar: str = "channel", -# ): + plt.legend() + plt.autoscale() + save_figure(fig_path, f"{quantity}", save_fig=save_fig) + + +def plot_data(data, quantity: str, pulse: int, tplot: float, key="raw", color=None): + str_to_add = "" + instr, quant = quantity.split(":") + if key == "raw": + marker = None + else: + marker = "o" + + _data = data[key][quantity][pulse] + tslice = slice(_data.t.min().values, _data.t.max().values) + if "error" not in _data.attrs: + _data.attrs["error"] = xr.full_like(_data, 0.0) + if "stdev" not in _data.attrs: + _data.attrs["stdev"] = xr.full_like(_data, 0.0) + _err = np.sqrt(_data.error**2 + _data.stdev**2) + _err = xr.where(_err / _data.values < 1.0, _err, 0.0) + if len(_data.dims) > 1: + str_to_add = f" @ {tplot:.3f} s" + tslice = _data.t.sel(t=tplot, method="nearest") + + _data = _data.sel(t=tslice) + _err = _err.sel(t=tslice) + if instr in "xrcs" and quant == "spectra": + bgnd = _data.sel(wavelength=slice(0.393, 0.388)).mean("wavelength") + _data -= bgnd + label = None + if key == "raw": + label = str(pulse) + alpha = 0.5 + plt.fill_between( + _data.coords[_data.dims[0]].values, + _data.values - _err.values, + _data.values + _err.values, + color=color, + alpha=alpha, + ) + if key == "binned": + alpha = 0.8 + plt.errorbar( + _data.coords[_data.dims[0]].values, + _data.values, + _err.values, + color=color, + alpha=alpha, + ) + _data.plot(label=label, color=color, alpha=alpha, marker=marker) + plt.title(f"{instr.upper()} {quant}" + str_to_add) + if __name__ == "__main__": plt.ioff() - plot_sxrc() + plot_st40_data([11226, 11225]) plt.show() diff --git a/indica/readers/st40reader.py b/indica/readers/st40reader.py index dc92d5ad..3ba45951 100644 --- a/indica/readers/st40reader.py +++ b/indica/readers/st40reader.py @@ -317,7 +317,7 @@ def __init__( pulse: int, tstart: float, tend: float, - server: str = "192.168.1.7", # 192.168.1.7 10.0.40.13 + server: str = "smaug", tree: str = "ST40", default_error: float = 0.05, max_freq: float = 1e6, diff --git a/indica/writers/bda_tree.py b/indica/writers/bda_tree.py index 5a5878a8..6dfa20ae 100644 --- a/indica/writers/bda_tree.py +++ b/indica/writers/bda_tree.py @@ -231,7 +231,7 @@ def check_analysis_run( which_run, ): # Checker function to see if data already exists in a run - IP_address_smaug = "192.168.1.7:8000" + IP_address_smaug = "smaug" conn = Connection(IP_address_smaug) conn.openTree("BDA", pulseNo) From 8ea96911fdba883bebffa2e11daae479f0df3e6d Mon Sep 17 00:00:00 2001 From: marcosertoli <74655352+marcosertoli@users.noreply.github.com> Date: Thu, 12 Oct 2023 13:36:55 +0100 Subject: [PATCH 7/8] Marcosertoli/dev (#292) * First commit new dev branch * Fixed bug to read desired revision in ReadST40 wrapper that was affecting the assigned equilibrium. * Modified TS limits to physical range. Corrected bug in equilibrium reader time settings. * Still having issues with TS MDS+ inconsistencies -> hardcoded x, y, z until this is solved. * Working implementation of simple spline fitting routine + example for TS data. * Minor corrections to load_modelling.., but still too complex. * Added TODO notes for things that need to be changed. Issues should be raised once these CAN be changed. * Ran pre-commit * Minor mods and to decrease complexity and pass linting --------- Co-authored-by: Marco Sertoli --- indica/converters/time.py | 2 +- indica/operators/spline_fit_easy.py | 95 ++++--- indica/readers/available_quantities.py | 1 + indica/readers/read_st40.py | 10 +- indica/readers/st40reader.py | 27 +- indica/workflows/load_modelling_plasma.py | 309 ++++++++++++++-------- 6 files changed, 277 insertions(+), 167 deletions(-) diff --git a/indica/converters/time.py b/indica/converters/time.py index cd4ced5f..7c35bd72 100644 --- a/indica/converters/time.py +++ b/indica/converters/time.py @@ -329,7 +329,7 @@ def bin_in_time_dt(tstart: float, tend: float, dt: float, data: DataArray) -> Da ------- : Array like the input, but binned along the time axis. - + TODO: add possibility of doing 50% overlap of time bins! """ check_bounds_bin(tstart, tend, dt, data) tlabels = get_tlabels_dt(tstart, tend, dt) diff --git a/indica/operators/spline_fit_easy.py b/indica/operators/spline_fit_easy.py index a012f5c0..d6549067 100644 --- a/indica/operators/spline_fit_easy.py +++ b/indica/operators/spline_fit_easy.py @@ -99,63 +99,86 @@ def residuals(yknots): def spline_fit_ts( - pulse: int, - tstart: float = 0.0, - tend: float = 0.2, - plot: bool = False, + pulse: int = 11314, + tstart: float = 0.03, + tend: float = 0.1, + dt: float = 0.01, + quantity: str = "te", + R_shift: float = 0.0, + knots: list = None, + plot: bool = True, ): + st40 = ReadST40(pulse, tstart=tstart, tend=tend, dt=dt) + st40(["ts"], R_shift=R_shift) + + if quantity == "te" and knots is None: + knots = [0, 0.3, 0.6, 0.8, 1.1] + if quantity == "ne" and knots is None: + knots = [0, 0.3, 0.6, 0.8, 0.95, 1.1] + data_all = st40.raw_data["ts"][quantity] + t = data_all.t + transform = data_all.transform + transform.convert_to_rho_theta(t=data_all.t) - ST40 = ReadST40(pulse, tstart=tstart, tend=tend) - ST40(["ts"]) - - Te_data_all = ST40.binned_data["ts"]["te"] - t = ST40.binned_data["ts"]["te"].t - transform = ST40.binned_data["ts"]["te"].transform R = transform.R Rmag = transform.equilibrium.rmag.interp(t=t) # Fit all available TS data - ind = np.full_like(Te_data_all, True) + ind = np.full_like(data_all, True) rho = xr.where(ind, transform.rho, np.nan) - Te_data = xr.where(ind, Te_data_all, np.nan) - Te_err = xr.where(ind, Te_data_all.error, np.nan) - Te_fit = fit_profile( - rho, Te_data, Te_err, knots=[0, 0.3, 0.5, 0.75, 0.95, 1.05], virtual_knots=False - ) + data = xr.where(ind, data_all, np.nan) + err = xr.where(ind, data_all.error, np.nan) + fit = fit_profile(rho, data, err, knots=knots, virtual_knots=False) # Use only HFS channels ind = R <= Rmag - rho = xr.where(ind, transform.rho, np.nan) - Te_data = xr.where(ind, Te_data_all, np.nan) - Te_err = xr.where(ind, Te_data_all.error, np.nan) - Te_fit_hfs = fit_profile(rho, Te_data, Te_err, virtual_knots=True) + rho_hfs = xr.where(ind, transform.rho, np.nan) + data_hfs = xr.where(ind, data_all, np.nan) + err_hfs = xr.where(ind, data_all.error, np.nan) + fit_hfs = fit_profile(rho_hfs, data_hfs, err_hfs, knots=knots, virtual_knots=True) + + # Use only LFS channels + ind = R >= Rmag + rho_lfs = xr.where(ind, transform.rho, np.nan) + data_lfs = xr.where(ind, data_all, np.nan) + err_lfs = xr.where(ind, data_all.error, np.nan) + fit_lfs = fit_profile(rho_lfs, data_lfs, err_lfs, knots=knots, virtual_knots=True) if plot: - for t in Te_data_all.t: + for t in data_all.t: plt.ioff() plt.errorbar( - Te_data_all.transform.rho.sel(t=t), - Te_data_all.sel(t=t), - Te_data_all.error.sel(t=t), + rho_hfs.sel(t=t), + data_hfs.sel(t=t), + err_hfs.sel(t=t), + marker="o", + label="data HFS", + color="blue", + ) + plt.errorbar( + rho_lfs.sel(t=t), + data_lfs.sel(t=t), + err_lfs.sel(t=t), marker="o", - label="data", + label="data LFS", + color="red", ) - Te_fit.sel(t=t).plot( - linewidth=5, alpha=0.5, color="orange", label="spline fit all" + fit.sel(t=t).plot( + linewidth=5, alpha=0.5, color="black", label="spline fit all" ) - Te_fit_hfs.sel(t=t).plot( - linewidth=5, alpha=0.5, color="red", label="spline fit HFS" + fit_lfs.sel(t=t).plot( + linewidth=5, alpha=0.5, color="red", label="spline fit LFS" + ) + fit_hfs.sel(t=t).plot( + linewidth=5, alpha=0.5, color="blue", label="spline fit HFS" ) plt.legend() plt.show() - return Te_data_all, Te_fit + return data_all, fit if __name__ == "__main__": - spline_fit_ts( - 10619, - tstart=0.0, - tend=0.2, - plot=True, - ) + plt.ioff() + spline_fit_ts(11089, quantity="ne") + plt.show() diff --git a/indica/readers/available_quantities.py b/indica/readers/available_quantities.py index bb5d9fb9..133d892e 100644 --- a/indica/readers/available_quantities.py +++ b/indica/readers/available_quantities.py @@ -1,5 +1,6 @@ """ Quantities that can be read with the current abstract reader implementation +TODO: change the tuple to DataArray (long_name, units) - see examples in abstractreader """ from typing import Dict diff --git a/indica/readers/read_st40.py b/indica/readers/read_st40.py index bcd2c577..4589dd00 100644 --- a/indica/readers/read_st40.py +++ b/indica/readers/read_st40.py @@ -45,7 +45,7 @@ "sxrc_xy1": {"brightness": (0, np.inf)}, "sxrc_xy2": {"brightness": (0, np.inf)}, "blom_xy1": {"brightness": (0, np.inf)}, - "ts": {"te": (0, 10.0e3), "ne": (0, 1.0e21)}, + "ts": {"te": (0, np.inf), "ne": (0, np.inf)}, "pi": {"spectra": (0, np.inf)}, "tws_c": {"spectra": (0, np.inf)}, } @@ -91,8 +91,12 @@ def __init__( self.tend = tend self.dt = dt - self.reader = ST40Reader(pulse, tstart - 0.05, tend + 0.05, tree=tree) - self.reader_equil = ST40Reader(pulse, tstart - 0.1, tend + 0.1, tree=tree) + _tend = tend + dt * 2 + _tstart = tstart - dt * 2 + if _tstart < 0: + _tstart = 0.0 + self.reader = ST40Reader(pulse, _tstart, _tend, tree=tree) + self.reader_equil = ST40Reader(pulse, _tstart, _tend, tree=tree) self.equilibrium: Equilibrium self.raw_data: dict = {} diff --git a/indica/readers/st40reader.py b/indica/readers/st40reader.py index 3ba45951..92f567c2 100644 --- a/indica/readers/st40reader.py +++ b/indica/readers/st40reader.py @@ -3,6 +3,7 @@ """ +from copy import deepcopy from typing import Any from typing import Dict from typing import List @@ -96,6 +97,7 @@ class ST40Reader(DataReader): "tws_c": "get_spectrometer", "ts": "get_thomson_scattering", } + # TODO: this will not be necessary once the MDS+ standardisation is complete UIDS_MDS = { "xrcs": "sxr", "princeton": "spectrom", @@ -296,7 +298,8 @@ class ST40Reader(DataReader): }, } - _IMPLEMENTATION_QUANTITIES = { # TODO: these will be different diagnostics!!!!!!! + # TODO: this can be deleted once MDS+ standardisation is complete + _IMPLEMENTATION_QUANTITIES = { "diode_arrays": { # GETTING THE DATA OF THE SXR CAMERA "sxr_camera_1": ("brightness", "total"), "sxr_camera_2": ("brightness", "50_Al_filtered"), @@ -305,6 +308,7 @@ class ST40Reader(DataReader): }, } + # TODO: this can be deleted once MDS+ standardisation is complete _RADIATION_RANGES = { "sxr_camera_1": (1, 20), "sxr_camera_2": (21, 40), @@ -768,6 +772,7 @@ def _get_charge_exchange( z, z_path = self._get_signal(uid, instrument, ":z", revision) R, R_path = self._get_signal(uid, instrument, ":R", revision) + # TODO: temporary fix until geometry sorted (especially pulse if statement..) try: location, location_path = self._get_signal( uid, instrument, ".geometry:location", revision @@ -779,7 +784,6 @@ def _get_charge_exchange( location = np.array([location]) direction = np.array([direction]) - # TODO: temporary fix until geometry sorted if location.shape[0] != x.shape[0]: if self.pulse > 10200: index = np.arange(18, 36) @@ -787,7 +791,6 @@ def _get_charge_exchange( index = np.arange(21, 36) location = location[index] direction = direction[index] - except TreeNNF: location = None direction = None @@ -812,10 +815,8 @@ def _get_charge_exchange( ) except TreeNNF: qval_err = np.full_like(qval, 0.0) - # q_path_err = "" dimensions, _ = self._get_signal_dims(q_path, len(qval.shape)) - results[q + "_records"] = q_path results[q] = qval results[f"{q}_error"] = qval_err @@ -868,13 +869,6 @@ def _get_spectrometer( location = np.array([location]) direction = np.array([direction]) - # if self.pulse > 10200: - # index = np.arange(18, 36) - # else: - # index = np.arange(21, 36) - # location = location[index] - # direction = direction[index] - for q in quantities: qval, q_path = self._get_signal( uid, @@ -892,7 +886,6 @@ def _get_spectrometer( ) except TreeNNF: qval_err = np.full_like(qval, 0.0) - # q_path_err = "" dimensions, _ = self._get_signal_dims(q_path, len(qval.shape)) @@ -953,6 +946,8 @@ def _get_diode_filters( _labels, _ = self._get_signal(uid, instrument, ":label", revision) if type(_labels[0]) == np.bytes_: labels = np.array([label.decode("UTF-8") for label in _labels]) + else: + labels = _labels results["times"] = times results["labels"] = labels @@ -1085,6 +1080,8 @@ def _get_thomson_scattering( revision = results["revision"] times, times_path = self._get_signal(uid, instrument, ":time", revision) + # TODO: hardcoded correction to TS coordinates to be fixed in MDS+ + print("\n Hardcoded correction to TS coordinates to be fixed in MDS+ \n") # location, location_path = self._get_signal( # uid, instrument, ".geometry:location", revision # ) @@ -1095,6 +1092,9 @@ def _get_thomson_scattering( y, y_path = self._get_signal(uid, instrument, ":y", revision) z, z_path = self._get_signal(uid, instrument, ":z", revision) R, R_path = self._get_signal(uid, instrument, ":R", revision) + z = R * 0.0 + x = deepcopy(R) + y = 0 for q in quantities: qval, q_path = self._get_signal( @@ -1170,6 +1170,7 @@ def get_revision_name(self, revision) -> str: """Return string defining RUN## or BEST if revision = 0""" if type(revision) == int: + rev_str = "" if revision < 0: rev_str = "" elif revision == 0: diff --git a/indica/workflows/load_modelling_plasma.py b/indica/workflows/load_modelling_plasma.py index 957577bc..d3167b0b 100644 --- a/indica/workflows/load_modelling_plasma.py +++ b/indica/workflows/load_modelling_plasma.py @@ -11,6 +11,7 @@ from indica.converters.line_of_sight import LineOfSightTransform from indica.converters.transect import TransectCoordinates from indica.equilibrium import Equilibrium +from indica.models.bolometer_camera import Bolometer from indica.models.charge_exchange import ChargeExchange from indica.models.diode_filters import BremsstrahlungDiode from indica.models.helike_spectroscopy import HelikeSpectrometer @@ -29,6 +30,7 @@ CMAP, COLORS = set_plot_colors() DIAGNOSTIC_MODELS = { + "smmh": Interferometry, "smmh1": Interferometry, "nirh1": Interferometry, "xrcs": HelikeSpectrometer, @@ -37,10 +39,23 @@ "cxff_tws_c": ChargeExchange, "cxqf_tws_c": ChargeExchange, "brems": BremsstrahlungDiode, - "sxr_diode_1": SXRcamera, - "sxr_camera_4": SXRcamera, + "sxrc_xy1": Bolometer, + "sxrc_xy2": SXRcamera, + "sxr_spd": SXRcamera, + "blom_xy1": Bolometer, } -INSTRUMENTS = ["smmh1", "nirh1", "xrcs", "sxr_diode_1", "efit", "brems"] +INSTRUMENTS: list = [ + "smmh", + "smmh1", + "nirh1", + "xrcs", + "blom_xy1", + "sxrc_xy1", + "sxrc_xy2", + "efit", + "sxr_spd", +] +REVISIONS: dict = {instr: 0 for instr in INSTRUMENTS} FIG_PATH = f"/home/{getpass.getuser()}/figures/Indica/load_modelling_examples/" plt.ion() @@ -201,9 +216,9 @@ def initialize_diagnostic_models( def example_run( dt: float = 0.01, verbose: bool = True, - example: str = "predictive", + example: str = "alsu_11314", all_runs: bool = False, - plot: bool = False, + plot: bool = True, save_fig: bool = False, fig_style="profiles", alpha: float = 1.0, @@ -216,6 +231,7 @@ def example_run( pulse_code, pulse, equil, + equil_run, code, runs, comment, @@ -237,7 +253,8 @@ def example_run( if pulse is not None: print("Reading ST40 data") st40 = ReadST40(pulse, tstart, tend, dt=dt, tree="st40") - st40(instruments=INSTRUMENTS, map_diagnostics=False) + REVISIONS["efit"] = equil_run + st40(instruments=INSTRUMENTS, map_diagnostics=False, revisions=REVISIONS) if equil != code: equilibrium = {run: Equilibrium(st40.raw_data[equil]) for run in runs} @@ -283,7 +300,10 @@ def example_run( models[instrument].set_plasma(_plasma) - _bckc = models[instrument]() + if instrument == "xrcs": + _bckc = models[instrument](moment_analysis=True) + else: + _bckc = models[instrument]() bckc[run][instrument] = deepcopy(_bckc) if plot or save_fig: @@ -501,11 +521,13 @@ def plot_data_bckc_comparison( run_str += f"-{runs[1]}" fig_path = f"{FIG_PATH}{pulse_code}_{time}_{code}_{run_str}/" - norm = {} - norm["brems"] = True - norm["sxr_camera_4"] = True - norm["sxr_diode_1"] = True - norm["xrcs"] = True + norm: dict = {} + norm["xrcs"] = {} + norm["xrcs"]["spectra"] = True + # norm["brems"] = True + # norm["sxr_camera_4"] = True + # norm["sxrc_xy2"] = {} + # norm["sxrc_xy2"]["brightness"] = True y0 = {} y0["nirh1"] = True y0["smmh1"] = True @@ -518,39 +540,38 @@ def plot_data_bckc_comparison( for run in runs: bckc[run]["efit"] = {"wp": plasma[run].wp} - for instrument in st40.raw_data.keys(): - for quantity in bckc[run][instrument].keys(): - print(instrument) - print(f" {quantity}") - run = runs[0] - + instruments = st40.raw_data.keys() + for instrument in instruments: + quantities = bckc[runs[0]][instrument].keys() + for quantity in quantities: if ( quantity not in st40.binned_data[instrument].keys() or quantity not in st40.raw_data[instrument].keys() ): continue - plt.figure() + print(instrument) + print(f" {quantity}") _raw = st40.raw_data[instrument][quantity] _binned = st40.binned_data[instrument][quantity] - _bckc = bckc[run][instrument][quantity] + _bckc = bckc[runs[0]][instrument][quantity] + tslice = slice(_bckc.t.min().values, _bckc.t.max().values) + str_to_add = "" + tslice_raw = tslice + tslice_binned = tslice + if "error" not in _binned.attrs: _binned.attrs["error"] = xr.full_like(_binned, 0.0) if "stdev" not in _binned.attrs: _binned.attrs["stdev"] = xr.full_like(_binned, 0.0) - err = np.sqrt(_binned.error**2 + _binned.stdev**2) err = xr.where(err / _binned.values < 1.0, err, 0.0) - if len(_bckc.dims) > 1: + if len(_binned.dims) > 1: str_to_add = f" @ {time:.3f} s" tslice_binned = _binned.t.sel(t=time, method="nearest") tslice_raw = _raw.t.sel(t=time, method="nearest") - else: - str_to_add = "" - tslice_raw = tslice - tslice_binned = tslice _raw = _raw.sel(t=tslice_raw) _binned = _binned.sel(t=tslice_binned) @@ -562,6 +583,7 @@ def plot_data_bckc_comparison( _binned -= bgnd _raw -= bgnd + plt.figure() _raw.plot( label="Raw", color=COLORS["raw"], @@ -583,17 +605,14 @@ def plot_data_bckc_comparison( markersize=markersize, ) - label = None + label: str = "Model" for run in runs: - if run == runs[-1]: - label = "Model" - _bckc = bckc[run][instrument][quantity].sel(t=tslice_binned) if instrument in norm.keys(): - _bckc = _bckc * _binned.max() / _bckc.max() - - if label is not None: - label += " (scaled)" + if quantity in norm[instrument].keys(): + _bckc = _bckc / _bckc.max() * _binned.max() + if label is not None: + label += " (scaled)" (_bckc).plot( label=label, @@ -601,91 +620,151 @@ def plot_data_bckc_comparison( linewidth=rcParams["lines.linewidth"] * 1.5, alpha=alpha, ) - set_axis_sci() - plt.title(f"{instrument.upper()} {quantity}" + str_to_add) - if instrument in y0.keys(): - plt.ylim( - 0, - ) + del label - if quantity == "spectra": - # TODO: wavelength axis is sorted from max to min... - plt.xlim(_bckc.wavelength.min(), _bckc.wavelength.max()) + set_axis_sci() + plt.title(f"{instrument.upper()} {quantity}" + str_to_add) + if instrument in y0.keys(): + plt.ylim( + 0, + ) + if quantity == "spectra": + plt.xlim(_bckc.wavelength.min(), _bckc.wavelength.max()) - plt.legend() - save_figure(fig_path, f"{instrument}_{quantity}", save_fig=save_fig) + plt.legend() + save_figure(fig_path, f"{instrument}_{quantity}", save_fig=save_fig) def example_params(example: str, all_runs: bool = False): - comment: str - pulse_code: int - pulse: int - equil: str - code: str - tstart: float - tend: float - tplot: float - runs = [f"RUN{run}" for run in (500 + np.arange(61, 77))] - - if example == "predictive": - comment = "Tests using fixed-boundary predictive ASTRA" - pulse_code = 13110009 - pulse = 10009 - equil = "astra" - code = "astra" - if not all_runs: - runs = ["RUN2621"] - tstart = 0.02 - tend = 0.08 - tplot = 0.06 - elif example == "interpretative_10009": - comment = "interpretative ASTRA using HDA/EFIT" - pulse_code = 13110009 - pulse = 10009 - equil = "efit" - code = "astra" - if not all_runs: - runs = ["RUN564"] - tstart = 0.03 - tend = 0.1 - tplot = 0.06 - elif example == "interpretative_9850": - comment = "ASTRA interpretative using HDA/EFIT" - pulse = 9850 - pulse_code = 13109850 - equil = "efit" - code = "astra" - if not all_runs: - runs = ["RUN564"] - tstart = 0.02 - tend = 0.1 - tplot = 0.08 - elif example == "interpretative_9229": - comment = "ASTRA interpretative using HDA/EFIT" - pulse = 9229 - pulse_code = 13109229 - equil = "efit" - code = "astra" - if not all_runs: - runs = ["RUN572"] - tstart = 0.03 - tend = 0.11 - tplot = 0.06 - elif example == "diverted": - comment = "predictive ASTRA using for diverted scenario" - pulse_code = 13000040 - pulse = 10009 - equil = "astra" - code = "astra" - if not all_runs: - runs = ["RUN292"] - tstart = 0.03 - tend = 0.11 - tplot = 0.1 - else: - raise ValueError(f"No parameters for example {example}") - - return pulse_code, pulse, equil, code, runs, comment, tstart, tend, tplot + runs_all: list = [f"RUN{run}" for run in (500 + np.arange(61, 77))] + + params = { + "predictive": dict( + comment="Tests using fixed-boundary predictive ASTRA", + pulse_code=13110009, + pulse=10009, + equil="astra", + equil_run=0, + code="astra", + runs=["RUN2621"], + tstart=0.02, + tend=0.08, + tplot=0.06, + ), + "interpretative_10009": dict( + comment="interpretative ASTRA using HDA/EFIT", + pulse_code=13110009, + pulse=10009, + equil="efit", + equil_run=0, + code="astra", + runs=["RUN564"], + tstart=0.03, + tend=0.1, + tplot=0.06, + ), + "interpretative_9850": dict( + comment="ASTRA interpretative using HDA/EFIT", + pulse=9850, + pulse_code=13109850, + equil="efit", + equil_run=0, + code="astra", + runs=["RUN564"], + tstart=0.02, + tend=0.1, + tplot=0.08, + ), + "interpretative_9229": dict( + comment="ASTRA interpretative using HDA/EFIT", + pulse=9229, + pulse_code=13109229, + equil="efit", + equil_run=0, + code="astra", + runs=["RUN572"], + tstart=0.03, + tend=0.11, + tplot=0.06, + ), + "diverted": dict( + comment="predictive ASTRA using for diverted scenario", + pulse_code=13000040, + pulse=10009, + equil="astra", + equil_run=0, + code="astra", + runs=["RUN292"], + tstart=0.03, + tend=0.11, + tplot=0.1, + ), + "aleksei_11228": dict( + comment="ASTRA using TS and invented Ti shapes", + pulse=11228, + pulse_code=13011228, + equil="efit", # "astra" + equil_run=1, + code="astra", + runs=["RUN610", "RUN611", "RUN612"], + tstart=0.03, + tend=0.11, + tplot=0.08, + ), + "alsu_11312": dict( + comment="ASTRA using TS and peaked Ti scaled to CXRS", + pulse=11312, + pulse_code=33011312, + equil="astra", + equil_run=0, + code="astra", + runs=["RUN21"], + tstart=0.065, + tend=0.095, + tplot=0.075, + ), + "alsu_11314": dict( + comment="ASTRA using TS and peaked Ti scaled to CXRS", + pulse=11314, + pulse_code=33011314, + equil="astra", + equil_run=0, + code="astra", + runs=["RUN12"], + tstart=0.065, + tend=0.095, + tplot=0.075, + ), + "alsu_11317": dict( + comment="ASTRA using TS and peaked Ti scaled to CXRS", + pulse=11317, + pulse_code=33011317, + equil="astra", + equil_run=0, + code="astra", + runs=["RUN9"], + tstart=0.065, + tend=0.095, + tplot=0.075, + ), + } + + _params = params[example] + if all_runs: + _params["runs"] = runs_all + + return ( + _params["pulse_code"], + _params["pulse"], + _params["equil"], + _params["equil_run"], + _params["code"], + _params["runs"], + _params["comment"], + _params["tstart"], + _params["tend"], + _params["tplot"], + ) # def add_gacode_data( @@ -723,4 +802,6 @@ def example_params(example: str, all_runs: bool = False): if __name__ == "__main__": + plt.ioff() example_run() + plt.show() From b7b7364d8036e7fcaecd16b5187abfd8bb4c624a Mon Sep 17 00:00:00 2001 From: Marco Sertoli Date: Thu, 26 Oct 2023 11:20:17 +0100 Subject: [PATCH 8/8] Minor bug fixes --- indica/plotters/plot_time_evolution.py | 11 ++++++++--- indica/workflows/load_modelling_plasma.py | 16 ++++++++++++++-- 2 files changed, 22 insertions(+), 5 deletions(-) diff --git a/indica/plotters/plot_time_evolution.py b/indica/plotters/plot_time_evolution.py index 2c41cb2c..7bf950fa 100644 --- a/indica/plotters/plot_time_evolution.py +++ b/indica/plotters/plot_time_evolution.py @@ -30,8 +30,8 @@ # "cxff_pi:vtor", # "cxff_tws_c:ti", # "cxff_tws_c:vtor", - # "ts:te", - # "ts:ne", + "ts:te", + "ts:ne", # "lines", ] @@ -69,7 +69,10 @@ def plot_st40_data( set_plot_rcparams(fig_style) xr.set_options(keep_attrs=True) - colors = CMAP(np.linspace(0.75, 0.1, len(pulses), dtype=float)) + if len(pulses) > 1: + colors = CMAP(np.linspace(0.75, 0.1, len(pulses), dtype=float)) + else: + colors = ["blue"] raw: dict = {quant: {} for quant in quantities} binned = deepcopy(raw) @@ -104,6 +107,8 @@ def plot_st40_data( plt.autoscale() save_figure(fig_path, f"{quantity}", save_fig=save_fig) + return data, st40 + def plot_data(data, quantity: str, pulse: int, tplot: float, key="raw", color=None): str_to_add = "" diff --git a/indica/workflows/load_modelling_plasma.py b/indica/workflows/load_modelling_plasma.py index d3167b0b..2b530a42 100644 --- a/indica/workflows/load_modelling_plasma.py +++ b/indica/workflows/load_modelling_plasma.py @@ -617,7 +617,7 @@ def plot_data_bckc_comparison( (_bckc).plot( label=label, color=COLORS["bckc"], - linewidth=rcParams["lines.linewidth"] * 1.5, + linewidth=rcParams["lines.linewidth"], alpha=alpha, ) del label @@ -699,6 +699,18 @@ def example_params(example: str, all_runs: bool = False): tend=0.11, tplot=0.1, ), + "michail_10014": dict( + comment="predictive ASTRA", + pulse_code=36010014, + pulse=10014, + equil="efit", + equil_run=0, + code="astra", + runs=["RUN24"], + tstart=0.03, + tend=0.1, + tplot=0.07, + ), "aleksei_11228": dict( comment="ASTRA using TS and invented Ti shapes", pulse=11228, @@ -803,5 +815,5 @@ def example_params(example: str, all_runs: bool = False): if __name__ == "__main__": plt.ioff() - example_run() + example_run(example="alsu_11314") plt.show()