diff --git a/tofu/__init__.py b/tofu/__init__.py index b2defb789..4ef97c747 100644 --- a/tofu/__init__.py +++ b/tofu/__init__.py @@ -91,7 +91,7 @@ lsubout = ['tofu.{0}'.format(ss) for ss in lsubout] msg = "\nThe following subpackages are not available:" msg += "\n - " + "\n - ".join(lsubout) - msg += "\n => see tofu.dsub[] for details." + msg += "\n => see print(tofu.dsub[]) for details." warnings.warn(msg) # ------------------------------------- diff --git a/tofu/data/_comp.py b/tofu/data/_comp.py index 2bbe96d4e..9af8e8117 100644 --- a/tofu/data/_comp.py +++ b/tofu/data/_comp.py @@ -432,6 +432,8 @@ def func(pts, vect=None, t=None, ntall=ntall, vquant[indtq[indtu[ii]], :], trifinder=trifind )(r, z).filled(fill_value) + if np.any(np.isnan(val)) and not np.isnan(fill_value): + val[np.isnan(val)] = fill_value return val, t # -------------------- @@ -463,6 +465,8 @@ def func(pts, vect=None, t=None, ntall=ntall, ind = indt == indtu[ii] val[ind, indok] = vquant[indtq[indtu[ii]], indpts[indok]] + if np.any(np.isnan(val)) and not np.isnan(fill_value): + val[np.isnan(val)] = fill_value return val, t @@ -518,14 +522,13 @@ def func(pts, vect=None, t=None, ntall=ntall, # interpolate 1d ind = indt == indtu[ii] val[ind, ...] = scpinterp.interp1d( - vr1[indtr1[indtu[ii]], :], + vr1[indtr1[ii], :], vquant[indtq[indtu[ii]], :], kind='linear', bounds_error=False, fill_value=fill_value )(np.asarray(vii)) val[np.isnan(val)] = fill_value - return val, t else: @@ -563,12 +566,16 @@ def func(pts, vect=None, t=None, ntall=ntall, # interpolate 1d ind = indt == indtu[ii] val[ind, indok] = scpinterp.interp1d( - vr1[indtr1[indtu[ii]], :], + vr1[indtr1[ii], :], vquant[indtq[indtu[ii]], :], kind='linear', bounds_error=False, fill_value=fill_value )(vr2[indtr2[ii], indpts[indok]]) + + # Double check nan in case vr2 is itself nan on parts of mesh + if np.any(np.isnan(val)) and not np.isnan(fill_value): + val[np.isnan(val)] = fill_value return val, t return func @@ -672,13 +679,13 @@ def func(pts, vect=None, t=None, ntall=ntall, tall=tall, tbinall=tbinall): # Get pts in (r,z,phi) - r, z = np.hypot(pts[0,:],pts[1,:]), pts[2,:] - phi = np.arctan2(pts[1,:],pts[0,:]) + r, z = np.hypot(pts[0, :], pts[1, :]), pts[2, :] + phi = np.arctan2(pts[1, :], pts[0, :]) # Deduce vect in (r,z,phi) - vR = np.cos(phi)*vect[0,:] + np.sin(phi)*vect[1,:] - vphi = -np.sin(phi)*vect[0,:] + np.cos(phi)*vect[1,:] - vZ = vect[2,:] + vR = np.cos(phi)*vect[0, :] + np.sin(phi)*vect[1, :] + vPhi = -np.sin(phi)*vect[0, :] + np.cos(phi)*vect[1, :] + vZ = vect[2, :] # Prepare output shapeval = list(pts.shape) @@ -689,6 +696,7 @@ def func(pts, vect=None, t=None, ntall=ntall, # Interpolate indpts = trifind(r, z) + indok = indpts > -1 if t is None: for ii in range(0,ntall): valR[ii, ...] = vq2dR[indtq[ii], indpts] @@ -697,9 +705,8 @@ def func(pts, vect=None, t=None, ntall=ntall, t = tall else: ntall, indt, indtu = plasma._get_indtu(t=t, tall=tall, - tbinall=tbinall, - idref1d=idref1d, - idref2d=idref2d)[1:] + tbinall=tbinall)[1:-2] + for ii in range(0, ntall): ind = indt == indtu[ii] valR[ind, ...] = vq2dR[indtq[indtu[ii]], indpts] diff --git a/tofu/data/_core.py b/tofu/data/_core.py index dae5917c3..a4f45c47e 100644 --- a/tofu/data/_core.py +++ b/tofu/data/_core.py @@ -2026,7 +2026,11 @@ def _extract_common_params(obj0, obj1=None): @staticmethod def _recreatefromoperator(d0, other, opfunc): - if type(other) in [int,float,np.int64,np.float64]: + if type(other) in [int, float, np.int64, np.float64]: + data = opfunc(d0.data, other) + dcom = d0._extract_common_params(d0) + + elif isinstance(other, np.ndarray): data = opfunc(d0.data, other) dcom = d0._extract_common_params(d0) @@ -2515,6 +2519,9 @@ def trifind(r, z, indpts = indR*nZ + indZ else: indpts = indZ*nR + indR + indout = ((r < R[0]) | (r > R[-1]) + | (z < Z[0]) | (z > Z[-1])) + indpts[indout] = -1 return indpts dd[dk][k0]['R'] = R @@ -3735,7 +3742,7 @@ def _get_finterp(self, idquant=None, idref1d=None, idref2d=None, idq2dR=None, idq2dPhi=None, idq2dZ=None, interp_t=None, interp_space=None, - fill_value=np.nan, ani=False, Type=None): + fill_value=None, ani=False, Type=None): if interp_t is None: interp_t = 'nearest' @@ -3853,7 +3860,7 @@ def _checkformat_qr12RPZ(self, quant=None, ref1d=None, ref2d=None, def get_finterp2d(self, quant=None, ref1d=None, ref2d=None, q2dR=None, q2dPhi=None, q2dZ=None, interp_t=None, interp_space=None, - fill_value=np.nan, Type=None): + fill_value=None, Type=None): """ Return the function interpolating (X,Y,Z) pts on a 1d/2d profile Can be used as input for tf.geom.CamLOS1D/2D.calc_signal() @@ -3880,7 +3887,7 @@ def interp_pts2profile(self, pts=None, vect=None, t=None, quant=None, ref1d=None, ref2d=None, q2dR=None, q2dPhi=None, q2dZ=None, interp_t=None, interp_space=None, - fill_value=np.nan, Type=None): + fill_value=None, Type=None): """ Return the value of the desired profiles_1d quantity For the desired inputs points (pts): @@ -3911,8 +3918,23 @@ def interp_pts2profile(self, pts=None, vect=None, t=None, else: idmesh = [id_ for id_ in self._ddata[idref2d]['depend'] if self._dindref[id_]['group'] == 'mesh'][0] - pts = self.dmesh[idmesh]['data']['nodes'] - pts = np.array([pts[:,0], np.zeros((pts.shape[0],)), pts[:,1]]) + if self.dmesh[idmesh]['data']['type'] == 'rect': + if self.dmesh[idmesh]['data']['shapeRZ'] == ('R', 'Z'): + R = np.repeat(self.dmesh[idmesh]['data']['R'], + self.dmesh[idmesh]['data']['nZ']) + Z = np.tile(self.dmesh[idmesh]['data']['Z'], + self.dmesh[idmesh]['data']['nR']) + else: + R = np.tile(self.dmesh[idmesh]['data']['R'], + self.dmesh[idmesh]['data']['nZ']) + Z = np.repeat(self.dmesh[idmesh]['data']['Z'], + self.dmesh[idmesh]['data']['nR']) + pts = np.array( + [R, np.zeros((self.dmesh[idmesh]['data']['size'],)), Z]) + else: + pts = self.dmesh[idmesh]['data']['nodes'] + pts = np.array( + [pts[:, 0], np.zeros((pts.shape[0],)), pts[:, 1]]) pts = np.atleast_2d(pts) if pts.shape[0] != 3: @@ -3945,7 +3967,7 @@ def calc_signal_from_Cam(self, cam, t=None, quant=None, ref1d=None, ref2d=None, q2dR=None, q2dPhi=None, q2dZ=None, Brightness=True, interp_t=None, - interp_space=None, fill_value=np.nan, + interp_space=None, fill_value=None, res=0.005, DL=None, resMode='abs', method='sum', ind=None, out=object, plot=True, dataname=None, fs=None, dmargin=None, wintit=None, invert=True, diff --git a/tofu/data/_plot.py b/tofu/data/_plot.py index d99cb0abe..4101ba8d9 100644 --- a/tofu/data/_plot.py +++ b/tofu/data/_plot.py @@ -364,9 +364,9 @@ def _DataCam12D_plot(lData, key=None, nchMax=_nchMax, ntMax=_ntMax, c0 = [all([dd.dlabels[kk] == lData[0].dlabels[kk] for dd in lData[1:]]) for kk in ['t','X','data']] if not all(c0): - msg = "All Data objects must have the same:\n" + msg = "All Data objects do not have the same:\n" msg += " dlabels['t'], dlabels['X'] and dlabels['data'] !" - raise Exception(msg) + warnings.warn(msg) # --------- @@ -387,8 +387,9 @@ def _DataCam12D_plot(lData, key=None, nchMax=_nchMax, ntMax=_ntMax, # Check nch and X c0 = [dd.nch == lData[0].nch for dd in lData[1:]] if not all(c0): - msg = "All Data objects must have the same number of channels (self.nch)" - msg += "\nYou can set the indices of the channels with self.set_indch()" + msg = ("All Data objects must have the same nb. of channels\n" + + "\t- self.nch: {}\n".format([dd.nch for dd in lData]) + + "\n => use self.set_indch()") raise Exception(msg) nch = lData[0].nch diff --git a/tofu/geom/_core.py b/tofu/geom/_core.py index acfa98a4c..736a05477 100644 --- a/tofu/geom/_core.py +++ b/tofu/geom/_core.py @@ -5968,6 +5968,9 @@ def _calc_signal_postformat( if Brightness is False: if dataname is None: dataname = r"LOS-integral x Etendue" + if E is None or np.all(np.isnan(E)): + msg = "Cannot use etendue, it was not set properly !" + raise Exception(msg) if t is None or len(t) == 1 or E.size == 1: sig = sig * E else: @@ -6207,7 +6210,7 @@ def calc_signal_from_Plasma2D( self, plasma2d, t=None, - newcalc=False, + newcalc=True, quant=None, ref1d=None, ref2d=None, @@ -6252,7 +6255,9 @@ def calc_signal_from_Plasma2D( if newcalc: # Get time vector - if t is None: + lc = [t is None, type(t) is str, type(t) is np.ndarray] + assert any(lc) + if lc[0]: out = plasma2d._checkformat_qr12RPZ( quant=quant, ref1d=ref1d, @@ -6262,6 +6267,8 @@ def calc_signal_from_Plasma2D( q2dZ=q2dZ, ) t = plasma2d._get_tcom(*out[:4])[0] + elif lc[1]: + t = plasma2d._ddata[t]['data'] else: t = np.atleast_1d(t).ravel() @@ -6359,7 +6366,7 @@ def funcbis(*args, **kwdargs): nbrep = np.r_[ indpts[0], np.diff(indpts), pts.shape[1] - indpts[-1] ] - vect = np.repeat(self.u, nbrep, axis=1) + vect = -np.repeat(self.u, nbrep, axis=1) if fill_value is None: fill_value = 0. diff --git a/tofu/imas2tofu/_core.py b/tofu/imas2tofu/_core.py index 612a17b67..6e941f810 100644 --- a/tofu/imas2tofu/_core.py +++ b/tofu/imas2tofu/_core.py @@ -154,8 +154,14 @@ class MultiIDSLoader(object): 'dim':'B', 'quant':'BT', 'units':'T'}, '2dBZ':{'str':'time_slice[time].ggd[0].b_field_z[0].values', 'dim':'B', 'quant':'BZ', 'units':'T'}, - '2dmeshNodes':{'str':'grids_ggd[0].grid[0].space[0].objects_per_dimension[0].object[].geometry'}, - '2dmeshFaces':{'str':'grids_ggd[0].grid[0].space[0].objects_per_dimension[2].object[].nodes'}}, + '2dmeshNodes': {'str': ('grids_ggd[0].grid[0].space[0]' + + '.objects_per_dimension[0]' + + '.object[].geometry')}, + '2dmeshFaces': {'str': ('grids_ggd[0].grid[0].space[0]' + + '.objects_per_dimension[2]' + + '.object[].nodes')}, + '2dmeshR': {'str': 'time_slice[0].profiles_2d[0].r'}, + '2dmeshZ': {'str': 'time_slice[0].profiles_2d[0].z'}}, 'core_profiles': {'t':{'str':'time'}, @@ -251,15 +257,16 @@ class MultiIDSLoader(object): 'floop_flux':{'str':'flux_loop[chan].flux.data', 'dim':'B flux', 'quant':'B flux', 'units':'Wb'}, 'floop_name':{'str':'flux_loop[chan].name'}, - 'floop_R':{'str':'flux_loop[chan].position.r', - 'dim':'distance', 'quant':'R', 'units':'m'}, - 'floop_Z':{'str':'flux_loop[chan].position.z', - 'dim':'distance', 'quant':'Z', 'units':'m'}}, + 'floop_R': {'str': 'flux_loop[chan].position.r', + 'dim': 'distance', 'quant': 'R', 'units': 'm'}, + 'floop_Z': {'str': 'flux_loop[chan].position.z', + 'dim': 'distance', 'quant': 'Z', 'units': 'm'}}, 'barometry': - {'t':{'str':'gauge[chan].pressure.time'}, - 'p':{'str':'gauge[chan].pressure.data', - 'dim':'pressure', 'quant':'p', 'units':'Pa?'}}, + {'t': {'str': 'gauge[chan].pressure.time'}, + 'names': {'str': 'gauge[chan].name'}, + 'p': {'str': 'gauge[chan].pressure.data', + 'dim': 'pressure', 'quant': 'p', 'units': 'Pa?'}}, 'neutron_diagnostic': {'t':{'str':'time', 'units':'s'}, @@ -282,6 +289,7 @@ class MultiIDSLoader(object): 'tau1keV':{'str':'channel[chan].optical_depth.data', 'dim':'optical_depth', 'quant':'tau', 'units':'adim.'}, 'validity_timed': {'str':'channel[chan].t_e.validity_timed'}, + 'names': {'str': 'channel[chan].name'}, 'Te0': {'str':'t_e_central.data', 'dim':'temperature', 'quant':'Te', 'units':'eV'}}, @@ -295,62 +303,82 @@ class MultiIDSLoader(object): 'dim':'distance', 'quant':'Z', 'units':'m'}, 'phi':{'str':'channel[chan].position.phi.data', 'dim':'angle', 'quant':'phi', 'units':'rad'}, - 'mode':{'str':'mode'}, - 'sweep':{'str':'sweep_time'}}, + 'names': {'str': 'channel[chan].name'}, + 'mode': {'str': 'mode'}, + 'sweep': {'str': 'sweep_time'}}, 'interferometer': - {'t':{'str':'time', - 'quant':'t', 'units':'s'}, - 'ne_integ':{'str':'channel[chan].n_e_line.data', - 'dim':'ne_integ', 'quant':'ne_integ', 'units':'/m2'}}, + {'t': {'str': 'time', + 'quant': 't', 'units': 's'}, + 'names': {'str': 'channel[chan].name'}, + 'ne_integ': {'str': 'channel[chan].n_e_line.data', + 'dim': 'ne_integ', 'quant': 'ne_integ', + 'units': '/m2', 'Brightness': True}}, 'polarimeter': - {'t':{'str':'time', - 'quant':'t', 'units':'s'}, - 'lamb':{'str':'channel[chan].wavelength', - 'dim':'distance', 'quant':'wavelength', 'units':'m'}, - 'fangle':{'str':'channel[chan].faraday_angle.data', - 'dim':'angle', 'quant':'faraday angle', 'units':'rad'}}, + {'t': {'str': 'time', + 'quant': 't', 'units': 's'}, + 'lamb': {'str': 'channel[chan].wavelength', + 'dim': 'distance', 'quant': 'wavelength', + 'units': 'm'}, + 'fangle': {'str': 'channel[chan].faraday_angle.data', + 'dim': 'angle', 'quant': 'faraday angle', + 'units': 'rad', 'Brightness': True}, + 'names': {'str': 'channel[chan].name'}}, 'bolometer': - {'t':{'str':'channel[chan].power.time', - 'quant':'t', 'units':'s'}, - 'power':{'str':'channel[chan].power.data', - 'dim':'power', 'quant':'power radiative', 'units':'W'}, - 'etendue':{'str':'channel[chan].etendue', - 'dim':'etendue', 'quant':'etendue', - 'units':'m2.sr'}, - 'tpower':{'str':'time','quant':'t', 'units':'s'}, - 'prad':{'str':'power_radiated_total', - 'dim':'power', 'quant':'power radiative', 'units':'W'}, - 'pradbulk':{'str':'power_radiated_inside_lcfs', - 'dim':'power', 'quant':'power radiative', - 'units':'W'}}, + {'t': {'str': 'channel[chan].power.time', + 'quant': 't', 'units': 's'}, + 'power': {'str': 'channel[chan].power.data', + 'dim': 'power', 'quant': 'power radiative', + 'units': 'W', 'Brightness': False}, + 'etendue': {'str': 'channel[chan].etendue', + 'dim': 'etendue', 'quant': 'etendue', + 'units': 'm2.sr'}, + 'names': {'str': 'channel[chan].name'}, + 'tpower': {'str': 'time', 'quant': 't', 'units': 's'}, + 'prad': {'str': 'power_radiated_total', + 'dim': 'power', 'quant': 'power radiative', + 'units': 'W'}, + 'pradbulk': {'str': 'power_radiated_inside_lcfs', + 'dim': 'power', 'quant': 'power radiative', + 'units': 'W'}}, 'soft_x_rays': - {'t':{'str':'time', - 'quant':'t', 'units':'s'}, - 'power':{'str':'channel[chan].power.data', - 'dim':'power', 'quant':'power radiative', 'units':'W'}, - 'brightness':{'str':'channel[chan].brightness.data', - 'dim':'brightness', 'quant':'brightness', 'units':'W/(m2.sr)'}, - 'etendue':{'str':'channel[chan].etendue', - 'dim':'etendue', 'quant':'etendue', 'units':'m2.sr'}}, + {'t': {'str': 'time', + 'quant': 't', 'units': 's'}, + 'power': {'str': 'channel[chan].power.data', + 'dim': 'power', 'quant': 'power radiative', + 'units': 'W', 'Brightness': False}, + 'brightness': {'str': 'channel[chan].brightness.data', + 'dim': 'brightness', 'quant': 'brightness', + 'units': 'W/(m2.sr)', 'Brightness': True}, + 'names': {'str': 'channel[chan].name'}, + 'etendue': {'str': 'channel[chan].etendue', + 'dim': 'etendue', 'quant': 'etendue', + 'units': 'm2.sr'}}, 'spectrometer_visible': {'t':{'str':'channel[chan].grating_spectrometer.radiance_spectral.time', 'quant':'t', 'units':'s'}, - 'spectra':{'str':'channel[chan].grating_spectrometer.radiance_spectral.data', - 'dim':'radiance_spectral', 'quant':'radiance_spectral', 'units':'ph/s/(m2.sr)/m'}, + 'spectra': {'str': ('channel[chan].grating_spectrometer' + + '.radiance_spectral.data'), + 'dim': 'radiance_spectral', + 'quant': 'radiance_spectral', + 'units': 'ph/s/(m2.sr)/m', 'Brightness': True}, + 'names': {'str': 'channel[chan].name'}, 'lamb':{'str':'channel[chan].grating_spectrometer.wavelengths', 'dim':'wavelength', 'quant':'wavelength', 'units':'m'}}, 'bremsstrahlung_visible': - {'t':{'str':'time', - 'quant':'t', 'units':'s'}, - 'radiance':{'str':'channel[chan].radiance_spectral.data', - 'dim':'radiance_spectral', 'quant':'radiance_spectral', - 'units':'ph/s/(m2.sr)/m'}, + {'t': {'str': 'time', + 'quant': 't', 'units': 's'}, + 'radiance': {'str': 'channel[chan].radiance_spectral.data', + 'dim': 'radiance_spectral', + 'quant': 'radiance_spectral', + 'units': 'ph/s/(m2.sr)/m', + 'Brightness': True}, + 'names': {'str': 'channel[chan].name'}, 'lamb_up': {'str':'channel[chan].filter.wavelength_upper'}, 'lamb_lo': {'str':'channel[chan].filter.wavelength_lower'}}, } @@ -479,11 +507,12 @@ class MultiIDSLoader(object): _icmod = lambda al, ar, axis=0: np.sum(al - ar, axis=axis) _eqB = lambda BT, BR, BZ: np.sqrt(BT**2 + BR**2 + BZ**2) def _rhopn1d(psi): - return np.sqrt( (psi - psi[:,0:1]) / (psi[:,-1] - psi[:,0])[:,None] ) + return np.sqrt((psi - psi[:, 0:1]) / (psi[:, -1] - psi[:, 0])[:, None]) def _rhopn2d(psi, psi0, psisep): - return np.sqrt( (psi - psi0[:,None]) / (psisep[:,None] - psi0[:,None]) ) + return np.sqrt( + (psi - psi0[:, None]) / (psisep[:, None] - psi0[:, None])) def _rhotn2d(phi): - return np.sqrt(phi / np.nanmax(phi, axis=1)[:,None]) + return np.sqrt(np.abs(phi) / np.nanmax(np.abs(phi), axis=1)[:, None]) def _eqSep(sepR, sepZ, npts=100): nt = len(sepR) @@ -642,7 +671,9 @@ def __init__(self, preset=None, dids=None, ids=None, occ=None, idd=None, idd = lidd[0] if len(lidd) > 0 else None self.add_ids(preset=preset, ids=ids, occ=occ, idd=idd, get=False) if ids_base is None: - ids_base = True + if not all([iids in self._IDS_BASE + for iids in self._dids.keys()]): + ids_base = True if ids_base is True: self.add_ids_base(get=False) if synthdiag is None: @@ -2218,77 +2249,173 @@ def _get_t0(self, t0=None): warnings.warn(msg) return t0 - - def to_Config(self, Name=None, occ=None, indDescription=None, plot=True): + def to_Config(self, Name=None, occ=None, + description_2d=None, mobile=None, plot=True): lidsok = ['wall'] - if indDescription is None: - indDescription = 0 + if description_2d is None: + description_2d = 0 # --------------------------- # Preliminary checks on data source consistency lids, lidd, shot, Exp = self._get_lidsidd_shotExp(lidsok, errshot=True, errExp=True, upper=True) - # ------------- - # Input dicts - - # config - config = None - if 'wall' in lids: - ids = 'wall' - - # occ = np.ndarray of valid int - occ = self._checkformat_getdata_occ(occ, ids) - assert occ.size == 1, "Please choose one occ only !" - occ = occ[0] - indoc = np.nonzero(self._dids[ids]['occ'] == occ)[0][0] - - wall = self._dids[ids]['ids'][indoc] - units = wall.description_2d[indDescription].limiter.unit - nunits = len(units) - - if nunits == 0: - msg = "There is no limiter unit stored !\n" - msg += "The required 2d description is empty:\n" - ms = "len(idd.%s[occ=%s].description_2d"%(ids,str(occ)) - msg += "%s[%s].limiter.unit) = 0"%(ms,str(indDescription)) + # ---------------- + # Trivial case + if 'wall' not in lids: + if plot: + msg = "ids 'wall' has not been loaded => Config not available!" raise Exception(msg) + return None - if Name is None: - Name = wall.description_2d[indDescription].type.name - if Name == '': - Name = 'imas wall' + # ---------------- + # Relevant case - import tofu.geom as mod + # Get relevant occ and description_2d + ids = 'wall' + # occ = np.ndarray of valid int + occ = self._checkformat_getdata_occ(occ, ids) + assert occ.size == 1, "Please choose one occ only !" + occ = occ[0] + indoc = np.nonzero(self._dids[ids]['occ'] == occ)[0][0] + + wall = self._dids[ids]['ids'][indoc].description_2d[description_2d] + kwargs = dict(Exp=Exp, Type='Tor') + + import tofu.geom as mod + + # Get vessel + nlim = len(wall.limiter.unit) + nmob = len(wall.mobile.unit) + # onelimonly = False + + # ---------------------------------- + # Relevant only if vessel is filled + # try: + # if len(wall.vessel.unit) != 1: + # msg = "There is no / several vessel.unit!" + # raise Exception(msg) + # if len(wall.vessel.unit[0].element) != 1: + # msg = "There is no / several vessel.unit[0].element!" + # raise Exception(msg) + # if len(wall.vessel.unit[0].element[0].outline.r) < 3: + # msg = "wall.vessel polygon has less than 3 points!" + # raise Exception(msg) + # name = wall.vessel.unit[0].element[0].name + # poly = np.array([wall.vessel.unit[0].element[0].outline.r, + # wall.vessel.unit[0].element[0].outline.z]) + # except Exception as err: + # # If vessel not in vessel, sometimes stored a a single limiter + # if nlim == 1: + # name = wall.limiter.unit[0].name + # poly = np.array([wall.limiter.unit[0].outline.r, + # wall.limiter.unit[0].outline.z]) + # onelimonly = True + # else: + # msg = ("There does not seem to be any vessel, " + # + "not in wall.vessel nor in wall.limiter!") + # raise Exception(msg) + # cls = None + # if name == '': + # name = 'ImasVessel' + # if '_' in name: + # ln = name.split('_') + # if len(ln) == 2: + # cls, name = ln + # else: + # name = name.replace('_', '') + # if cls is None: + # cls = 'Ves' + # assert cls in ['Ves', 'PlasmaDomain'] + # ves = getattr(mod, cls)(Poly=poly, Name=name, **kwargs) + + # Determine if mobile or not + lS = [] + # if onelimonly is False: + if mobile is None: + if nlim == 0 and nmob > 0: + mobile = True + elif nmob == 0 and nlim > 0: + mobile = False + elif nmob == nlim: + msgw = 'wall.description_2[{}]'.format(description_2d) + msg = ("\nids wall has same number of limiter / mobile units\n" + + "\t- len({}.limiter.unit) = {}\n".format(msgw, nlim) + + "\t- len({}.mobile.unit) = {}\n".format(msgw, nmob) + + " => Choosing limiter by default") + warnings.warn(msg) + mobile = False + else: + msgw = 'wall.description_2[{}]'.format(description_2d) + msg = ("Can't decide automatically whether to choose" + + " limiter or mobile!\n" + + "\t- len({}.limiter.unit) = {}\n".format(msgw, nlim) + + "\t- len({}.mobile.unit) = {}".format(msgw, nmob)) + raise Exception(msg) + assert isinstance(mobile, bool) - lS = [None for _ in units] - kwargs = dict(Exp=Exp, Type='Tor') - for ii in range(0,nunits): - poly = np.array([units[ii].outline.r, units[ii].outline.z]) + # Get PFC + if mobile is True: + units = wall.mobile.unit + else: + units = wall.limiter.unit + nunits = len(units) + + if nunits == 0: + msg = ("There is no unit stored !\n" + + "The required 2d description is empty:\n") + ms = "len(idd.{}[occ={}].description_2d".format(ids, occ) + msg += "{}[{}].limiter.unit) = 0".format(ms, + description_2d) + raise Exception(msg) + + lS = [None for _ in units] + for ii in range(0, nunits): + try: + if mobile is True: + outline = units[ii].outline[0] + else: + outline = units[ii].outline + poly = np.array([outline.r, outline.z]) if units[ii].phi_extensions.size > 0: - pos, extent = units[ii].phi_extensions.T + pos, extent = units[ii].phi_extensions.T else: pos, extent = None, None name = units[ii].name - cls = None + cls, mobi = None, None if name == '': name = 'unit{:02.0f}'.format(ii) if '_' in name: ln = name.split('_') if len(ln) == 2: cls, name = ln + elif len(ln) == 3: + cls, name, mobi = ln else: - name = name.replace('_','') + name = name.replace('_', '') if cls is None: - if ii == nunits-1: + if ii == nunits - 1: cls = 'Ves' else: cls = 'PFC' - lS[ii] = getattr(mod,cls)(Poly=poly, pos=pos, extent=extent, - Name=name, **kwargs) + mobi = mobi == 'mobile' + lS[ii] = getattr(mod, cls)(Poly=poly, pos=pos, + extent=extent, + Name=name, mobile=mobi, + **kwargs) + except Exception as err: + msg = ("PFC unit[{}] named {} ".format(ii, name) + + "could not be loaded!\n" + + str(err)) + raise Exception(msg) - config = mod.Config(lStruct=lS, Name=Name, **kwargs) + if Name is None: + Name = wall.type.name + if Name == '': + Name = 'imas wall' + + config = mod.Config(lStruct=lS, Name=Name, **kwargs) # Output if plot: @@ -2300,7 +2427,8 @@ def _checkformat_Plasma2D_dsig(self, dsig=None): lidsok = set(self._lidsplasma).intersection(self._dids.keys()) lscom = ['t'] - lsmesh = ['2dmeshNodes','2dmeshFaces'] + lsmesh = ['2dmeshNodes', '2dmeshFaces', + '2dmeshR', '2dmeshZ'] lc = [dsig is None, type(dsig) is str, @@ -2489,7 +2617,8 @@ def _checkformat_mesh_Rect(R, Z, datashape=None, assert R.ndim in [1, 2] assert Z.ndim in [1, 2] shapeu = np.unique(np.r_[R.shape, Z.shape]) - shapeRZ = [None, None] + if shapeRZ is None: + shapeRZ = [None, None] if R.ndim == 1: assert np.all(np.diff(R) > 0.) else: @@ -2526,37 +2655,25 @@ def _checkformat_mesh_Rect(R, Z, datashape=None, if datashape is not None: if None in shapeRZ: pass + shapeRZ = tuple(shapeRZ) - if shapeRZ == ['R', 'Z']: + if shapeRZ == ('R', 'Z'): assert datashape == (R.size, Z.size) - elif shapeRZ == ['Z', 'R']: + elif shapeRZ == ('Z', 'R'): assert datashape == (Z.size, R.size) else: msg = "Inconsistent data shape !" raise Exception(msg) - shapeRZ = tuple(shapeRZ) - assert shapeRZ in [('R', 'Z'), ('Z', 'R')] + if None not in shapeRZ: + shapeRZ = tuple(shapeRZ) + assert shapeRZ in [('R', 'Z'), ('Z', 'R')] return R, Z, shapeRZ, 0 - - - - - - - - - # TBF def get_mesh_from_ggd(path_to_ggd, ggdindex=0): pass - - - - - def _get_dextra(self, dextra=None, fordata=False, nan=True, pos=None): lc = [dextra == False, dextra is None, type(dextra) is str, type(dextra) is list, type(dextra) is dict] @@ -2841,7 +2958,7 @@ def to_Plasma2D(self, tlim=None, dsig=None, t0=None, npts, datashape = None, None keym = '{}.mesh'.format(ids) if cmesh else None for ss in set(out_.keys()).difference(lsigmesh): - assert out_[ss].ndim in [1,2] + assert out_[ss].ndim in [1, 2, 3] if out_[ss].ndim == 1: out_[ss] = np.atleast_2d(out_[ss]) shape = out_[ss].shape @@ -2864,10 +2981,10 @@ def to_Plasma2D(self, tlim=None, dsig=None, t0=None, shape = out_[ss].shape if len(shape) == 3: assert nt == shape[0] - datashape = tuple(shape[1], shape[2]) + datashape = (shape[1], shape[2]) if shapeRZ is None: - msg = ("Please provide shapeRZ" - + "indexing is ambiguous") + msg = ("Please provide shapeRZ," + + " indexing is ambiguous") raise Exception(msg) size = shape[1]*shape[2] if shapeRZ == ('R', 'Z'): @@ -2988,13 +3105,129 @@ def _checkformat_Cam_geom(self, ids=None, geomcls=None, indch=None): lgeom = [kk for kk in dir(tfg) if 'Cam' in kk] if geomcls not in [False] + lgeom: - msg = "Arg geomcls must be in %s"%str([False]+lgeom) + msg = "Arg geomcls must be in {}".format([False]+lgeom) raise Exception(msg) return geomcls + def _get_indch_geomtdata(self, indch=None, indch_auto=None, + dgeom=None, t=None, + ids=None, out=None, dsig=None, kk=None): + nch = 0 if indch is None else len(indch) + # Get from geometry of LOS consistency + if dgeom is not None: + indnan = np.logical_or(np.any(np.isnan(dgeom[0]), axis=0), + np.any(np.isnan(dgeom[1]), axis=0)) + if np.any(indnan) and not np.all(indnan): + if indch_auto is not True: + dmsg = {True: 'not available', False: 'ok'} + ls = ['index {} los {}'.format(ii, dmsg[indnan[ii]]) + for ii in range(0, dgeom[0].shape[1])] + msg = ("The geometry of all channels is not available !\n" + + "Please choose indch to get all LOS!\n" + + "Currently:\n" + + "\n ".join(ls) + + "\n\n => Solution: choose indch accordingly !") + raise Exception(msg) + else: + msg = ("Geometry missing for some los !\n" + + " => indch automatically set to:\n" + + " {}".format(indch)) + warnings.warn(msg) + if indch is None: + indch = (~indnan).nonzero()[0] + else: + indch = set(indch).intersection((~indnan).nonzero()[0]) + indch = np.array(indch, dtype=int) + + # Get from time vectors consistency + if t is not None: + if indch_auto is not True: + if indch is None: + ls = [ + 'index {} {}.shape {}'.format(ii, kk, + out[dsig[kk]][ii].shape) + for ii in range(0, len(out[dsig[kk]])) + ] + else: + ls = [ + 'index {} {}.shape {}'.format(indch[ii], kk, + out[dsig[kk]][ii].shape) + for ii in range(0, len(out[dsig[kk]])) + ] + msg = ("The following is supposed to be a np.ndarray:\n" + + " - diag: {}\n".format(ids) + + " - shortcut: {}\n".format(dsig[kk]) + + " - used as: {} input\n".format(kk) + + " Observed type: {}\n".format(type(out[dsig[kk]])) + + " Probable cause: non-uniform shape (vs channels)\n" + + " => shapes :\n " + + "\n ".join(ls) + + "\n => Solution: choose indch accordingly !") + raise Exception(msg) + ls = [t[ii].shape for ii in range(0, len(t))] + lsu = list(set([ssu for ssu in ls if 0 not in ssu])) + su = lsu[np.argmax([ls.count(ssu) for ssu in lsu])] + msg = ("indch set automatically for {}\n".format(ids) + + " (due to inhomogenous time shapes)\n" + + " - main shape: {}\n".format(su) + + " - nb. chan. selected: {}\n".format(len(indch)) + + " - indch: {}".format(indch)) + warnings.warn(msg) + + if indch is None: + indch = [ii for ii in range(0, len(t)) if ls[ii] == su] + else: + indch = [indch[ii] for ii in range(0, len(indch)) + if ls[indch[ii]] == su] + + # Get from data consistency + if all([ss is not None for ss in [out, kk, dsig]]): + if indch_auto: + ls = [out[dsig[kk]][ii].shape + for ii in range(0, len(out[dsig[kk]]))] + lsu = list(set([ssu for ssu in ls if 0 not in ssu])) + su = lsu[np.argmax([ls.count(ssu) for ssu in lsu])] + if indch is None: + indch = [ii for ii in range(0, len(out[dsig[kk]])) + if ls[ii] == su] + else: + indch = [indch[ii] for ii in range(0, len(out[dsig[kk]])) + if ls[ii] == su] + msg = ("indch set automatically for {}\n".format(ids) + + " (due to inhomogeneous data shapes)\n" + + " - main shape: {}\n".format(su) + + " - nb. chan. selected: {}\n".format(len(indch)) + + " - indch: {}".format(indch)) + warnings.warn(msg) + else: + if indch is None: + ls = [ + 'index {} {}.shape {}'.format(ii, kk, + out[dsig[kk]][ii].shape) + for ii in range(0, len(out[dsig[kk]])) + ] + else: + ls = [ + 'index {} {}.shape {}'.format(indch[ii], kk, + out[dsig[kk]][ii].shape) + for ii in range(0, len(out[dsig[kk]])) + ] + msg = ("The following is supposed to be a np.ndarray:\n" + + " - diag: {}\n".format(ids) + + " - shortcut: {}\n".format(dsig[kk]) + + " - used as: {} input\n".format(kk) + + " Observed type: {}\n".format(type(out[dsig[kk]])) + + " Probable cause: non-uniform shape (vs channels)\n" + + " => shapes :\n " + + "\n ".join(ls) + + "\n => Solution: choose indch accordingly !") + raise Exception(msg) + nchout = 0 if indch is None else len(indch) + return indch, nchout != nch + def _to_Cam_Du(self, ids, lk, indch, nan=None, pos=None): - Etendues, Surfaces = None, None + Etendues, Surfaces, names = None, None, None out = self.get_data(ids, sig=list(lk), indch=indch, nan=nan, pos=pos) if 'los_ptsRZPhi' in out.keys() and out['los_ptsRZPhi'].size > 0: @@ -3019,7 +3252,9 @@ def _to_Cam_Du(self, ids, lk, indch, nan=None, pos=None): Etendues = out['etendue'] if 'surface' in out.keys() and len(out['surface']) > 0: Surfaces = out['surface'] - return dgeom, Etendues, Surfaces + if 'names' in out.keys() and len(out['names']) > 0: + names = out['names'] + return dgeom, Etendues, Surfaces, names @@ -3043,10 +3278,9 @@ def to_Cam(self, ids=None, indch=None, indch_auto=False, config = self.to_Config(Name=Name, occ=occ, plot=False) # dchans + dchans = {} if indch is not None: - dchans = {'ind':indch} - else: - dchans = None + dchans['ind'] = indch # cam cam = None @@ -3057,35 +3291,26 @@ def to_Cam(self, ids=None, indch=None, indch_auto=False, raise Exception(msg) if 'LOS' in geom: - lk = ['los_ptsRZPhi','etendue','surface'] + lk = ['los_ptsRZPhi', 'etendue', 'surface', 'names'] lkok = set(self._dshort[ids].keys()) lkok = lkok.union(self._dcomp[ids].keys()) lk = list(set(lk).intersection(lkok)) - dgeom, Etendues, Surfaces = self._to_Cam_Du(ids, lk, indch, - nan=nan, pos=pos) - - indnan = np.logical_or(np.any(np.isnan(dgeom[0]),axis=0), - np.any(np.isnan(dgeom[1]),axis=0)) - if np.any(indnan) and not np.all(indnan): - indch_sug = (~indnan).nonzero()[0] - if indch_auto != True: - dmsg = {True: 'not available', False:'ok'} - msg = "The geometry of all channels is not available !\n" - msg += "Please choose indch to get all channels geomery !\n" - msg += "Currently:\n" - ls = ['index %s los %s'%(ii,dmsg[indnan[ii]]) - for ii in range(0,dgeom[0].shape[1])] - msg += "\n ".join(ls) - msg += "\n\n => Solution: choose indch accordingly !" - raise Exception(msg) - else: - indch = indch_sug - dgeom, Etendues, Surfaces = self._to_Cam_Du(ids, lk, indch, - nan=nan, pos=pos) - msg = "Geometry missing for some los !\n" - msg += " => indch automatically set to:\n" - msg += " %s"%str(indch) - warnings.warn(msg) + dgeom, Etendues, Surfaces, names = self._to_Cam_Du(ids, lk, indch, + nan=nan, + pos=pos) + + # Check all channels can be used, reset indch if necessary + indch, modif = self._get_indch_geomtdata(indch=indch, + indch_auto=indch_auto, + dgeom=dgeom) + if modif is True: + dgeom, Etendues, Surfaces, names = self._to_Cam_Du(ids, lk, + indch, + nan=nan, + pos=pos) + + if names is not None: + dchans['names'] = names import tofu.geom as tfg cam = getattr(tfg, geom)(dgeom=dgeom, config=config, @@ -3166,7 +3391,8 @@ def _checkformat_Data_dsig(self, ids=None, dsig=None, data=None, X=None, def to_Data(self, ids=None, dsig=None, data=None, X=None, tlim=None, indch=None, indch_auto=False, Name=None, occ=None, config=None, dextra=None, t0=None, datacls=None, geomcls=None, - plot=True, bck=True, fallback_X=None, nan=True, pos=None): + plot=True, bck=True, fallback_X=None, nan=True, pos=None, + return_indch=False): # dsig datacls, geomcls, dsig = self._checkformat_Data_dsig(ids, dsig, @@ -3193,63 +3419,29 @@ def to_Data(self, ids=None, dsig=None, data=None, X=None, tlim=None, else: dchans = None - # cam + # ----------- + # Get geom cam = None indchanstr = self._dshort[ids][dsig['data']]['str'].index('[chan]') chanstr = self._dshort[ids][dsig['data']]['str'][:indchanstr] nchMax = len(getattr(self._dids[ids]['ids'][0], chanstr)) + dgeom = None if geomcls != False: Etendues, Surfaces = None, None if config is None: msg = "A config must be provided to compute the geometry !" raise Exception(msg) - dgeom = None if 'LOS' in geomcls: - lk = ['los_ptsRZPhi','etendue','surface'] + lk_geom = ['los_ptsRZPhi', 'etendue', 'surface'] lkok = set(self._dshort[ids].keys()) lkok = lkok.union(self._dcomp[ids].keys()) - lk = list(set(lk).intersection(lkok)) - dgeom, Etendues, Surfaces = self._to_Cam_Du(ids, lk, indch, - nan=nan, pos=pos) - - indnan = np.logical_or(np.any(np.isnan(dgeom[0]),axis=0), - np.any(np.isnan(dgeom[1]),axis=0)) - if np.any(indnan) and not np.all(indnan): - indch_sug = (~indnan).nonzero()[0] - if indch_auto != True: - dmsg = {True: 'not available', False:'ok'} - msg = "The geometry of all channels is not available !\n" - msg += "Please de-activate geometry loading (geomcls=False)\n" - msg += " or choose indch to get all channels geometry !\n" - msg += "Currently:\n" - ls = ['index %s los %s'%(ii,dmsg[indnan[ii]]) - for ii in range(0,dgeom[0].shape[1])] - msg += "\n ".join(ls) - msg += "\n\n => Solution: choose indch accordingly !" - msg += " Suggested indch (los %s):\n"%dmsg[True] - msg += " %s"%str(indch_sug) - raise Exception(msg) - else: - indch = indch_sug - dgeom, Etendues, Surfaces = self._to_Cam_Du(ids, lk, indch, - nan=nan, pos=pos) - msg = "Geometry missing for some los !\n" - msg += " => indch automatically set to:\n" - msg += " %s"%str(indch) - warnings.warn(msg) - - if dgeom is not None: - import tofu.geom as tfg - cam = getattr(tfg, geomcls)(dgeom=dgeom, config=config, - Etendues=Etendues, Surfaces=Surfaces, - Name=Name, Diag=ids, Exp=Exp, - dchans=dchans) - cam.Id.set_dUSR( {'imas-nchMax': nchMax} ) + lk_geom = list(set(lk_geom).intersection(lkok)) + dgeom, Etendues, Surfaces, names = self._to_Cam_Du( + ids, lk_geom, indch, nan=nan, pos=pos) - - # ----------------------- - # data + # ---------- + # Get time lk = sorted(dsig.keys()) dins = dict.fromkeys(lk) t = self.get_data(ids, sig=dsig.get('t', 't'), indch=indch)['t'] @@ -3259,89 +3451,47 @@ def to_Data(self, ids=None, dsig=None, data=None, X=None, tlim=None, msg += " - 't' = %s"%str(t) raise Exception(msg) + # ----------- + # Check indch if type(t) is list: - if indch_auto == True: - ls = [t[ii].shape for ii in range(0,len(t))] - lsu = list(set([ssu for ssu in ls if 0 not in ssu])) - su = lsu[np.argmax([ls.count(ssu) for ssu in lsu])] - if indch is None: - indch = [ii for ii in range(0,len(t)) if ls[ii] == su] - else: - indchcam = [ii for ii in range(0,len(t)) if ls[ii] == su] - indch = [indch[ii] for ii in range(0,len(t)) if ls[ii] == su] - t = self.get_data(ids, sig='t', indch=indch)['t'] - if cam is not None: - cam = cam.get_subset(indch=indchcam) - msg = "indch set automatically for %s\n"%ids - msg += " (due to inhomogenous time shapes)\n" - msg += " - main shape: %s\n"%str(su) - msg += " - nb. chan. selected: %s\n"%len(indch) - msg += " - indch: %s"%str(indch) - warnings.warn(msg) - - else: - msg = "The time vector does not seem to be homogeneous !\n" - msg += "Please choose indch such that all channels have same t !\n" - msg += "Currently:\n" - if indch is None: - ls = ['index %s t.shape %s'%(ii,str(t[ii].shape)) - for ii in range(0,len(t))] - else: - ls = ['index %s t.shape %s'%(indch[ii],str(t[ii].shape)) - for ii in range(0,len(t))] - msg += "\n ".join(ls) - msg += "\n => Solution: choose indch accordingly !" - raise Exception(msg) + indch, modif = self._get_indch_geomtdata(indch=indch, + indch_auto=indch_auto, + dgeom=dgeom, t=t) + assert modif is True + else: + indch, modif = self._get_indch_geomtdata(indch=indch, + indch_auto=indch_auto, + dgeom=dgeom) + if modif is True: + if geomcls is not False: + dgeom, Etendues, Surfaces, names = self._to_Cam_Du( + ids, lk_geom, indch, nan=nan, pos=pos) + t = self.get_data(ids, sig='t', indch=indch)['t'] + modif = False + + if names is not None: + dchans['names'] = names if t.ndim == 2: - assert np.all(np.isclose(t, t[0:1,:])) - t = t[0,:] + assert np.all(np.isclose(t, t[0:1, :])) + t = t[0, :] dins['t'] = t indt = self._checkformat_tlim(t, tlim=tlim)['indt'] + # ----------- + # Get data out = self.get_data(ids, sig=[dsig[k] for k in lk], indt=indt, indch=indch, nan=nan, pos=pos) for kk in set(lk).difference('t'): if not isinstance(out[dsig[kk]], np.ndarray): - if indch_auto: - ls = [out[dsig[kk]][ii].shape - for ii in range(0,len(out[dsig[kk]]))] - lsu = list(set([ssu for ssu in ls if 0 not in ssu])) - su = lsu[np.argmax([ls.count(ssu) for ssu in lsu])] - if indch is None: - indch = [ii for ii in range(0,len(out[dsig[kk]])) - if ls[ii] == su] - else: - indch = [indch[ii] for ii in range(0,len(out[dsig[kk]])) - if ls[ii] == su] + indch, modifk = self._get_indch_geomtdata( + indch=indch, indch_auto=indch_auto, + out=out, dsig=dsig, kk=kk) + if modifk is True: out = self.get_data(ids, sig=[dsig[k] for k in lk], - indt=indt, indch=indch, nan=nan, - pos=pos) - if cam is not None: - cam = cam.get_subset(indch=indch) - msg = "indch set automatically for %s\n"%ids - msg += " (due to inhomogeneous data shapes)\n" - msg += " - main shape: %s\n"%str(su) - msg += " - nb. chan. selected: %s\n"%len(indch) - msg += " - indch: %s"%str(indch) - warnings.warn(msg) - else: - msg = "The following is supposed to be a np.ndarray:\n" - msg += " - diag: %s\n"%ids - msg += " - shortcut: %s\n"%dsig[kk] - msg += " - used as: %s input\n"%kk - msg += " Observed type: %s\n"%str(type(out[dsig[kk]])) - msg += " Probable cause: non-uniform shape (vs channels)\n" - msg += " => shapes :\n " - if indch is None: - ls = ['index %s %s.shape %s'%(ii,kk,str(out[dsig[kk]][ii].shape)) - for ii in range(0,len(out[dsig[kk]]))] - else: - ls = ['index %s %s.shape %s'%(indch[ii],kk,str(out[dsig[kk]][ii].shape)) - for ii in range(0,len(out[dsig[kk]]))] - msg += "\n ".join(ls) - msg += "\n => Solution: choose indch accordingly !" - raise Exception(msg) + indt=indt, indch=indch, + nan=nan, pos=pos) + modif = True # Arrange depending on shape and field if type(out[dsig[kk]]) is not np.ndarray: @@ -3350,7 +3500,7 @@ def to_Data(self, ids=None, dsig=None, data=None, X=None, tlim=None, ipdb.set_trace() raise Exception(msg) - assert out[dsig[kk]].ndim in [1,2,3] + assert out[dsig[kk]].ndim in [1, 2, 3] if out[dsig[kk]].ndim == 1: out[dsig[kk]] = np.atleast_2d(out[dsig[kk]]) @@ -3368,6 +3518,15 @@ def to_Data(self, ids=None, dsig=None, data=None, X=None, tlim=None, assert kk == 'data' dins[kk] = np.swapaxes(out[dsig[kk]].T, 1,2) + # Update dgeom if necessary + if modif is True and geomcls is not False: + dgeom, Etendues, Surfaces, names = self._to_Cam_Du( + ids, lk_geom, indch, + nan=nan, pos=pos) + modif = False + + # -------------------------- + # Format special ids cases if ids == 'reflectometer_profile': dins['X'] = np.fliplr(dins['X']) dins['data'] = np.fliplr(dins['data']) @@ -3382,7 +3541,6 @@ def to_Data(self, ids=None, dsig=None, data=None, X=None, tlim=None, fallback_X = 1.1*np.nanmax(dins['X']) dins['X'][np.isnan(dins['X'])] = fallback_X - # Apply indt if was not done in get_data for kk,vv in dins.items(): if (vv.ndim == 2 or kk == 't') and vv.shape[0] > indt.size: @@ -3408,17 +3566,29 @@ def to_Data(self, ids=None, dsig=None, data=None, X=None, tlim=None, for tt in dextra.keys(): dextra[tt]['t'] = dextra[tt]['t'] - t0 + # -------------- + # Create objects + if geomcls is not False and dgeom is not None: + import tofu.geom as tfg + cam = getattr(tfg, geomcls)(dgeom=dgeom, config=config, + Etendues=Etendues, Surfaces=Surfaces, + Name=Name, Diag=ids, Exp=Exp, + dchans=dchans) + cam.Id.set_dUSR({'imas-nchMax': nchMax}) + import tofu.data as tfd conf = None if cam is not None else config Data = getattr(tfd, datacls)(Name=Name, Diag=ids, Exp=Exp, shot=shot, lCam=cam, config=conf, dextra=dextra, dchans=dchans, **dins) - Data.Id.set_dUSR( {'imas-nchMax': nchMax} ) if plot: Data.plot(draw=True, bck=bck) - return Data + if return_indch is True: + return Data, indch + else: + return Data def _get_synth(self, ids, dsig=None, @@ -3488,13 +3658,30 @@ def _get_synth(self, ids, dsig=None, def calc_signal(self, ids=None, dsig=None, tlim=None, t=None, res=None, quant=None, ref1d=None, ref2d=None, q2dR=None, q2dPhi=None, q2dZ=None, - Brightness=None, interp_t=None, + Brightness=None, interp_t=None, newcalc=True, indch=None, indch_auto=False, Name=None, occ_cam=None, occ_plasma=None, config=None, dextra=None, t0=None, datacls=None, geomcls=None, bck=True, fallback_X=None, nan=True, pos=None, plot=True, plot_compare=None, plot_plasma=None): + # Check / format inputs + if plot is None: + plot = True + + if plot: + if plot_compare is None: + plot_compare = True + if plot_plasma is None: + plot_plasma = True + + # Get experimental data first if relevant + # to get correct indch for comparison + if plot and plot_compare: + data, indch = self.to_Data(ids, indch=indch, + indch_auto=indch_auto, t0=t0, + return_indch=True, plot=False) + # Get camera cam = self.to_Cam(ids=ids, indch=indch, Name=None, occ=occ_cam, config=config, @@ -3585,8 +3772,10 @@ def calc_signal(self, ids=None, dsig=None, tlim=None, t=None, res=None, # Calculate synthetic signal if Brightness is None: Brightness = self._didsdiag[ids]['synth'].get('Brightness', None) + dq['fill_value'] = 0. sig, units = cam.calc_signal_from_Plasma2D(plasma, res=res, t=t, Brightness=Brightness, + newcalc=newcalc, plot=False, **dq) sig._dextra = plasma.get_dextra(dextra) @@ -3594,16 +3783,35 @@ def calc_signal(self, ids=None, dsig=None, tlim=None, t=None, res=None, if ids == 'interferometer': sig = 2.*sig elif ids == 'polarimeter': - sig = 2.*sig + # For polarimeter, the vect is along the LOS + # it is not the direction of + sig = -2.*sig + + # Safety check regarding Brightness + _, _, dsig_exp = self._checkformat_Data_dsig(ids) + kdata = dsig_exp['data'] + B_exp = self._dshort[ids][kdata].get('Brightness', None) + err_comp = False + if Brightness != B_exp: + u_exp = self._dshort[ids][kdata].get('units') + msg = ("\nCalculated synthetic and chosen experimental data " + + "do not seem directly comparable !\n" + + "\t- chosen experimental data: " + + "{}, ({}), Brightness = {}\n".format(kdata, + u_exp, B_exp) + + "\t- calculated synthetic data: " + + "int({}), ({}), Brightness = {}\n".format(dq['quant'], + units, + Brightness) + + "\n => Consider changing data or Brigthness value") + err_comp = True + warnings.warn(msg) # plot if plot: - if plot_compare is None: - plot_compare = True - if plot_plasma is None: - plot_plasma = True if plot_compare: - data = self.to_Data(ids, indch=indch, t0=t0, plot=False) + if err_comp: + raise Exception(msg) sig._dlabels = data.dlabels data.plot_compare(sig) else: @@ -3624,7 +3832,7 @@ def calc_signal(self, ids=None, dsig=None, tlim=None, t=None, res=None, def load_Config(shot=None, run=None, user=None, tokamak=None, version=None, - Name=None, occ=0, indDescription=0, plot=True): + Name=None, occ=0, description_2d=0, plot=True): didd = MultiIDSLoader() didd.add_idd(shot=shot, run=run, @@ -3632,7 +3840,7 @@ def load_Config(shot=None, run=None, user=None, tokamak=None, version=None, didd.add_ids('wall', get=True) return didd.to_Config(Name=Name, occ=occ, - indDescription=indDescription, plot=plot) + description_2d=description_2d, plot=plot) # occ ? @@ -3842,10 +4050,10 @@ def _save_to_imas(obj, shot=None, run=None, refshot=None, refrun=None, occ=None, user=None, tokamak=None, version=None, dryrun=False, tfversion=None, verb=True, **kwdargs): - dfunc = {'Struct':_save_to_imas_Struct, - 'Config':_save_to_imas_Config, - 'CamLOS1D':_save_to_imas_CamLOS1D, - 'DataCam1D':_save_to_imas_DataCam1D} + dfunc = {'Struct': _save_to_imas_Struct, + 'Config': _save_to_imas_Config, + 'CamLOS1D': _save_to_imas_CamLOS1D, + 'DataCam1D': _save_to_imas_DataCam1D} # Preliminary check on object class @@ -3897,14 +4105,23 @@ def _save_to_imas(obj, shot=None, run=None, refshot=None, refrun=None, # Class-specific functions #-------------------------------- -def _save_to_imas_Struct( obj, +def _save_to_imas_Struct(obj, shot=None, run=None, refshot=None, refrun=None, occ=None, user=None, tokamak=None, version=None, dryrun=False, tfversion=None, verb=True, - description_2d=0, unit=0): + description_2d=None, description_typeindex=None, + unit=None): if occ is None: occ = 0 + if description_2d is None: + description_2d = 0 + if description_typeindex is None: + description_typeindex = 2 + description_typeindex = int(description_typeindex) + if unit is None: + unit = 0 + # Create or open IDS # ------------------ idd, shotfile = _open_create_idd(shot=shot, run=run, @@ -3918,8 +4135,19 @@ def _save_to_imas_Struct( obj, # data # -------- idd.wall.description_2d.resize( description_2d + 1 ) - idd.wall.description_2d[description_2d].limiter.unit.resize(1) - node = idd.wall.description_2d[description_2d].limiter.unit[0] + idd.wall.description_2d[description_2d].type.index = ( + description_typeindex) + idd.wall.description_2d[description_2d].type.name = ( + '{}_{}'.format(obj.__class__.__name__, obj.Id.Name)) + idd.wall.description_2d[description_2d].type.description = ( + "tofu-generated wall. Each PFC is represented independently as a" + + " closed polygon in tofu, which saves them as disjoint PFCs") + if obj._dgeom['mobile'] is True: + idd.wall.description_2d[description_2d].mobile.unit.resize(unit+1) + node = idd.wall.description_2d[description_2d].mobile.unit[unit] + else: + idd.wall.description_2d[description_2d].limiter.unit.resize(unit+1) + node = idd.wall.description_2d[description_2d].limiter.unit[unit] node.outline.r = obj._dgeom['Poly'][0,:] node.outline.z = obj._dgeom['Poly'][1,:] if obj.noccur > 0: @@ -3949,14 +4177,17 @@ def _save_to_imas_Struct( obj, err=err0, dryrun=dryrun, verb=verb) -def _save_to_imas_Config( obj, idd=None, shotfile=None, +def _save_to_imas_Config(obj, idd=None, shotfile=None, shot=None, run=None, refshot=None, refrun=None, occ=None, user=None, tokamak=None, version=None, dryrun=False, tfversion=None, close=True, verb=True, - description_2d=None): + description_2d=None, description_typeindex=None): if occ is None: occ = 0 + if description_2d is None: + description_2d = 0 + # Create or open IDS # ------------------ if idd is None: @@ -3974,20 +4205,22 @@ def _save_to_imas_Config( obj, idd=None, shotfile=None, nS = len(lS) if len(lclsIn) != 1: - msg = "One StructIn subclass is allowed / necessary !" + msg = "One (and only one) StructIn subclass is allowed / necessary !" raise Exception(msg) - if description_2d is None: + if description_typeindex is None: if nS == 1 and lcls[0] in ['Ves', 'PlasmaDomain']: - description_2d = 0 + description_typeindex = 0 else: - description_2d = 2 - assert description_2d in [0, 2] + description_typeindex = 1 + assert description_typeindex in [0, 1] - # Make sure StructIn is last (IMAS requirement) - ind = lcls.index(lclsIn[0]) - lS[-1], lS[ind] = lS[ind], lS[-1] + # Check whether there is any mobile element + ismobile = any([ss._dgeom['mobile'] for ss in lS]) + # Isolate StructIn and take out from lS + ves = lS.pop(lcls.index(lclsIn[0])) + nS = len(lS) # Fill in data # ------------------ @@ -3995,23 +4228,83 @@ def _save_to_imas_Config( obj, idd=None, shotfile=None, # data # -------- idd.wall.description_2d.resize( description_2d + 1 ) - idd.wall.description_2d[description_2d].type.name = obj.Id.Name - idd.wall.description_2d[description_2d].limiter.unit.resize(nS) - for ii in range(0,nS): - node = idd.wall.description_2d[description_2d].limiter.unit[ii] - node.outline.r = lS[ii].Poly_closed[0,:] - node.outline.z = lS[ii].Poly_closed[1,:] - if lS[ii].noccur > 0: - node.phi_extensions = np.array([lS[ii].pos, lS[ii].extent]).T - node.closed = True - node.name = '%s_%s'%(lS[ii].__class__.__name__, lS[ii].Id.Name) + wall = idd.wall.description_2d[description_2d] + wall.type.name = obj.Id.Name + wall.type.index = description_typeindex + wall.type.description = ( + "tofu-generated wall. Each PFC is represented independently as a" + + " closed polygon in tofu, which saves them as disjoint PFCs") + + # Fill limiter / mobile + if ismobile: + # resize nS + 1 for vessel + wall.mobile.unit.resize(nS + 1) + units = wall.mobile.unit + for ii in range(0, nS): + units[ii].outline.resize(1) + units[ii].outline[0].r = lS[ii].Poly[0, :] + units[ii].outline[0].z = lS[ii].Poly[1, :] + if lS[ii].noccur > 0: + units[ii].phi_extensions = np.array([lS[ii].pos, + lS[ii].extent]).T + units[ii].closed = True + name = '{}_{}'.format(lS[ii].__class__.__name__, + lS[ii].Id.Name) + if lS[ii]._dgeom['mobile'] is True: + name = name + '_mobile' + units[ii].name = name + else: + # resize nS + 1 for vessel + wall.limiter.unit.resize(nS + 1) + units = wall.limiter.unit + for ii in range(0, nS): + units[ii].outline.r = lS[ii].Poly[0, :] + units[ii].outline.z = lS[ii].Poly[1, :] + if lS[ii].noccur > 0: + units[ii].phi_extensions = np.array([lS[ii].pos, + lS[ii].extent]).T + units[ii].closed = True + name = '{}_{}'.format(lS[ii].__class__.__name__, + lS[ii].Id.Name) + if lS[ii]._dgeom['mobile'] is True: + name = name + '_mobile' + units[ii].name = name + + # Add Vessel at the end + ii = nS + if ismobile: + units[ii].outline.resize(1) + units[ii].outline[0].r = ves.Poly[0, :] + units[ii].outline[0].z = ves.Poly[1, :] + else: + units[ii].outline.r = ves.Poly[0, :] + units[ii].outline.z = ves.Poly[1, :] + units[ii].closed = True + units[ii].name = '{}_{}'.format(ves.__class__.__name__, + ves.Id.Name) + + # ---------------------------------- + # Fill vessel if needed + # vesname = '{}_{}'.format(ves.__class__.__name__, ves.Id.Name) + # wall.vessel.name = vesname + # wall.vessel.index = 1 + # wall.vessel.description = ( + # "tofu-generated vessel outline, with a unique unit / element") + + # wall.vessel.unit.resize(1) + # wall.vessel.unit[0].element.resize(1) + # element = wall.vessel.unit[0].element[0] + # element.name = vesname + # element.outline.r = ves.Poly[0, :] + # element.outline.z = ves.Poly[1, :] + # ---------------------------------- # IDS properties # -------------- com = "PFC contour generated:\n" - com += " - from %s"%obj.Id.SaveName - com += " - by tofu %s"%tfversion + com += " - from {}".format(obj.Id.SaveName) + com += " - by tofu {}".format(tfversion) _fill_idsproperties(idd.wall, com, tfversion) err0 = None diff --git a/tofu/utils.py b/tofu/utils.py index 4b6289428..7868d9b39 100644 --- a/tofu/utils.py +++ b/tofu/utils.py @@ -648,7 +648,7 @@ def _get_exception(q, ids, qtype='quantity'): def load_from_imas(shot=None, run=None, user=None, tokamak=None, version=None, ids=None, Name=None, returnas=None, tlim=None, config=None, - occ=None, indch=None, indDescription=None, equilibrium=None, + occ=None, indch=None, description_2d=None, equilibrium=None, dsig=None, data=None, X=None, t0=None, dextra=None, plot=True, plot_sig=None, plot_X=None, sharex=False, invertx=None, @@ -972,9 +972,9 @@ def load_from_imas(shot=None, run=None, user=None, tokamak=None, version=None, # export to instances for ii in range(0,nids): if returnas[ii] == 'Config': - dout[ss]['Config'].append(multi.to_Config(Name=Name, occ=occ, - indDescription=indDescription, - plot=False)) + dout[ss]['Config'].append(multi.to_Config( + Name=Name, occ=occ, + description_2d=description_2d, plot=False)) elif returnas[ii] == 'Plasma2D': dout[ss]['Plasma2D'].append(multi.to_Plasma2D(Name=Name, occ=occ, @@ -1037,7 +1037,7 @@ def load_from_imas(shot=None, run=None, user=None, tokamak=None, version=None, def calc_from_imas(shot=None, run=None, user=None, tokamak=None, version=None, ids=None, Name=None, out=None, tlim=None, config=None, - occ=None, indch=None, indDescription=None, equilibrium=None, + occ=None, indch=None, description_2d=None, equilibrium=None, dsig=None, data=None, X=None, t0=None, dextra=None, Brightness=None, res=None, interp_t=None, plot=True, plot_compare=True, sharex=False, diff --git a/tofu/version.py b/tofu/version.py index 4fc937bc2..7b58261e4 100644 --- a/tofu/version.py +++ b/tofu/version.py @@ -1,2 +1,2 @@ # Do not edit, pipeline versioning governed by git tags! -__version__ = '1.4.2-a5-48-gac0016b' +__version__ = '1.4.2-a5-91-g5e9b601a'