diff --git a/tofu/geom/_core.py b/tofu/geom/_core.py index 1bd4a0f98..1f5fd37e0 100644 --- a/tofu/geom/_core.py +++ b/tofu/geom/_core.py @@ -3074,37 +3074,21 @@ def _get_phithetaproj_dist( return dist, indStruct - def plot_phithetaproj_dist( - self, - refpt=None, - ntheta=None, - nphi=None, - theta=None, - phi=None, - cmap=None, - ax=None, - fs=None, - tit=None, - wintit=None, - draw=None, - ): - dist, indStruct = self._get_phithetaproj_dist( - refpt=refpt, ntheta=ntheta, nphi=nphi, theta=theta, phi=phi - ) - return _plot.Config_phithetaproj_dist( - self, - refpt, - dist, - indStruct, - cmap=cmap, - ax=ax, - fs=fs, - tit=tit, - wintit=wintit, - draw=draw, - ) + def plot_phithetaproj_dist(self, refpt=None, ntheta=None, nphi=None, + theta=None, phi=None, cmap=None, invertx=None, + ax=None, fs=None, tit=None, wintit=None, + draw=None): + dist, indStruct = self._get_phithetaproj_dist(refpt=refpt, + ntheta=ntheta, nphi=nphi, + theta=theta, phi=phi) + return _plot.Config_phithetaproj_dist(self, refpt, dist, indStruct, + cmap=cmap, ax=ax, fs=fs, + tit=tit, wintit=wintit, + invertx=invertx, draw=draw) + def isInside(self, pts, In="(X,Y,Z)", log="any"): + """ Return a 2D array of bool Equivalent to applying isInside to each Struct diff --git a/tofu/geom/_plot.py b/tofu/geom/_plot.py index bba7a02e3..88235b1a3 100644 --- a/tofu/geom/_plot.py +++ b/tofu/geom/_plot.py @@ -692,13 +692,15 @@ def Plot_Impact_3DPoly(T, Leg="", ax=None, Ang=_def.TorPAng, def Config_phithetaproj_dist(config, refpt, dist, indStruct, distonly=False, - cmap=None, vmin=None, vmax=None, + cmap=None, vmin=None, vmax=None, invertx=None, ax=None, fs=None, cbck=(0.8,0.8,0.8,0.8), tit=None, wintit=None, legend=None, draw=None): if cmap is None: cmap = 'touch' lS = config.lStruct indsu = np.unique(indStruct) + if invertx is None: + invertx = True # set extent ratio = refpt[0] / np.nanmin(dist) @@ -745,6 +747,8 @@ def Config_phithetaproj_dist(config, refpt, dist, indStruct, # labels.append( '%s_%s'%(lS[ii].Id.Cls, lS[ii].Id.Name) ) # dax['cross'][0].legend(handles, labels, frameon=False, # bbox_to_anchor=(1.01,1.), loc=2, borderaxespad=0.) + if invertx is True: + dax['dist'][0].invert_xaxis() if draw: fig.canvas.draw() diff --git a/tofu/mag/equimap.py b/tofu/mag/equimap.py index 5b063e037..e9ae396df 100644 --- a/tofu/mag/equimap.py +++ b/tofu/mag/equimap.py @@ -20,19 +20,19 @@ import imas_west except ImportError as err: pass -# print(err) + print(err) try: import pywed as pw except ImportError as err: pass -# print(err) + print(err) # Project modules try: import mag_ripple as mr except ImportError as err: pass -# print(err) + print(err) __all__ = ['get'] diff --git a/tofu/mag/magFieldLines.py b/tofu/mag/magFieldLines.py index a50999e29..39c495437 100644 --- a/tofu/mag/magFieldLines.py +++ b/tofu/mag/magFieldLines.py @@ -160,7 +160,7 @@ def __init__(self, shot, run=0, occ=0, user='imas_public', machine='west'): bounds_error=False) def trace_mline(self, init_state, time, direction='FWD', \ - length_line=None, stp=None): + length_line=None, stp=None, ripple=True): ''' Traces the field line given a starting point. Integration step defined by stp and maximum length of the field line @@ -174,6 +174,7 @@ def trace_mline(self, init_state, time, direction='FWD', \ - time : array of times requested (list or numpy array) - direction : direction of integration 'FWD' forward or 'REV' reverse (string) + - ripple : take into account magnetic ripple or not (boolean) output: returns a dictionary containing: @@ -209,7 +210,8 @@ def trace_mline(self, init_state, time, direction='FWD', \ # !!!!!!!!!!!!!!!!!!!!! # HARD CODED CORRECTION - bt_intp_t *= -1 + if (np.nanmean(bt_intp_t) > 0): + bt_intp_t *= -1 # !!!!!!!!!!!!!!!!!!!!! # Interpolate current @@ -236,7 +238,7 @@ def trace_mline(self, init_state, time, direction='FWD', \ for jj in range(init_state.shape[1]): out_prime = self.integrate_solve_ivp(init_state[:, jj], \ - direction, s, ds) + direction, s, ds, ripple) out_prime['init_point'] = init_state[:, jj] out_prime['time'] = self.ar_time[ii] # Return list of dict @@ -246,7 +248,7 @@ def trace_mline(self, init_state, time, direction='FWD', \ return outMagLine - def integrate_solve_ivp(self, init_state, direction, s, ds): + def integrate_solve_ivp(self, init_state, direction, s, ds, ripple=True): ''' Integration step defined by stp and maximum length of the field line defined by s. @@ -258,6 +260,7 @@ def integrate_solve_ivp(self, init_state, direction, s, ds): (string) - s : field line maximum length - ds : vector of steps + - ripple : boolean, take into account magnetic ripple or not output: returns a dictionary containing: @@ -268,14 +271,26 @@ def integrate_solve_ivp(self, init_state, direction, s, ds): - y : y coordinate (np.array) - cp : collision point with the wall (list) ''' - if (direction=='FWD'): - sol = spode.solve_ivp(self.mfld3dcylfwd, [0, s], init_state, - method='RK23', t_eval=ds, - events=self.hit_wall_circ) - elif (direction=='REV'): - sol = spode.solve_ivp(self.mfld3dcylrev, [0, s], init_state, - method='RK23', t_eval=ds, - events=self.hit_wall_circ) + if (ripple): + if direction == 'FWD': + sol = spode.solve_ivp(self.mfld3dcylfwd, [0, s], init_state, + method='RK23', t_eval=ds, + events=self.hit_wall_circ) + elif direction == 'REV': + sol = spode.solve_ivp(self.mfld3dcylrev, [0, s], init_state, + method='RK23', t_eval=ds, + events=self.hit_wall_circ) + else: + if direction == 'FWD': + sol = spode.solve_ivp(self.mfld3dcylfwd_no_ripple, + [0, s], init_state, + method='RK23', t_eval=ds, + events=self.hit_wall_circ) + elif direction == 'REV': + sol = spode.solve_ivp(self.mfld3dcylrev_no_ripple, + [0, s], init_state, + method='RK23', t_eval=ds, + events=self.hit_wall_circ) sgf = sol.t rgf = sol.y[0] pgf = sol.y[1] @@ -321,9 +336,43 @@ def mfld3dcylrev(self, s, state): ''' R, P, Z = state Br, Bt, Bz = self.b_field_interp(R, P, Z) - #Br=self.fBp_r(R,Z)[0]+self.fBt_r(P,self.ftheta(R,Z),R) - #Bt=self.fBt_t(P,self.ftheta(R,Z),R) - #Bz=self.fBp_z(R,Z)[0] + # Br=self.fBp_r(R,Z)[0]+self.fBt_r(P,self.ftheta(R,Z),R) + # Bt=self.fBt_t(P,self.ftheta(R,Z),R) + # Bz=self.fBp_z(R,Z)[0] + B = np.sqrt(Br*Br + Bz*Bz + Bt*Bt) + d_R = -Br/B + d_P = -Bt/B*1/R + d_Z = -Bz/B + return [d_R, d_P, d_Z] + + def mfld3dcylfwd_no_ripple(self, s, state): + ''' + Returns the right end side of the field line system of equations in + cyclindrical (R, Phi, Z) coord. + This is the case for forward integration. + ''' + R, P, Z = state + Br, Bt, Bz = self.b_field_interp_no_ripple(R, P, Z) + # Br=self.fBp_r(R,Z)[0]+self.fBt_r(P,self.ftheta(R,Z),R) + # Bt=self.fBt_t(P,self.ftheta(R,Z),R) + # Bz=self.fBp_z(R,Z)[0] + B = np.sqrt(Br*Br + Bz*Bz + Bt*Bt) + d_R = Br/B + d_P = Bt/B*1/R + d_Z = Bz/B + return [d_R, d_P, d_Z] + + def mfld3dcylrev_no_ripple(self, s, state): + ''' + Returns the right end side of the field line system of equations in + cyclindrical (R, Phi, Z) coord. + This is the case for backward integration. + ''' + R, P, Z = state + Br, Bt, Bz = self.b_field_interp_no_ripple(R, P, Z) + # Br=self.fBp_r(R,Z)[0]+self.fBt_r(P,self.ftheta(R,Z),R) + # Bt=self.fBt_t(P,self.ftheta(R,Z),R) + # Bz=self.fBp_z(R,Z)[0] B = np.sqrt(Br*Br + Bz*Bz + Bt*Bt) d_R = -Br/B d_P = -Bt/B*1/R @@ -340,14 +389,17 @@ def hit_wall_circ(self, s, state): mat file ''' R, P, Z = state - #if np.abs(Z)<=0.5 and R>3.01 and np.deg2rad(89.5)<=P<=np.deg2rad(90.5): - # return 0 - #if np.abs(Z)<=0.5 and R>3.01 and np.deg2rad(179.5)<=P<=np.deg2rad(180.5): - # print('Got 2') - # return 0 - #elif np.abs(Z)<=0.5 and R>3.01 and np.deg2rad(269.5)<=P<=np.deg2rad(270.5): - # print('Got 3') - # return 0 + # if (np.abs(Z) <= 0.5 and R > 3.01 + # and np.deg2rad(89.5) <= P <= np.deg2rad(90.5)): + # return 0 + # if (np.abs(Z) <= 0.5 and R > 3.01 + # and np.deg2rad(179.5) <= P <= np.deg2rad(180.5)): + # print('Got 2') + # return 0 + # elif (np.abs(Z) <= 0.5 and R > 3.01 + # and np.deg2rad(269.5) <= P <= np.deg2rad(270.5)): + # print('Got 3') + # return 0 if self.wall_ck==False: Rc=2.460 Zc=0.0 @@ -362,10 +414,9 @@ def hit_wall_circ(self, s, state): hit_wall_circ.terminal=True - def b_field_interp(self, R, Phi, Z, no_ripple=False): + def b_field_interp(self, R, Phi, Z): ''' Linear interpolation of B vector components at R, Phi, Z positions ''' - interp_points = np.vstack((R, Z)).transpose() br_intp = self.br_lin_intp.__call__(interp_points) @@ -384,6 +435,17 @@ def b_field_interp(self, R, Phi, Z, no_ripple=False): return br_intp[0], bt_intp[0], bz_intp[0] + def b_field_interp_no_ripple(self, R, Phi, Z): + ''' Linear interpolation of B vector components at R, Phi, Z positions + ''' + interp_points = np.vstack((R, Z)).transpose() + + br_intp = self.br_lin_intp.__call__(interp_points) + bt_intp = self.bt_lin_intp.__call__(interp_points) + bz_intp = self.bz_lin_intp.__call__(interp_points) + + return br_intp[0], bt_intp[0], bz_intp[0] + def plot_trace(self, trace): ''' Plots a summary of the magnetic field line trace. @@ -511,3 +573,27 @@ def plot_trace_3D(self, trace): ax.set_xlabel('x [m]'); ax.set_ylabel('y [m]'); ax.set_zlabel('z [m]') plt.title('Magnetic field line trace') plt.show() + + def plot_poincare(self, trace_plt, tor_angle=0): + ''' + Plots Poincare plot at given toroidal angle + + input: + - trace_plt : the magneticl field line data for several points + one single time (dictionary) + ''' + + fig, ax = plt.subplots() + + for ii in range(len(trace_plt)): + loc_cross = trace_plt[ii]['p'] % (tor_angle + 2*np.pi) + mask_cross = np.r_[False, loc_cross[1:] + < loc_cross[:-1]] & np.r_[loc_cross[:-1] + < loc_cross[1:], + False] + ax.plot(trace_plt[ii]['r'][mask_cross], + trace_plt[ii]['z'][mask_cross], 'o', markersize=3) + + ax.set_aspect('equal') + ax.set_xlabel('R [m]') + ax.set_ylabel('Z [m]') diff --git a/tofu/mag/mag_ripple/__init__.py b/tofu/mag/mag_ripple/__init__.py index bc5c08174..7575a4fde 100644 --- a/tofu/mag/mag_ripple/__init__.py +++ b/tofu/mag/mag_ripple/__init__.py @@ -3,3 +3,20 @@ ''' Magnetic ripple ''' + +import warnings +import traceback + +try: + try: + from tofu.mag.mag_ripple import * + except Exception: + from .mag_ripple import * + del warnings, traceback +except Exception as err: + msg = str(traceback.format_exc()) + msg += "\n\n => the optional sub-package is not usable\n" + warnings.warn(msg) + del msg, err + +__all__ = ['mag_ripple'] diff --git a/tofu/mag/test_equimap.py b/tofu/mag/test_equimap.py index a88b114b5..08755ade6 100644 --- a/tofu/mag/test_equimap.py +++ b/tofu/mag/test_equimap.py @@ -11,288 +11,289 @@ import os import sys -#print('path 1 =', sys.path) -sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), '..'))) -#print('path 2 =', sys.path) - -# Local modules -import equimap -import imas - -interp_points = 60 -shot = 53221 - -approx_time_in = 39.9881 -eps_time = 1E-3 -ctr_plot = [0] - -# CALL EQUIMAP -oute = equimap.get(shot, time=[35.9249267578125, 39.98], \ - R=[2., 2.2, 2.4, 2.7, 3.], Phi=[0, 0, 0, 0, 0], Z=[0, 0, 0, 0, 0], \ - quantity='psi') - -print('') -print('oute.shape =', oute.shape) -# oute should be: -# array([[-2.27165658, -3.01387112, -3.98753036, -3.92317787, -2.43752966], -# [ 1.09062263, 0.34659524, -0.63331812, -0.47944498, 1.01023176]]) -print('oute =', oute) -print('') - -idd = imas.ids(shot, 0) -idd.open_env('imas_public', 'west', '3') -idd.equilibrium.get() -out = idd.equilibrium - -equiDict = {} - -# Declaration of arrays 2d plots -equi_grid = idd.equilibrium.grids_ggd[0].grid[0] -NbrPoints = len(equi_grid.space[0].objects_per_dimension[0].object) -equiDict['r'] = np.full(NbrPoints, np.nan) -equiDict['z'] = np.full(NbrPoints, np.nan) -for ii in range(NbrPoints): - equiDict['r'][ii] = equi_grid.space[0].objects_per_dimension[0]. \ - object[ii].geometry[0] - equiDict['z'][ii] = equi_grid.space[0].objects_per_dimension[0]. \ - object[ii].geometry[1] - -arg_time = np.argmin(np.abs(out.time - approx_time_in)) - -print('\nDistance between requested and calculated time =', \ - np.abs(out.time[arg_time] - approx_time_in), '\n') - -time_in = out.time[arg_time] - eps_time - -print('\nNew distance between requested and calculated time =', \ - np.abs(out.time[arg_time] - time_in), '\n', \ - 'with eps_time =', eps_time, '\n') - -R_all = np.linspace(np.min(equiDict['r']), np.max(equiDict['r']), interp_points) -Z_all = np.linspace(np.min(equiDict['z']), np.max(equiDict['z']), interp_points) - -R_all_tot = np.repeat(R_all, interp_points) -Z_all_tot = np.tile(Z_all, interp_points) - -# CALL EQUIMAP -outa = equimap.get(shot, time=time_in, \ - R=R_all_tot, Phi=np.zeros(R_all_tot.shape), Z=Z_all_tot, \ - quantity='psi') - -outar = outa.reshape((interp_points, interp_points)) -Rr = R_all_tot.reshape((interp_points, interp_points)) -Zr = Z_all_tot.reshape((interp_points, interp_points)) - -print('') -print('outa.shape =', outa.shape) -print('outar.shape =', outar.shape) -#print('oute =', oute) -print('') - -#plt.figure() -#cs_ref = plt.contour(out.interp2D.r, out.interp2D.z, \ -# out.interp2D.psi[np.argmin(np.abs(out.time - time_in))], \ -# ctr_plot, colors='k', linestyle='solid') -#cs_ref_tri = plt.tricontour(equiDict['r'], equiDict['z'], \ -# out.ggd[0].grid.space[0].objects_per_dimension[1].nodes, \ -# out.ggd[0].psi[np.argmin(np.abs(out.time - time_in))], \ -# ctr_plot, colors='red', linestyles='dashdot') -#cs = plt.contour(Rr, Zr, outar, ctr_plot, \ -# colors='orange', linestyles='dashed') -#plt.plot(np.nan, np.nan, color='k', linestyle='solid', label='ref interp2D') -#plt.plot(np.nan, np.nan, color='red', linestyle='dashdot', label='ref ggd') -#plt.plot(np.nan, np.nan, color='orange', linestyle='dashed', label='test interp') -#plt.legend() -#plt.xlabel('R [m]') -#plt.ylabel('Z [m]') -#plt.title('Contour psi, level ' + str(ctr_plot) + ' at time ' + str(time_in) + 's') -#plt.axes().set(aspect=1) -# -#p = cs.collections[0].get_paths()[0] -#v = p.vertices -#x_cs = v[:, 0] -#y_cs = v[:, 1] - -#p = cs_ref.collections[0].get_paths()[0] -#v = p.vertices -#x_cs_ref = v[:, 0] -#y_cs_ref = v[:, 1] - -plt.figure() -plt.contourf(Rr, Zr, outar) -plt.plot(np.squeeze(out.time_slice[arg_time].boundary.outline.r), \ - np.squeeze(out.time_slice[arg_time].boundary.outline.z), \ - linewidth=2, color='red') -plt.plot(out.time_slice[arg_time].global_quantities.magnetic_axis.r, \ - out.time_slice[arg_time].global_quantities.magnetic_axis.z, \ - marker='+', color='red', markersize=20) -plt.colorbar() -plt.xlabel('R [m]') -plt.ylabel('Z [m]') -plt.title('psi interpolated at time ' + str(time_in) + 's') -plt.axes().set(aspect=1) - -#plt.figure() -#plt.plot(x_cs_ref, y_cs_ref, label='Ref') -#plt.plot(x_cs, y_cs, '--', label='Interp') -#plt.xlabel('R [m]') -#plt.ylabel('Z [m]') -#plt.legend() - -#pts_cs_ref = np.linspace(0, 1, x_cs_ref.shape[0]) -#pts_cs = np.linspace(0, 1, x_cs.shape[0]) - -#x_cs_intp = np.interp(pts_cs_ref, pts_cs, x_cs) -#y_cs_intp = np.interp(pts_cs_ref, pts_cs, y_cs) - -#error_x = np.abs((x_cs_intp - x_cs_ref)/x_cs_ref)*100 -#error_y = np.abs((y_cs_intp - y_cs_ref)/y_cs_ref)*100 - -#plt.figure() -#plt.plot(error_x, label='err. x') -#plt.plot(error_y, label='err. y') -#plt.legend() -#plt.ylabel('Error interp psi contour level ' + str(ctr_plot) + ' in [%]' \ -# + '\nwith time difference, eps_time = ' + str(eps_time*100) + '%', \ -# multialignment='center') - -# CALL EQUIMAP -out_rho_pol_norm = equimap.get(shot, time=time_in, \ - R=R_all_tot, Phi=np.zeros(R_all_tot.shape), Z=Z_all_tot, \ - quantity='rho_pol_norm') - -print('') -print('out_rho_pol_norm.shape =', out_rho_pol_norm.shape) -#print('out_rho_pol_norm =', out_rho_pol_norm) -print('') - -plt.figure() -out_rho_pol_norm_r = out_rho_pol_norm.reshape((interp_points, interp_points)) -cs = plt.contourf(Rr, Zr, out_rho_pol_norm_r) -plt.colorbar() -plt.plot(np.squeeze(out.time_slice[arg_time].boundary.outline.r), \ - np.squeeze(out.time_slice[arg_time].boundary.outline.z), \ - linewidth=2, color='red') -plt.plot(out.time_slice[arg_time].global_quantities.magnetic_axis.r, \ - out.time_slice[arg_time].global_quantities.magnetic_axis.z, \ - marker='+', color='red', markersize=20) -plt.xlabel('R [m]') -plt.ylabel('Z [m]') -plt.title('rho_pol_norm') - -# CALL EQUIMAP -out_rho_tor_norm = equimap.get(shot, time=time_in, \ - R=R_all_tot, Phi=np.zeros(R_all_tot.shape), Z=Z_all_tot, \ - quantity='rho_tor_norm') -print('') -print('out_rho_tor_norm.shape =', out_rho_tor_norm.shape) -#print('out_rho_tor_norm =', out_rho_tor_norm) -print('') -plt.figure() -out_rho_tor_norm_r = out_rho_tor_norm.reshape((interp_points, interp_points)) -cs = plt.contourf(Rr, Zr, out_rho_tor_norm_r) -plt.colorbar() -arg_time = np.argmin(np.abs(out.time - time_in)) -plt.plot(np.squeeze(out.time_slice[arg_time].boundary.outline.r), \ - np.squeeze(out.time_slice[arg_time].boundary.outline.z), \ - linewidth=2, color='red') -plt.plot(out.time_slice[arg_time].global_quantities.magnetic_axis.r, \ - out.time_slice[arg_time].global_quantities.magnetic_axis.z, \ - marker='+', color='red', markersize=20) -plt.xlabel('R [m]') -plt.ylabel('Z [m]') -plt.title('rho_tor_norm') - -# CALL EQUIMAP -out_rho_tor = equimap.get(shot, time=time_in, \ - R=R_all_tot, Phi=np.zeros(R_all_tot.shape), Z=Z_all_tot, \ - quantity='rho_tor') - -print('') -print('out_rho_tor.shape =', out_rho_tor.shape) -#print('out_rho_tor =', out_rho_tor) -print('') -plt.figure() -out_rho_tor_r = out_rho_tor.reshape((interp_points, interp_points)) -cs = plt.contourf(Rr, Zr, out_rho_tor_r) -plt.colorbar() -arg_time = np.argmin(np.abs(out.time - time_in)) -plt.plot(np.squeeze(out.time_slice[arg_time].boundary.outline.r), \ - np.squeeze(out.time_slice[arg_time].boundary.outline.z), \ - linewidth=2, color='red') -plt.plot(out.time_slice[arg_time].global_quantities.magnetic_axis.r, \ - out.time_slice[arg_time].global_quantities.magnetic_axis.z, \ - marker='+', color='red', markersize=20) -plt.xlabel('R [m]') -plt.ylabel('Z [m]') -plt.title('rho_tor') - -plt.figure() -plt.plot(out.time_slice[arg_time].profiles_1d.psi, out.time_slice[arg_time].profiles_1d.rho_tor) -plt.xlabel('Psi') -plt.ylabel('rho_tor') - -#plt.figure() -#plt.tricontourf(equiDict['r'], equiDict['z'], \ -# out.ggd[0].grid.space[0].objects_per_dimension[1].nodes, \ -# out.ggd[0].psi[arg_time]) -#plt.plot(np.squeeze(out.boundary.outline.r[arg_time]), \ -# np.squeeze(out.boundary.outline.z[arg_time]), \ -# linewidth=2, color='red') -#plt.plot(out.global_quantities.magnetic_axis.r[arg_time], \ -# out.global_quantities.magnetic_axis.z[arg_time], \ -# marker='+', color='red', markersize=20) -#plt.colorbar() -#plt.xlabel('R [m]') -#plt.ylabel('Z [m]') -#plt.title('psi ggd') - -# CALL EQUIMAP -out_b_norm = equimap.get(shot, time=time_in, \ - R=R_all_tot, Phi=np.zeros(R_all_tot.shape), Z=Z_all_tot, \ - quantity='b_field_norm') - -print('') -print('out_b_norm.shape =', out_b_norm.shape) -#print('out_b_norm =', out_b_norm) -print('') -plt.figure() -out_b_norm_r = out_b_norm.reshape((interp_points, interp_points)) -cs = plt.contourf(Rr, Zr, out_b_norm_r) -plt.colorbar() -arg_time = np.argmin(np.abs(out.time - time_in)) -plt.plot(np.squeeze(out.time_slice[arg_time].boundary.outline.r), \ - np.squeeze(out.time_slice[arg_time].boundary.outline.z), \ - linewidth=2, color='red') -plt.plot(out.time_slice[arg_time].global_quantities.magnetic_axis.r, \ - out.time_slice[arg_time].global_quantities.magnetic_axis.z, \ - marker='+', color='red', markersize=20) -plt.xlabel('R [m]') -plt.ylabel('Z [m]') -plt.title('b_norm') - -# CALL EQUIMAP -out_theta = equimap.get(shot, time=time_in, \ - R=R_all_tot, Phi=np.zeros(R_all_tot.shape), Z=Z_all_tot, \ - quantity='theta') - -print('') -print('out_theta.shape =', out_theta.shape) -print('out_theta =', out_theta) -print('') -plt.figure() -out_theta_r = out_theta.reshape((interp_points, interp_points)) -cs = plt.contourf(Rr, Zr, out_theta_r, 100) -plt.colorbar() -arg_time = np.argmin(np.abs(out.time - time_in)) -plt.plot(np.squeeze(out.time_slice[arg_time].boundary.outline.r), \ - np.squeeze(out.time_slice[arg_time].boundary.outline.z), \ - linewidth=2, color='red') -plt.plot(out.time_slice[arg_time].global_quantities.magnetic_axis.r, \ - out.time_slice[arg_time].global_quantities.magnetic_axis.z, \ - marker='+', color='red', markersize=20) -plt.xlabel('R [m]') -plt.ylabel('Z [m]') -plt.title('theta') - -plt.show() +if __name__ == '__main__': + #print('path 1 =', sys.path) + sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), '..'))) + #print('path 2 =', sys.path) + + # Local modules + import equimap + import imas + + interp_points = 60 + shot = 53221 + + approx_time_in = 39.9881 + eps_time = 1E-3 + ctr_plot = [0] + + # CALL EQUIMAP + oute = equimap.get(shot, time=[35.9249267578125, 39.98], \ + R=[2., 2.2, 2.4, 2.7, 3.], Phi=[0, 0, 0, 0, 0], Z=[0, 0, 0, 0, 0], \ + quantity='psi') + + print('') + print('oute.shape =', oute.shape) + # oute should be: + # array([[-2.27165658, -3.01387112, -3.98753036, -3.92317787, -2.43752966], + # [ 1.09062263, 0.34659524, -0.63331812, -0.47944498, 1.01023176]]) + print('oute =', oute) + print('') + + idd = imas.ids(shot, 0) + idd.open_env('imas_public', 'west', '3') + idd.equilibrium.get() + out = idd.equilibrium + + equiDict = {} + + # Declaration of arrays 2d plots + equi_grid = idd.equilibrium.grids_ggd[0].grid[0] + NbrPoints = len(equi_grid.space[0].objects_per_dimension[0].object) + equiDict['r'] = np.full(NbrPoints, np.nan) + equiDict['z'] = np.full(NbrPoints, np.nan) + for ii in range(NbrPoints): + equiDict['r'][ii] = equi_grid.space[0].objects_per_dimension[0]. \ + object[ii].geometry[0] + equiDict['z'][ii] = equi_grid.space[0].objects_per_dimension[0]. \ + object[ii].geometry[1] + + arg_time = np.argmin(np.abs(out.time - approx_time_in)) + + print('\nDistance between requested and calculated time =', \ + np.abs(out.time[arg_time] - approx_time_in), '\n') + + time_in = out.time[arg_time] - eps_time + + print('\nNew distance between requested and calculated time =', \ + np.abs(out.time[arg_time] - time_in), '\n', \ + 'with eps_time =', eps_time, '\n') + + R_all = np.linspace(np.min(equiDict['r']), np.max(equiDict['r']), interp_points) + Z_all = np.linspace(np.min(equiDict['z']), np.max(equiDict['z']), interp_points) + + R_all_tot = np.repeat(R_all, interp_points) + Z_all_tot = np.tile(Z_all, interp_points) + + # CALL EQUIMAP + outa = equimap.get(shot, time=time_in, \ + R=R_all_tot, Phi=np.zeros(R_all_tot.shape), Z=Z_all_tot, \ + quantity='psi') + + outar = outa.reshape((interp_points, interp_points)) + Rr = R_all_tot.reshape((interp_points, interp_points)) + Zr = Z_all_tot.reshape((interp_points, interp_points)) + + print('') + print('outa.shape =', outa.shape) + print('outar.shape =', outar.shape) + #print('oute =', oute) + print('') + + #plt.figure() + #cs_ref = plt.contour(out.interp2D.r, out.interp2D.z, \ + # out.interp2D.psi[np.argmin(np.abs(out.time - time_in))], \ + # ctr_plot, colors='k', linestyle='solid') + #cs_ref_tri = plt.tricontour(equiDict['r'], equiDict['z'], \ + # out.ggd[0].grid.space[0].objects_per_dimension[1].nodes, \ + # out.ggd[0].psi[np.argmin(np.abs(out.time - time_in))], \ + # ctr_plot, colors='red', linestyles='dashdot') + #cs = plt.contour(Rr, Zr, outar, ctr_plot, \ + # colors='orange', linestyles='dashed') + #plt.plot(np.nan, np.nan, color='k', linestyle='solid', label='ref interp2D') + #plt.plot(np.nan, np.nan, color='red', linestyle='dashdot', label='ref ggd') + #plt.plot(np.nan, np.nan, color='orange', linestyle='dashed', label='test interp') + #plt.legend() + #plt.xlabel('R [m]') + #plt.ylabel('Z [m]') + #plt.title('Contour psi, level ' + str(ctr_plot) + ' at time ' + str(time_in) + 's') + #plt.axes().set(aspect=1) + # + #p = cs.collections[0].get_paths()[0] + #v = p.vertices + #x_cs = v[:, 0] + #y_cs = v[:, 1] + + #p = cs_ref.collections[0].get_paths()[0] + #v = p.vertices + #x_cs_ref = v[:, 0] + #y_cs_ref = v[:, 1] + + plt.figure() + plt.contourf(Rr, Zr, outar) + plt.plot(np.squeeze(out.time_slice[arg_time].boundary.outline.r), \ + np.squeeze(out.time_slice[arg_time].boundary.outline.z), \ + linewidth=2, color='red') + plt.plot(out.time_slice[arg_time].global_quantities.magnetic_axis.r, \ + out.time_slice[arg_time].global_quantities.magnetic_axis.z, \ + marker='+', color='red', markersize=20) + plt.colorbar() + plt.xlabel('R [m]') + plt.ylabel('Z [m]') + plt.title('psi interpolated at time ' + str(time_in) + 's') + plt.axes().set(aspect=1) + + #plt.figure() + #plt.plot(x_cs_ref, y_cs_ref, label='Ref') + #plt.plot(x_cs, y_cs, '--', label='Interp') + #plt.xlabel('R [m]') + #plt.ylabel('Z [m]') + #plt.legend() + + #pts_cs_ref = np.linspace(0, 1, x_cs_ref.shape[0]) + #pts_cs = np.linspace(0, 1, x_cs.shape[0]) + + #x_cs_intp = np.interp(pts_cs_ref, pts_cs, x_cs) + #y_cs_intp = np.interp(pts_cs_ref, pts_cs, y_cs) + + #error_x = np.abs((x_cs_intp - x_cs_ref)/x_cs_ref)*100 + #error_y = np.abs((y_cs_intp - y_cs_ref)/y_cs_ref)*100 + + #plt.figure() + #plt.plot(error_x, label='err. x') + #plt.plot(error_y, label='err. y') + #plt.legend() + #plt.ylabel('Error interp psi contour level ' + str(ctr_plot) + ' in [%]' \ + # + '\nwith time difference, eps_time = ' + str(eps_time*100) + '%', \ + # multialignment='center') + + # CALL EQUIMAP + out_rho_pol_norm = equimap.get(shot, time=time_in, \ + R=R_all_tot, Phi=np.zeros(R_all_tot.shape), Z=Z_all_tot, \ + quantity='rho_pol_norm') + + print('') + print('out_rho_pol_norm.shape =', out_rho_pol_norm.shape) + #print('out_rho_pol_norm =', out_rho_pol_norm) + print('') + + plt.figure() + out_rho_pol_norm_r = out_rho_pol_norm.reshape((interp_points, interp_points)) + cs = plt.contourf(Rr, Zr, out_rho_pol_norm_r) + plt.colorbar() + plt.plot(np.squeeze(out.time_slice[arg_time].boundary.outline.r), \ + np.squeeze(out.time_slice[arg_time].boundary.outline.z), \ + linewidth=2, color='red') + plt.plot(out.time_slice[arg_time].global_quantities.magnetic_axis.r, \ + out.time_slice[arg_time].global_quantities.magnetic_axis.z, \ + marker='+', color='red', markersize=20) + plt.xlabel('R [m]') + plt.ylabel('Z [m]') + plt.title('rho_pol_norm') + + # CALL EQUIMAP + out_rho_tor_norm = equimap.get(shot, time=time_in, \ + R=R_all_tot, Phi=np.zeros(R_all_tot.shape), Z=Z_all_tot, \ + quantity='rho_tor_norm') + print('') + print('out_rho_tor_norm.shape =', out_rho_tor_norm.shape) + #print('out_rho_tor_norm =', out_rho_tor_norm) + print('') + plt.figure() + out_rho_tor_norm_r = out_rho_tor_norm.reshape((interp_points, interp_points)) + cs = plt.contourf(Rr, Zr, out_rho_tor_norm_r) + plt.colorbar() + arg_time = np.argmin(np.abs(out.time - time_in)) + plt.plot(np.squeeze(out.time_slice[arg_time].boundary.outline.r), \ + np.squeeze(out.time_slice[arg_time].boundary.outline.z), \ + linewidth=2, color='red') + plt.plot(out.time_slice[arg_time].global_quantities.magnetic_axis.r, \ + out.time_slice[arg_time].global_quantities.magnetic_axis.z, \ + marker='+', color='red', markersize=20) + plt.xlabel('R [m]') + plt.ylabel('Z [m]') + plt.title('rho_tor_norm') + + # CALL EQUIMAP + out_rho_tor = equimap.get(shot, time=time_in, \ + R=R_all_tot, Phi=np.zeros(R_all_tot.shape), Z=Z_all_tot, \ + quantity='rho_tor') + + print('') + print('out_rho_tor.shape =', out_rho_tor.shape) + #print('out_rho_tor =', out_rho_tor) + print('') + plt.figure() + out_rho_tor_r = out_rho_tor.reshape((interp_points, interp_points)) + cs = plt.contourf(Rr, Zr, out_rho_tor_r) + plt.colorbar() + arg_time = np.argmin(np.abs(out.time - time_in)) + plt.plot(np.squeeze(out.time_slice[arg_time].boundary.outline.r), \ + np.squeeze(out.time_slice[arg_time].boundary.outline.z), \ + linewidth=2, color='red') + plt.plot(out.time_slice[arg_time].global_quantities.magnetic_axis.r, \ + out.time_slice[arg_time].global_quantities.magnetic_axis.z, \ + marker='+', color='red', markersize=20) + plt.xlabel('R [m]') + plt.ylabel('Z [m]') + plt.title('rho_tor') + + plt.figure() + plt.plot(out.time_slice[arg_time].profiles_1d.psi, out.time_slice[arg_time].profiles_1d.rho_tor) + plt.xlabel('Psi') + plt.ylabel('rho_tor') + + #plt.figure() + #plt.tricontourf(equiDict['r'], equiDict['z'], \ + # out.ggd[0].grid.space[0].objects_per_dimension[1].nodes, \ + # out.ggd[0].psi[arg_time]) + #plt.plot(np.squeeze(out.boundary.outline.r[arg_time]), \ + # np.squeeze(out.boundary.outline.z[arg_time]), \ + # linewidth=2, color='red') + #plt.plot(out.global_quantities.magnetic_axis.r[arg_time], \ + # out.global_quantities.magnetic_axis.z[arg_time], \ + # marker='+', color='red', markersize=20) + #plt.colorbar() + #plt.xlabel('R [m]') + #plt.ylabel('Z [m]') + #plt.title('psi ggd') + + # CALL EQUIMAP + out_b_norm = equimap.get(shot, time=time_in, \ + R=R_all_tot, Phi=np.zeros(R_all_tot.shape), Z=Z_all_tot, \ + quantity='b_field_norm') + + print('') + print('out_b_norm.shape =', out_b_norm.shape) + #print('out_b_norm =', out_b_norm) + print('') + plt.figure() + out_b_norm_r = out_b_norm.reshape((interp_points, interp_points)) + cs = plt.contourf(Rr, Zr, out_b_norm_r) + plt.colorbar() + arg_time = np.argmin(np.abs(out.time - time_in)) + plt.plot(np.squeeze(out.time_slice[arg_time].boundary.outline.r), \ + np.squeeze(out.time_slice[arg_time].boundary.outline.z), \ + linewidth=2, color='red') + plt.plot(out.time_slice[arg_time].global_quantities.magnetic_axis.r, \ + out.time_slice[arg_time].global_quantities.magnetic_axis.z, \ + marker='+', color='red', markersize=20) + plt.xlabel('R [m]') + plt.ylabel('Z [m]') + plt.title('b_norm') + + # CALL EQUIMAP + out_theta = equimap.get(shot, time=time_in, \ + R=R_all_tot, Phi=np.zeros(R_all_tot.shape), Z=Z_all_tot, \ + quantity='theta') + + print('') + print('out_theta.shape =', out_theta.shape) + print('out_theta =', out_theta) + print('') + plt.figure() + out_theta_r = out_theta.reshape((interp_points, interp_points)) + cs = plt.contourf(Rr, Zr, out_theta_r, 100) + plt.colorbar() + arg_time = np.argmin(np.abs(out.time - time_in)) + plt.plot(np.squeeze(out.time_slice[arg_time].boundary.outline.r), \ + np.squeeze(out.time_slice[arg_time].boundary.outline.z), \ + linewidth=2, color='red') + plt.plot(out.time_slice[arg_time].global_quantities.magnetic_axis.r, \ + out.time_slice[arg_time].global_quantities.magnetic_axis.z, \ + marker='+', color='red', markersize=20) + plt.xlabel('R [m]') + plt.ylabel('Z [m]') + plt.title('theta') + + plt.show() diff --git a/tofu/mag/test_magFieldLines.py b/tofu/mag/test_magFieldLines.py new file mode 100644 index 000000000..d1a44b66b --- /dev/null +++ b/tofu/mag/test_magFieldLines.py @@ -0,0 +1,20 @@ +# -*- coding: utf-8 -*- +# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 +# Also if needed: retab +''' + TEST magFieldLines +''' +from __future__ import (unicode_literals, absolute_import, \ + print_function, division) +import doctest +import os +import sys + +if __name__ == '__main__': + #print('path 1 =', sys.path) + sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), '..'))) + #print('path 2 =', sys.path) + #import mag + from mag import magFieldLines + path = sys.path.pop(0) + doctest.testmod(magFieldLines, verbose=True) diff --git a/tofu/mag/test_ripple.py b/tofu/mag/test_ripple.py index 8006aeded..1880dd861 100644 --- a/tofu/mag/test_ripple.py +++ b/tofu/mag/test_ripple.py @@ -13,435 +13,436 @@ import time #import warnings -#print('path 1 =', sys.path) -sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), '..'))) -#print('path 2 =', sys.path) - -# Local modules -import imas -import equimap -#import imas_west -#import pywed as pw - -shot = 53221 -tol_val = 1E-10 - -# For 2D plots -interp_points = 60 - -# FIRST POINT B_NORM -# ------------------ -time_in = np.linspace(36, 37, 10) -Phi_in = np.linspace(0, 2*np.pi/18, 100) -R_in = np.full(Phi_in.shape, 3) -Z_in = np.zeros(R_in.shape) - -# Read equilibrium data -idd = imas.ids(shot, 0) -idd.open_env('imas_public', 'west', '3') -idd.equilibrium.get() -out = idd.equilibrium - -equiDict = {} - -# Declaration of arrays 2d plots -equi_grid = idd.equilibrium.grids_ggd[0].grid[0] -NbrPoints = len(equi_grid.space[0].objects_per_dimension[0].object) -equiDict['r'] = np.full(NbrPoints, np.nan) -equiDict['z'] = np.full(NbrPoints, np.nan) -for ii in range(NbrPoints): - equiDict['r'][ii] = equi_grid.space[0].objects_per_dimension[0]. \ - object[ii].geometry[0] - equiDict['z'][ii] = equi_grid.space[0].objects_per_dimension[0]. \ - object[ii].geometry[1] - -# For 2D plots -R_all = np.linspace(np.min(equiDict['r']), np.max(equiDict['r']), interp_points) -Z_all = np.linspace(np.min(equiDict['z']), np.max(equiDict['z']), interp_points) - -R_all_tot = np.repeat(R_all, interp_points) -Z_all_tot = np.tile(Z_all, interp_points) - -Rr = R_all_tot.reshape((interp_points, interp_points)) -Zr = Z_all_tot.reshape((interp_points, interp_points)) - -# CALL EQUIMAP -start = time.time() -oute = equimap.get(shot, time=time_in, \ - R=R_in, Phi=Phi_in, Z=Z_in, \ - quantity='b_field_norm') -end = time.time() -print() -print('time in equimap.get b_norm =', end - start) -print() -print('oute.shape b_norm =', oute.shape) - -# CALL EQUIMAP -start = time.time() -oute_noR = equimap.get(shot, time=time_in, \ - R=R_in, Phi=Phi_in, Z=Z_in, \ - quantity='b_field_norm', no_ripple=True) -end = time.time() -print() -print('time in equimap.get b_norm no Ripple =', end - start) -print() -print('oute.shape b_norm no ripple =', oute_noR.shape) - -print() -print('Mean value B_norm ripple =', np.mean(oute[int(0.5*oute.shape[0]), :])) -print('Mean value B_norm NO ripple =', \ - np.mean(oute_noR[int(0.5*oute_noR.shape[0]), :])) -diff_mean_val = np.mean(oute[int(0.5*oute.shape[0]), :]) \ - - np.mean(oute_noR[int(0.5*oute_noR.shape[0]), :]) -print('Diff mean values =', diff_mean_val) -percent_diff = np.abs(100*diff_mean_val \ - / np.mean(oute[int(0.5*oute.shape[0]), :])) -print('Percent diff mean values =', percent_diff) - -# CHECK -# ----- -if (np.abs(percent_diff - 0.011052598088) > tol_val): - print() - print('ERROR: Higher than tolerance percent difference ' \ - + str(np.abs(percent_diff - 0.011052598088))) - print() - #raise RuntimeError -# FOR: -# shot = 53221 -# time_in = np.linspace(36, 37, 10) -# Phi_in = np.linspace(0, 2*np.pi/18, 100) -# R_in = np.full(Phi_in.shape, 3) -# Z_in = np.zeros(R_in.shape) -# RESULTS: -# Mean value B_norm ripple = 3.05593472975 -# Mean value B_norm NO ripple = 3.05627248994 -# Diff mean values = -0.000337760183512 -# Percent diff mean values = 0.011052598088 -print() - -# PLOTS -plt.figure() -plt.plot(time_in, oute[:, -1], label='B_norm at R={0}, Phi=Z=0'.format(R_in[-1])) -plt.plot(time_in, oute_noR[:, -1], label='B_norm no ripple at R={0}, Phi=Z=0'.format(R_in[-1])) -plt.legend() -plt.xlabel('Time [s]') -plt.ylabel('B_norm [T]') -plt.figure() -plt.plot(Phi_in, oute[int(0.5*oute.shape[0]), :], \ - label='B_norm at t={0:.2f}, R={1}, Z=0'.format( \ - time_in[int(0.5*oute.shape[0])], R_in[-1])) -plt.plot(Phi_in, oute_noR[int(0.5*oute.shape[0]), :], \ - label='B_norm no ripple at t={0:.2f}, R={1}, Z=0'.format( \ - time_in[int(0.5*oute.shape[0])], R_in[-1])) -plt.legend() -plt.xlabel('Phi [rad]') -plt.ylabel('B_norm [T]') - -# SECOND POSITION B_NORM -# ---------------------- -t_ignitron = [] -t_ignitron.append(32) -print() -print('t_igni =', t_ignitron[0]) -print() -time_in = np.linspace(t_ignitron[0], 38, 10) -Phi_in = np.linspace(0, 2*np.pi/18, 100) -R_in = np.full(Phi_in.shape, 2.43) -Z_in = np.full(Phi_in.shape, 0.57) - -# CALL EQUIMAP -start = time.time() -oute = equimap.get(shot, time=time_in, \ - R=R_in, Phi=Phi_in, Z=Z_in, \ - quantity='b_field_norm') -end = time.time() -print() -print('time in equimap.get 2 b_norm =', end - start) -print() -print('oute.shape 2 b_norm =', oute.shape) - -# CALL EQUIMAP -start = time.time() -oute_noR = equimap.get(shot, time=time_in, \ - R=R_in, Phi=Phi_in, Z=Z_in, \ - quantity='b_field_norm', no_ripple=True) -end = time.time() -print() -print('time in equimap.get 2 b_norm no ripple =', end - start) -print() -print('oute.shape 2 b_norm no ripple =', oute_noR.shape) - -# PLOTS -plt.figure() -plt.plot(time_in, oute[:, -1], \ - label='B_norm at R={0}, Phi={1:.2f}, Z={2}'.format( \ - R_in[-1], Phi_in[-1], Z_in[-1])) -plt.plot(time_in, oute_noR[:, -1], \ - label='B_norm no ripple at R={0}, Phi={1:.2f}, Z={2}'.format( \ - R_in[-1], Phi_in[-1], Z_in[-1])) -plt.legend() -plt.xlabel('Time [s]') -plt.ylabel('B_norm [T]') -plt.figure() -plt.plot(Phi_in, oute[int(0.5*oute.shape[0]), :], \ - label='B_norm at t={0:.2f}, R={1}, Z={2}'.format( \ - time_in[int(0.5*oute.shape[0])], R_in[-1], Z_in[-1])) -plt.plot(Phi_in, oute_noR[int(0.5*oute.shape[0]), :], \ - label='B_norm no ripple at t={0:.2f}, R={1}, Z={2}'.format( \ - time_in[int(0.5*oute.shape[0])], R_in[-1], Z_in[-1])) -plt.legend() -plt.xlabel('Phi [rad]') -plt.ylabel('B_norm [T]') - -# B_NORM 2D -# --------- -# CALL EQUIMAP -start = time.time() -outa = equimap.get(shot, time=time_in, \ - R=R_all_tot, Phi=np.zeros(R_all_tot.shape), Z=Z_all_tot, \ - quantity='b_field_norm') -end = time.time() -print() -print('time in equimap.get b_norm 2D =', end - start) -print() -outar = outa[int(0.5*outa.shape[0])].reshape((interp_points, interp_points)) -plt.figure() -plt.contourf(Rr, Zr, outar) -plt.colorbar() -arg_time = np.argmin(np.abs(out.time - time_in[int(0.5*outa.shape[0])])) -plt.plot(np.squeeze(out.time_slice[arg_time].boundary.outline.r), \ - np.squeeze(out.time_slice[arg_time].boundary.outline.z), \ - linewidth=2, color='red') -plt.plot(out.time_slice[arg_time].global_quantities.magnetic_axis.r, \ - out.time_slice[arg_time].global_quantities.magnetic_axis.z, \ - marker='+', color='red', markersize=20) -plt.xlabel('R [m]') -plt.ylabel('Z [m]') -plt.title('B_norm t={0:.2f}'.format(time_in[int(0.5*outa.shape[0])])) - -# B_R TEST -# -------- -Phi_in = np.linspace(0, 2*np.pi/18, 100) -R_in = np.full(Phi_in.shape, 3) -Z_in = np.full(Phi_in.shape, 0) - -# CALL EQUIMAP -start = time.time() -oute = equimap.get(shot, time=time_in, \ - R=R_in, Phi=Phi_in, Z=Z_in, \ - quantity='b_field_r') -end = time.time() -print() -print('time in equimap.get br =', end - start) -print() -print('oute.shape br =', oute.shape) - -# CALL EQUIMAP -start = time.time() -oute_noR = equimap.get(shot, time=time_in, \ - R=R_in, Phi=Phi_in, Z=Z_in, \ - quantity='b_field_r', no_ripple=True) -end = time.time() -print() -print('time in equimap.get br no ripple =', end - start) -print() -print('oute.shape br no ripple =', oute_noR.shape) - -# PLOTS -plt.figure() -plt.plot(time_in, oute[:, -1], \ - label='B_r at R={0}, Phi={1:.2f}, Z={2}'.format( \ - R_in[-1], Phi_in[-1], Z_in[-1])) -plt.plot(time_in, oute_noR[:, -1], \ - label='B_r no ripple at R={0}, Phi={1:.2f}, Z={2}'.format( \ - R_in[-1], Phi_in[-1], Z_in[-1])) -plt.legend() -plt.xlabel('Time [s]') -plt.ylabel('B_r [T]') - -plt.figure() -plt.plot(Phi_in, oute[int(0.5*oute.shape[0]), :], \ - label='B_r at t={0:.2f}, R={1}, Z={2}'.format( \ - time_in[int(0.5*oute.shape[0])], R_in[-1], Z_in[-1])) -plt.plot(Phi_in, oute_noR[int(0.5*oute.shape[0]), :], \ - label='B_r no ripple at t={0:.2f}, R={1}, Z={2}'.format( \ - time_in[int(0.5*oute.shape[0])], R_in[-1], Z_in[-1])) -plt.legend() -plt.xlabel('Phi [rad]') -plt.ylabel('B_r [T]') - -# CALL EQUIMAP -start = time.time() -outa = equimap.get(shot, time=time_in, \ - R=R_all_tot, Phi=np.zeros(R_all_tot.shape), Z=Z_all_tot, \ - quantity='b_field_r') -end = time.time() -print() -print('time in equimap.get br 2D =', end - start) -print() -outar = outa[int(0.5*outa.shape[0])].reshape((interp_points, interp_points)) -plt.figure() -plt.contourf(Rr, Zr, outar) -plt.colorbar() -arg_time = np.argmin(np.abs(out.time - time_in[int(0.5*outa.shape[0])])) -plt.plot(np.squeeze(out.time_slice[arg_time].boundary.outline.r), \ - np.squeeze(out.time_slice[arg_time].boundary.outline.z), \ - linewidth=2, color='red') -plt.plot(out.time_slice[arg_time].global_quantities.magnetic_axis.r, \ - out.time_slice[arg_time].global_quantities.magnetic_axis.z, \ - marker='+', color='red', markersize=20) -plt.xlabel('R [m]') -plt.ylabel('Z [m]') -plt.title('B_r t={0:.2f}'.format(time_in[int(0.5*outa.shape[0])])) - -# B_Z TEST -# -------- -Phi_in = np.linspace(0, 2*np.pi/18, 100) -R_in = np.full(Phi_in.shape, 3) -Z_in = np.full(Phi_in.shape, 0.2) - -# CALL EQUIMAP -start = time.time() -oute = equimap.get(shot, time=time_in, \ - R=R_in, Phi=Phi_in, Z=Z_in, \ - quantity='b_field_z') -end = time.time() -print() -print('time in equimap.get bz =', end - start) -print() -print('oute.shape bz =', oute.shape) - -# CALL EQUIMAP -start = time.time() -oute_noR = equimap.get(shot, time=time_in, \ - R=R_in, Phi=Phi_in, Z=Z_in, \ - quantity='b_field_z', no_ripple=True) -end = time.time() -print() -print('time in equimap.get bz no ripple =', end - start) -print() -print('oute.shape bz no ripple =', oute_noR.shape) - -# PLOTS -plt.figure() -plt.plot(time_in, oute[:, -1], \ - label='B_z at R={0}, Phi={1:.2f}, Z={2}'.format( \ - R_in[-1], Phi_in[-1], Z_in[-1])) -plt.plot(time_in, oute_noR[:, -1], \ - label='B_z no ripple at R={0}, Phi={1:.2f}, Z={2}'.format( \ - R_in[-1], Phi_in[-1], Z_in[-1])) -plt.legend() -plt.xlabel('Time [s]') -plt.ylabel('B_z [T]') - -plt.figure() -plt.plot(Phi_in, oute[int(0.5*oute.shape[0]), :], \ - label='B_z at t={0:.2f}, R={1}, Z={2}'.format( \ - time_in[int(0.5*oute.shape[0])], R_in[-1], Z_in[-1])) -plt.plot(Phi_in, oute_noR[int(0.5*oute.shape[0]), :], \ - label='B_z no ripple at t={0:.2f}, R={1}, Z={2}'.format( \ - time_in[int(0.5*oute.shape[0])], R_in[-1], Z_in[-1])) -plt.legend() -plt.xlabel('Phi [rad]') -plt.ylabel('B_z [T]') - -# CALL EQUIMAP -start = time.time() -outa = equimap.get(shot, time=time_in, \ - R=R_all_tot, Phi=np.zeros(R_all_tot.shape), Z=Z_all_tot, \ - quantity='b_field_z') -end = time.time() -print() -print('time in equimap.get bz 2D =', end - start) -print() -outar = outa[int(0.5*outa.shape[0])].reshape((interp_points, interp_points)) -plt.figure() -plt.contourf(Rr, Zr, outar) -plt.colorbar() -arg_time = np.argmin(np.abs(out.time - time_in[int(0.5*outa.shape[0])])) -plt.plot(np.squeeze(out.time_slice[arg_time].boundary.outline.r), \ - np.squeeze(out.time_slice[arg_time].boundary.outline.z), \ - linewidth=2, color='red') -plt.plot(out.time_slice[arg_time].global_quantities.magnetic_axis.r, \ - out.time_slice[arg_time].global_quantities.magnetic_axis.z, \ - marker='+', color='red', markersize=20) -plt.xlabel('R [m]') -plt.ylabel('Z [m]') -plt.title('B_z t={0:.2f}'.format(time_in[int(0.5*outa.shape[0])])) - -# B_TOR TEST -# ---------- -Phi_in = np.linspace(0, 2*np.pi/18, 100) -R_in = np.full(Phi_in.shape, 3) -Z_in = np.full(Phi_in.shape, 0) - -# CALL EQUIMAP -start = time.time() -oute = equimap.get(shot, time=time_in, \ - R=R_in, Phi=Phi_in, Z=Z_in, \ - quantity='b_field_tor') -end = time.time() -print() -print('time in equimap.get btor =', end - start) -print() -print('oute.shape btor =', oute.shape) - -# CALL EQUIMAP -start = time.time() -oute_noR = equimap.get(shot, time=time_in, \ - R=R_in, Phi=Phi_in, Z=Z_in, \ - quantity='b_field_tor', no_ripple=True) -end = time.time() -print() -print('time in equimap.get btor no ripple =', end - start) -print() -print('oute.shape btor no ripple =', oute_noR.shape) - -# PLOTS -plt.figure() -plt.plot(time_in, oute[:, -1], \ - label='B_tor at R={0}, Phi={1:.2f}, Z={2}'.format( \ - R_in[-1], Phi_in[-1], Z_in[-1])) -plt.plot(time_in, oute_noR[:, -1], \ - label='B_tor no ripple at R={0}, Phi={1:.2f}, Z={2}'.format( \ - R_in[-1], Phi_in[-1], Z_in[-1])) -plt.legend() -plt.xlabel('Time [s]') -plt.ylabel('B_tor [T]') - -plt.figure() -plt.plot(Phi_in, oute[int(0.5*oute.shape[0]), :], \ - label='B_tor at t={0:.2f}, R={1}, Z={2}'.format( \ - time_in[int(0.5*oute.shape[0])], R_in[-1], Z_in[-1])) -plt.plot(Phi_in, oute_noR[int(0.5*oute.shape[0]), :], \ - label='B_tor no ripple at t={0:.2f}, R={1}, Z={2}'.format( \ - time_in[int(0.5*oute.shape[0])], R_in[-1], Z_in[-1])) -plt.legend() -plt.xlabel('Phi [rad]') -plt.ylabel('B_tor [T]') - -# CALL EQUIMAP -start = time.time() -outa = equimap.get(shot, time=time_in, \ - R=R_all_tot, Phi=np.zeros(R_all_tot.shape), Z=Z_all_tot, \ - quantity='b_field_tor') -end = time.time() -print() -print('time in equimap.get btor 2D =', end - start) -print() -outar = outa[int(0.5*outa.shape[0])].reshape((interp_points, interp_points)) -plt.figure() -plt.contourf(Rr, Zr, outar) -plt.colorbar() -arg_time = np.argmin(np.abs(out.time - time_in[int(0.5*outa.shape[0])])) -plt.plot(np.squeeze(out.time_slice[arg_time].boundary.outline.r), \ - np.squeeze(out.time_slice[arg_time].boundary.outline.z), \ - linewidth=2, color='red') -plt.plot(out.time_slice[arg_time].global_quantities.magnetic_axis.r, \ - out.time_slice[arg_time].global_quantities.magnetic_axis.z, \ - marker='+', color='red', markersize=20) -plt.xlabel('R [m]') -plt.ylabel('Z [m]') -plt.title('B_tor t={0:.2f}'.format(time_in[int(0.5*outa.shape[0])])) - -plt.show() +if __name__ == '__main__': + #print('path 1 =', sys.path) + sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), '..'))) + #print('path 2 =', sys.path) + + # Local modules + import imas + import equimap + #import imas_west + #import pywed as pw + + shot = 53221 + tol_val = 1E-10 + + # For 2D plots + interp_points = 60 + + # FIRST POINT B_NORM + # ------------------ + time_in = np.linspace(36, 37, 10) + Phi_in = np.linspace(0, 2*np.pi/18, 100) + R_in = np.full(Phi_in.shape, 3) + Z_in = np.zeros(R_in.shape) + + # Read equilibrium data + idd = imas.ids(shot, 0) + idd.open_env('imas_public', 'west', '3') + idd.equilibrium.get() + out = idd.equilibrium + + equiDict = {} + + # Declaration of arrays 2d plots + equi_grid = idd.equilibrium.grids_ggd[0].grid[0] + NbrPoints = len(equi_grid.space[0].objects_per_dimension[0].object) + equiDict['r'] = np.full(NbrPoints, np.nan) + equiDict['z'] = np.full(NbrPoints, np.nan) + for ii in range(NbrPoints): + equiDict['r'][ii] = equi_grid.space[0].objects_per_dimension[0]. \ + object[ii].geometry[0] + equiDict['z'][ii] = equi_grid.space[0].objects_per_dimension[0]. \ + object[ii].geometry[1] + + # For 2D plots + R_all = np.linspace(np.min(equiDict['r']), np.max(equiDict['r']), interp_points) + Z_all = np.linspace(np.min(equiDict['z']), np.max(equiDict['z']), interp_points) + + R_all_tot = np.repeat(R_all, interp_points) + Z_all_tot = np.tile(Z_all, interp_points) + + Rr = R_all_tot.reshape((interp_points, interp_points)) + Zr = Z_all_tot.reshape((interp_points, interp_points)) + + # CALL EQUIMAP + start = time.time() + oute = equimap.get(shot, time=time_in, \ + R=R_in, Phi=Phi_in, Z=Z_in, \ + quantity='b_field_norm') + end = time.time() + print() + print('time in equimap.get b_norm =', end - start) + print() + print('oute.shape b_norm =', oute.shape) + + # CALL EQUIMAP + start = time.time() + oute_noR = equimap.get(shot, time=time_in, \ + R=R_in, Phi=Phi_in, Z=Z_in, \ + quantity='b_field_norm', no_ripple=True) + end = time.time() + print() + print('time in equimap.get b_norm no Ripple =', end - start) + print() + print('oute.shape b_norm no ripple =', oute_noR.shape) + + print() + print('Mean value B_norm ripple =', np.mean(oute[int(0.5*oute.shape[0]), :])) + print('Mean value B_norm NO ripple =', \ + np.mean(oute_noR[int(0.5*oute_noR.shape[0]), :])) + diff_mean_val = np.mean(oute[int(0.5*oute.shape[0]), :]) \ + - np.mean(oute_noR[int(0.5*oute_noR.shape[0]), :]) + print('Diff mean values =', diff_mean_val) + percent_diff = np.abs(100*diff_mean_val \ + / np.mean(oute[int(0.5*oute.shape[0]), :])) + print('Percent diff mean values =', percent_diff) + + # CHECK + # ----- + if (np.abs(percent_diff - 0.011052598088) > tol_val): + print() + print('ERROR: Higher than tolerance percent difference ' \ + + str(np.abs(percent_diff - 0.011052598088))) + print() + #raise RuntimeError + # FOR: + # shot = 53221 + # time_in = np.linspace(36, 37, 10) + # Phi_in = np.linspace(0, 2*np.pi/18, 100) + # R_in = np.full(Phi_in.shape, 3) + # Z_in = np.zeros(R_in.shape) + # RESULTS: + # Mean value B_norm ripple = 3.05593472975 + # Mean value B_norm NO ripple = 3.05627248994 + # Diff mean values = -0.000337760183512 + # Percent diff mean values = 0.011052598088 + print() + + # PLOTS + plt.figure() + plt.plot(time_in, oute[:, -1], label='B_norm at R={0}, Phi=Z=0'.format(R_in[-1])) + plt.plot(time_in, oute_noR[:, -1], label='B_norm no ripple at R={0}, Phi=Z=0'.format(R_in[-1])) + plt.legend() + plt.xlabel('Time [s]') + plt.ylabel('B_norm [T]') + plt.figure() + plt.plot(Phi_in, oute[int(0.5*oute.shape[0]), :], \ + label='B_norm at t={0:.2f}, R={1}, Z=0'.format( \ + time_in[int(0.5*oute.shape[0])], R_in[-1])) + plt.plot(Phi_in, oute_noR[int(0.5*oute.shape[0]), :], \ + label='B_norm no ripple at t={0:.2f}, R={1}, Z=0'.format( \ + time_in[int(0.5*oute.shape[0])], R_in[-1])) + plt.legend() + plt.xlabel('Phi [rad]') + plt.ylabel('B_norm [T]') + + # SECOND POSITION B_NORM + # ---------------------- + t_ignitron = [] + t_ignitron.append(32) + print() + print('t_igni =', t_ignitron[0]) + print() + time_in = np.linspace(t_ignitron[0], 38, 10) + Phi_in = np.linspace(0, 2*np.pi/18, 100) + R_in = np.full(Phi_in.shape, 2.43) + Z_in = np.full(Phi_in.shape, 0.57) + + # CALL EQUIMAP + start = time.time() + oute = equimap.get(shot, time=time_in, \ + R=R_in, Phi=Phi_in, Z=Z_in, \ + quantity='b_field_norm') + end = time.time() + print() + print('time in equimap.get 2 b_norm =', end - start) + print() + print('oute.shape 2 b_norm =', oute.shape) + + # CALL EQUIMAP + start = time.time() + oute_noR = equimap.get(shot, time=time_in, \ + R=R_in, Phi=Phi_in, Z=Z_in, \ + quantity='b_field_norm', no_ripple=True) + end = time.time() + print() + print('time in equimap.get 2 b_norm no ripple =', end - start) + print() + print('oute.shape 2 b_norm no ripple =', oute_noR.shape) + + # PLOTS + plt.figure() + plt.plot(time_in, oute[:, -1], \ + label='B_norm at R={0}, Phi={1:.2f}, Z={2}'.format( \ + R_in[-1], Phi_in[-1], Z_in[-1])) + plt.plot(time_in, oute_noR[:, -1], \ + label='B_norm no ripple at R={0}, Phi={1:.2f}, Z={2}'.format( \ + R_in[-1], Phi_in[-1], Z_in[-1])) + plt.legend() + plt.xlabel('Time [s]') + plt.ylabel('B_norm [T]') + plt.figure() + plt.plot(Phi_in, oute[int(0.5*oute.shape[0]), :], \ + label='B_norm at t={0:.2f}, R={1}, Z={2}'.format( \ + time_in[int(0.5*oute.shape[0])], R_in[-1], Z_in[-1])) + plt.plot(Phi_in, oute_noR[int(0.5*oute.shape[0]), :], \ + label='B_norm no ripple at t={0:.2f}, R={1}, Z={2}'.format( \ + time_in[int(0.5*oute.shape[0])], R_in[-1], Z_in[-1])) + plt.legend() + plt.xlabel('Phi [rad]') + plt.ylabel('B_norm [T]') + + # B_NORM 2D + # --------- + # CALL EQUIMAP + start = time.time() + outa = equimap.get(shot, time=time_in, \ + R=R_all_tot, Phi=np.zeros(R_all_tot.shape), Z=Z_all_tot, \ + quantity='b_field_norm') + end = time.time() + print() + print('time in equimap.get b_norm 2D =', end - start) + print() + outar = outa[int(0.5*outa.shape[0])].reshape((interp_points, interp_points)) + plt.figure() + plt.contourf(Rr, Zr, outar) + plt.colorbar() + arg_time = np.argmin(np.abs(out.time - time_in[int(0.5*outa.shape[0])])) + plt.plot(np.squeeze(out.time_slice[arg_time].boundary.outline.r), \ + np.squeeze(out.time_slice[arg_time].boundary.outline.z), \ + linewidth=2, color='red') + plt.plot(out.time_slice[arg_time].global_quantities.magnetic_axis.r, \ + out.time_slice[arg_time].global_quantities.magnetic_axis.z, \ + marker='+', color='red', markersize=20) + plt.xlabel('R [m]') + plt.ylabel('Z [m]') + plt.title('B_norm t={0:.2f}'.format(time_in[int(0.5*outa.shape[0])])) + + # B_R TEST + # -------- + Phi_in = np.linspace(0, 2*np.pi/18, 100) + R_in = np.full(Phi_in.shape, 3) + Z_in = np.full(Phi_in.shape, 0) + + # CALL EQUIMAP + start = time.time() + oute = equimap.get(shot, time=time_in, \ + R=R_in, Phi=Phi_in, Z=Z_in, \ + quantity='b_field_r') + end = time.time() + print() + print('time in equimap.get br =', end - start) + print() + print('oute.shape br =', oute.shape) + + # CALL EQUIMAP + start = time.time() + oute_noR = equimap.get(shot, time=time_in, \ + R=R_in, Phi=Phi_in, Z=Z_in, \ + quantity='b_field_r', no_ripple=True) + end = time.time() + print() + print('time in equimap.get br no ripple =', end - start) + print() + print('oute.shape br no ripple =', oute_noR.shape) + + # PLOTS + plt.figure() + plt.plot(time_in, oute[:, -1], \ + label='B_r at R={0}, Phi={1:.2f}, Z={2}'.format( \ + R_in[-1], Phi_in[-1], Z_in[-1])) + plt.plot(time_in, oute_noR[:, -1], \ + label='B_r no ripple at R={0}, Phi={1:.2f}, Z={2}'.format( \ + R_in[-1], Phi_in[-1], Z_in[-1])) + plt.legend() + plt.xlabel('Time [s]') + plt.ylabel('B_r [T]') + + plt.figure() + plt.plot(Phi_in, oute[int(0.5*oute.shape[0]), :], \ + label='B_r at t={0:.2f}, R={1}, Z={2}'.format( \ + time_in[int(0.5*oute.shape[0])], R_in[-1], Z_in[-1])) + plt.plot(Phi_in, oute_noR[int(0.5*oute.shape[0]), :], \ + label='B_r no ripple at t={0:.2f}, R={1}, Z={2}'.format( \ + time_in[int(0.5*oute.shape[0])], R_in[-1], Z_in[-1])) + plt.legend() + plt.xlabel('Phi [rad]') + plt.ylabel('B_r [T]') + + # CALL EQUIMAP + start = time.time() + outa = equimap.get(shot, time=time_in, \ + R=R_all_tot, Phi=np.zeros(R_all_tot.shape), Z=Z_all_tot, \ + quantity='b_field_r') + end = time.time() + print() + print('time in equimap.get br 2D =', end - start) + print() + outar = outa[int(0.5*outa.shape[0])].reshape((interp_points, interp_points)) + plt.figure() + plt.contourf(Rr, Zr, outar) + plt.colorbar() + arg_time = np.argmin(np.abs(out.time - time_in[int(0.5*outa.shape[0])])) + plt.plot(np.squeeze(out.time_slice[arg_time].boundary.outline.r), \ + np.squeeze(out.time_slice[arg_time].boundary.outline.z), \ + linewidth=2, color='red') + plt.plot(out.time_slice[arg_time].global_quantities.magnetic_axis.r, \ + out.time_slice[arg_time].global_quantities.magnetic_axis.z, \ + marker='+', color='red', markersize=20) + plt.xlabel('R [m]') + plt.ylabel('Z [m]') + plt.title('B_r t={0:.2f}'.format(time_in[int(0.5*outa.shape[0])])) + + # B_Z TEST + # -------- + Phi_in = np.linspace(0, 2*np.pi/18, 100) + R_in = np.full(Phi_in.shape, 3) + Z_in = np.full(Phi_in.shape, 0.2) + + # CALL EQUIMAP + start = time.time() + oute = equimap.get(shot, time=time_in, \ + R=R_in, Phi=Phi_in, Z=Z_in, \ + quantity='b_field_z') + end = time.time() + print() + print('time in equimap.get bz =', end - start) + print() + print('oute.shape bz =', oute.shape) + + # CALL EQUIMAP + start = time.time() + oute_noR = equimap.get(shot, time=time_in, \ + R=R_in, Phi=Phi_in, Z=Z_in, \ + quantity='b_field_z', no_ripple=True) + end = time.time() + print() + print('time in equimap.get bz no ripple =', end - start) + print() + print('oute.shape bz no ripple =', oute_noR.shape) + + # PLOTS + plt.figure() + plt.plot(time_in, oute[:, -1], \ + label='B_z at R={0}, Phi={1:.2f}, Z={2}'.format( \ + R_in[-1], Phi_in[-1], Z_in[-1])) + plt.plot(time_in, oute_noR[:, -1], \ + label='B_z no ripple at R={0}, Phi={1:.2f}, Z={2}'.format( \ + R_in[-1], Phi_in[-1], Z_in[-1])) + plt.legend() + plt.xlabel('Time [s]') + plt.ylabel('B_z [T]') + + plt.figure() + plt.plot(Phi_in, oute[int(0.5*oute.shape[0]), :], \ + label='B_z at t={0:.2f}, R={1}, Z={2}'.format( \ + time_in[int(0.5*oute.shape[0])], R_in[-1], Z_in[-1])) + plt.plot(Phi_in, oute_noR[int(0.5*oute.shape[0]), :], \ + label='B_z no ripple at t={0:.2f}, R={1}, Z={2}'.format( \ + time_in[int(0.5*oute.shape[0])], R_in[-1], Z_in[-1])) + plt.legend() + plt.xlabel('Phi [rad]') + plt.ylabel('B_z [T]') + + # CALL EQUIMAP + start = time.time() + outa = equimap.get(shot, time=time_in, \ + R=R_all_tot, Phi=np.zeros(R_all_tot.shape), Z=Z_all_tot, \ + quantity='b_field_z') + end = time.time() + print() + print('time in equimap.get bz 2D =', end - start) + print() + outar = outa[int(0.5*outa.shape[0])].reshape((interp_points, interp_points)) + plt.figure() + plt.contourf(Rr, Zr, outar) + plt.colorbar() + arg_time = np.argmin(np.abs(out.time - time_in[int(0.5*outa.shape[0])])) + plt.plot(np.squeeze(out.time_slice[arg_time].boundary.outline.r), \ + np.squeeze(out.time_slice[arg_time].boundary.outline.z), \ + linewidth=2, color='red') + plt.plot(out.time_slice[arg_time].global_quantities.magnetic_axis.r, \ + out.time_slice[arg_time].global_quantities.magnetic_axis.z, \ + marker='+', color='red', markersize=20) + plt.xlabel('R [m]') + plt.ylabel('Z [m]') + plt.title('B_z t={0:.2f}'.format(time_in[int(0.5*outa.shape[0])])) + + # B_TOR TEST + # ---------- + Phi_in = np.linspace(0, 2*np.pi/18, 100) + R_in = np.full(Phi_in.shape, 3) + Z_in = np.full(Phi_in.shape, 0) + + # CALL EQUIMAP + start = time.time() + oute = equimap.get(shot, time=time_in, \ + R=R_in, Phi=Phi_in, Z=Z_in, \ + quantity='b_field_tor') + end = time.time() + print() + print('time in equimap.get btor =', end - start) + print() + print('oute.shape btor =', oute.shape) + + # CALL EQUIMAP + start = time.time() + oute_noR = equimap.get(shot, time=time_in, \ + R=R_in, Phi=Phi_in, Z=Z_in, \ + quantity='b_field_tor', no_ripple=True) + end = time.time() + print() + print('time in equimap.get btor no ripple =', end - start) + print() + print('oute.shape btor no ripple =', oute_noR.shape) + + # PLOTS + plt.figure() + plt.plot(time_in, oute[:, -1], \ + label='B_tor at R={0}, Phi={1:.2f}, Z={2}'.format( \ + R_in[-1], Phi_in[-1], Z_in[-1])) + plt.plot(time_in, oute_noR[:, -1], \ + label='B_tor no ripple at R={0}, Phi={1:.2f}, Z={2}'.format( \ + R_in[-1], Phi_in[-1], Z_in[-1])) + plt.legend() + plt.xlabel('Time [s]') + plt.ylabel('B_tor [T]') + + plt.figure() + plt.plot(Phi_in, oute[int(0.5*oute.shape[0]), :], \ + label='B_tor at t={0:.2f}, R={1}, Z={2}'.format( \ + time_in[int(0.5*oute.shape[0])], R_in[-1], Z_in[-1])) + plt.plot(Phi_in, oute_noR[int(0.5*oute.shape[0]), :], \ + label='B_tor no ripple at t={0:.2f}, R={1}, Z={2}'.format( \ + time_in[int(0.5*oute.shape[0])], R_in[-1], Z_in[-1])) + plt.legend() + plt.xlabel('Phi [rad]') + plt.ylabel('B_tor [T]') + + # CALL EQUIMAP + start = time.time() + outa = equimap.get(shot, time=time_in, \ + R=R_all_tot, Phi=np.zeros(R_all_tot.shape), Z=Z_all_tot, \ + quantity='b_field_tor') + end = time.time() + print() + print('time in equimap.get btor 2D =', end - start) + print() + outar = outa[int(0.5*outa.shape[0])].reshape((interp_points, interp_points)) + plt.figure() + plt.contourf(Rr, Zr, outar) + plt.colorbar() + arg_time = np.argmin(np.abs(out.time - time_in[int(0.5*outa.shape[0])])) + plt.plot(np.squeeze(out.time_slice[arg_time].boundary.outline.r), \ + np.squeeze(out.time_slice[arg_time].boundary.outline.z), \ + linewidth=2, color='red') + plt.plot(out.time_slice[arg_time].global_quantities.magnetic_axis.r, \ + out.time_slice[arg_time].global_quantities.magnetic_axis.z, \ + marker='+', color='red', markersize=20) + plt.xlabel('R [m]') + plt.ylabel('Z [m]') + plt.title('B_tor t={0:.2f}'.format(time_in[int(0.5*outa.shape[0])])) + + plt.show() diff --git a/tofu/utils.py b/tofu/utils.py index 475b6a338..f3131bd68 100644 --- a/tofu/utils.py +++ b/tofu/utils.py @@ -643,8 +643,9 @@ def load_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, dsig=None, data=None, X=None, t0=None, dextra=None, - plot=True, plot_sig=None, plot_X=None, sharex=False, - bck=True, indch_auto=True, t=None, init=None): + plot=True, plot_sig=None, plot_X=None, + sharex=False, invertx=None, + bck=True, indch_auto=True, t=None, init=None, dR_sep=None): # ------------------- # import imas2tofu try: @@ -698,11 +699,36 @@ def load_from_imas(shot=None, run=None, user=None, tokamak=None, version=None, assert shot.size == 1 import tofu.mag as tfm plot = True + if invertx is None: + invertx = True + + multi = imas2tofu.MultiIDSLoader(shot=shot[0], run=run, user=user, + tokamak=tokamak, version=version, + ids='equilibrium') + equi = multi.to_Plasma2D() + if t is None: - t = np.r_[38] + # Get time in the middle of equilibrium time interval + t = equi.ddata['equilibrium.t']['data'][ + int(0.5*equi.ddata['equilibrium.t']['data'].size)] t = np.atleast_1d(t).ravel() - if init is None: - init = [[2.9],[0.],[0.]] + + equi_ind_t = np.abs(t - equi.ddata['equilibrium.t']['data']).argmin() + equi_ind_r_ext = np.argmax(equi.ddata['equilibrium.sep']['data'] + [equi_ind_t][0]) + equi_r_ext = equi.ddata['equilibrium.sep']['data'][ + equi_ind_t][0][equi_ind_r_ext] + equi_z_r_ext = equi.ddata['equilibrium.sep']['data'][ + equi_ind_t][1][equi_ind_r_ext] + + nbr_init = 10 + if dR_sep is not None: + r_init = [equi_r_ext + dR_sep]*nbr_init + else: + r_init = [equi_r_ext]*nbr_init + phi_init = [ii*2.*np.pi/nbr_init for ii in range(nbr_init)] + z_init = [equi_z_r_ext]*nbr_init + init_plt = [r_init, phi_init, z_init] if False: multi = imas2tofu.MultiIDSLoader(shot=shot[0], run=run, user=user, @@ -711,31 +737,92 @@ def load_from_imas(shot=None, run=None, user=None, tokamak=None, version=None, config = multi.to_Config(plot=False) else: import tofu.geom as tfg - config = tfg.utils.create_config('B2') + config = tfg.utils.create_config('B3') if config.nStruct > 1: config.set_colors_random() - trace = tfm.MagFieldLines(int(shot[0])).trace_mline(init, t, - direction='FWD', - length_line=None, - stp=None) + trace = tfm.MagFieldLines(int(shot[0])).trace_mline(init_plt, t, + direction='FWD', + length_line=35, + stp=None) + trace_rev = tfm.MagFieldLines( + int(shot[0])).trace_mline(init_plt, t, + direction='REV', + length_line=35, + stp=None) + refpt = np.r_[2.4,0.] - dax = config.plot_phithetaproj_dist(refpt) + dax = config.plot_phithetaproj_dist(refpt, invertx=invertx) + + if init is not None: + trace_init = tfm.MagFieldLines( + int(shot[0])).trace_mline(init, t, + direction='FWD', + length_line=35, + stp=None) + trace_init_rev = tfm.MagFieldLines( + int(shot[0])).trace_mline(init, t, + direction='REV', + length_line=35, + stp=None) + trace_init[0] = trace_init[0] + trace_init_rev[0] + + for kk in range(0, len(trace_init[0])): + phi_init = np.arctan2(np.sin(trace_init[0][kk]['p']), + np.cos(trace_init[0][kk]['p'])) + theta_init = np.arctan2(trace_init[0][kk]['z']-refpt[1], + trace_init[0][kk]['r']-refpt[0]) + indnan = ((np.abs(np.diff(phi_init)) > np.pi) + | (np.abs(np.diff(theta_init)) + > np.pi)).nonzero()[0] + 1 + dax['dist'][0].plot(np.insert(phi_init, indnan, np.nan), + np.insert(theta_init, indnan, np.nan), + linewidth=3, color='red') + dax['cross'][0].plot(trace_init[0][kk]['r'], + trace_init[0][kk]['z'], + linewidth=3, color='red') + x = trace_init[0][kk]['r']*np.cos(trace_init[0][kk]['p']) + y = trace_init[0][kk]['r']*np.sin(trace_init[0][kk]['p']) + dax['hor'][0].plot(x, y, linewidth=3, color='red') + alpha_mag_lines = 0.7 + else: + alpha_mag_lines = 1. + for ii in range(0,len(trace)): + # Concatenate trace lists + trace[ii] = trace[ii] + trace_rev[ii] for jj in range(0,len(trace[ii])): lab = r't = %s s'%str(t[ii]) - phi = np.arctan2(np.sin(trace[ii][jj]['p']), np.cos(trace[ii][jj]['p'])) - theta = np.arctan2(trace[ii][jj]['z']-refpt[1], trace[ii][jj]['r']-refpt[0]) + phi = np.arctan2(np.sin(trace[ii][jj]['p']), + np.cos(trace[ii][jj]['p'])) + theta = np.arctan2(trace[ii][jj]['z']-refpt[1], + trace[ii][jj]['r']-refpt[0]) # insert nans for clean periodicity indnan = ((np.abs(np.diff(phi)) > np.pi) | (np.abs(np.diff(theta)) > np.pi)).nonzero()[0] + 1 dax['dist'][0].plot(np.insert(phi, indnan, np.nan), np.insert(theta, indnan, np.nan), - label=lab) + label=lab, alpha=alpha_mag_lines) dax['cross'][0].plot(trace[ii][jj]['r'], trace[ii][jj]['z'], - label=lab) + label=lab, alpha=alpha_mag_lines) x = trace[ii][jj]['r']*np.cos(trace[ii][jj]['p']) y = trace[ii][jj]['r']*np.sin(trace[ii][jj]['p']) - dax['hor'][0].plot(x, y, label=lab) + dax['hor'][0].plot(x, y, label=lab, alpha=alpha_mag_lines) + + dax['cross'][0].plot(equi.ddata['equilibrium.sep'][ + 'data'][equi_ind_t][0], + equi.ddata['equilibrium.sep'][ + 'data'][equi_ind_t][1], + linestyle='-.', color='k', alpha=0.8) + dax['cross'][0].plot(multi.get_data('equilibrium')['strike0'][ + equi_ind_t][0], + multi.get_data('equilibrium')['strike0'][ + equi_ind_t][1], '+', color='k', markersize=10) + dax['cross'][0].plot(multi.get_data('equilibrium')['strike1'][ + equi_ind_t][0], + multi.get_data('equilibrium')['strike1'][ + equi_ind_t][1], '+', color='k', markersize=10) + dax['t'][0].figure.suptitle('Shot {0}, t = {1:6.3f} s' + .format(shot[0], t[0])) return dax diff --git a/tofu/version.py b/tofu/version.py index 31088c3e2..1bd359aa2 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.1-264-g1f49781' +__version__ = '1.4.1-303-gf837259' diff --git a/tofuplot.py b/tofuplot.py index f8e3a2fb6..7b8ed9906 100755 --- a/tofuplot.py +++ b/tofuplot.py @@ -19,7 +19,7 @@ if istofugit: # Make sure we load the corresponding tofu - sys.path.insert(1,_HERE) + sys.path.insert(1, _HERE) import tofu as tf from tofu.imas2tofu import MultiIDSLoader _ = sys.path.pop(1) @@ -49,7 +49,7 @@ _VERSION = '3' _LIDS_DIAG = MultiIDSLoader._lidsdiag _LIDS_PLASMA = tf.imas2tofu.MultiIDSLoader._lidsplasma -_LIDS = _LIDS_DIAG + _LIDS_PLASMA +_LIDS = _LIDS_DIAG + _LIDS_PLASMA + ['magfieldlines'] _T0 = 'IGNITRON' _SHAREX = False _BCK = True @@ -75,7 +75,7 @@ def call_tfloadimas(shot=None, run=_RUN, user=_USER, tokamak=_TOKAMAK, version=_VERSION, ids=None, quantity=None, X=None, t0=_T0, sharex=_SHAREX, indch=None, indch_auto=None, - background=_BCK): + background=_BCK, t=None, dR_sep=None, init=None): lidspla = [ids_ for ids_ in ids if ids_ in _LIDS_PLASMA] if t0.lower() == 'none': @@ -85,7 +85,8 @@ def call_tfloadimas(shot=None, run=_RUN, user=_USER, tokamak=tokamak, version=version, ids=ids, indch=indch, indch_auto=indch_auto, plot_sig=quantity, plot_X=X, - t0=t0, plot=True, sharex=sharex, bck=background) + t0=t0, plot=True, sharex=sharex, bck=background, + t=t, dR_sep=dR_sep, init=init) plt.show(block=True) @@ -125,11 +126,11 @@ def _str2bool(v): msg = 'username of the DB where the datafile is located' parser.add_argument('-u','--user',help=msg, required=False, default=_USER) msg = 'tokamak name of the DB where the datafile is located' - parser.add_argument('-t','--tokamak',help=msg, required=False, + parser.add_argument('-tok', '--tokamak', help=msg, required=False, default=_TOKAMAK) - parser.add_argument('-r','--run',help='run number', + parser.add_argument('-r', '--run', help='run number', required=False, type=int, default=_RUN) - parser.add_argument('-v','--version',help='version number', + parser.add_argument('-v', '--version', help='version number', required=False, type=str, default=_VERSION) msg = "ids from which to load diagnostics data, can be:\n%s"%repr(_LIDS) @@ -143,12 +144,21 @@ def _str2bool(v): nargs='+', default=None) parser.add_argument('-t0', '--t0', type=str, required=False, help='Reference time event setting t = 0', default=_T0) + parser.add_argument('-t', '--t', type=float, required=False, + help='Input time when needed') + parser.add_argument('-dR_sep', '--dR_sep', type=float, required=False, + help='Distance to separatrix from r_ext to plot' + + ' 10 magnetic field lines') + parser.add_argument('-init', '--init', type=float, required=False, nargs=3, + help='Manual coordinates of point that a RED magnetic' + + ' field line will cross on graphics,' + + ' give coordinates as: R [m], Phi [rad], Z [m]') parser.add_argument('-ich', '--indch', type=int, required=False, help='indices of channels to be loaded', nargs='+', default=None) parser.add_argument('-ichauto', '--indch_auto', type=bool, required=False, - help='automatically determine indices of channels to be loaded', - default=True) + help='automatically determine indices of' + + ' channels to be loaded', default=True) parser.add_argument('-sx', '--sharex', type=_str2bool, required=False, help='Should X axis be shared between diagnostics ids ?', default=_SHAREX, const=True, nargs='?')