Coverage for fiqus/plotters/PlotPythonConductorAC.py: 8%

1import os 

2import numpy as np 

3import re 

4import matplotlib.pyplot as plt 

5from matplotlib.ticker import FuncFormatter 

6from matplotlib.animation import FuncAnimation 

7import pandas as pd 

8from scipy import integrate, interpolate 

9from ruamel.yaml import YAML 

10import subprocess 

11 

12# from fiqus.utils.Utils import FilesAndFolders as Util 

13# from fiqus.data.DataFiQuSConductorAC_Strand import CACStrandSolve, CACStrandPostproc, CACStrandMesh, CACStrandGeometry 

14 

15def create_non_overwriting_filepath(folder_path, base_name, extension, overwrite): 

16 """ 

17 Creates a filepath that does not overwrite any existing files. 

18 

19 This function checks if a file already exists at the specified filepath. If the file exists and `overwrite` is False,  

20 it modifies the filepath to create a new file instead of overwriting the existing one.  

21 If `overwrite` is True or the file does not exist, it returns the filepath as it is. 

22 

23 Parameters 

24 ---------- 

25 folder_path : str 

26 The path to the folder where the file will be created. 

27 base_name : str 

28 The base name of the file. 

29 extension : str 

30 The extension of the file. 

31 overwrite : bool, optional 

32 If True, the function will overwrite an existing file. If False, the function will modify the filepath to avoid overwriting. Defaults to False. 

33 

34 Returns 

35 ------- 

36 str 

37 The final filepath. If `overwrite` is False and a file already exists at the original filepath, this will be a new filepath that does not overwrite any existing files. 

38 """ 

39 if os.path.exists(os.path.join(folder_path, base_name+extension)) and not overwrite: 

40 counter = 1 

41 new_name = base_name + f"_{counter}" + extension 

42 while os.path.exists(os.path.join(folder_path, new_name)): 

43 new_name = base_name + f"_{counter}" + extension 

44 counter += 1 

45 return os.path.join(folder_path, new_name) 

46 

47 return os.path.join(folder_path, base_name+extension) 

48 

49class YamlWrapper: 

50 """ 

51 A wrapper class for YAML data that allows accessing dictionary data using dot notation. 

52 """ 

53 def __init__(self, data): 

54 for key, value in data.items(): 

55 if isinstance(value, dict): 

56 value = YamlWrapper(value) 

57 self.__dict__[key] = value 

58 

59 def __getattr__(self, item): 

60 return self.__dict__.get(item) 

61 

62 def __setattr__(self, key, value): 

63 self.__dict__[key] = value 

64 

65def load_yaml(file_path): 

66 """ 

67 Load a YAML file and return the data as a YamlWrapper object. This enables accessing the data using dot notation (e.g., data.key.subkey), without needing a predefined pydantic model, allowing for better backwards compatibility. 

68 """ 

69 yaml = YAML() 

70 with open(file_path, 'r') as file: 

71 data = yaml.load(file) 

72 return YamlWrapper(data) 

73 

74def is_latex_installed(): 

75 """ 

76 Check if LaTeX is installed on the system. 

77 """ 

78 try: 

79 subprocess.run(['pdflatex', '--version'], stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL) 

80 return True 

81 except subprocess.CalledProcessError: 

82 return False 

83 except FileNotFoundError: 

84 return False 

85 

86class SimulationData: 

87 """ 

88 Class used to store and manage data from a single simulation. 

89 

90 This class is responsible for loading and organizing the data from a single simulation.  

91 It stores the data in various attributes and provides methods for retrieving and processing the data. 

92 

93 """ 

94 def __init__(self, model_data_output_path, geometry_name, mesh_name, solution_name) -> None: 

95 self.model_data_output_path = model_data_output_path # This is the path to the folder where the model output data is stored (e.g. geometries) 

96 self.geometry_name = geometry_name # Name of the geometry folder 

97 self.mesh_name = mesh_name # Name of the mesh folder  

98 self.solution_name = solution_name # Name of the solution folder 

99 

100 # Organize the folders: 

101 self.geometry_folder = os.path.join(self.model_data_output_path, geometry_name) # Path to the geometry folder 

102 self.mesh_folder = os.path.join(self.geometry_folder, mesh_name) # Path to the mesh folder 

103 self.solution_folder = os.path.join(self.mesh_folder, solution_name) # Path to the solution folder 

104 

105 # Store the YAML input-files in a data model, fdm: 

106 self.geometry, self.mesh, self.solve = self.retrieve_fiqusDataModel() 

107 

108 # Store losses, simulation time and check if the simulation crashed: 

109 temp_file_path = os.path.join(self.solution_folder, 'test_temporary') 

110 loss_file = [f for f in os.listdir(temp_file_path) if f.startswith('power') and f.endswith('.txt')][0] 

111 self.power_columns = ['Time', 'FilamentLoss', 'CouplingLoss', 'EddyLoss', 'TotalLoss', 'CouplingLoss_dyn', 'TotalLoss_dyn'] # Only in the case dynamic correction is used, must be changed later 

112 self.power = pd.read_csv(os.path.join(self.solution_folder, 'test_temporary', loss_file), sep = ' ', names=self.power_columns) # Store instantaneous losses as pandas dataframe 

113 self.crash = True if 'crash_report.txt' in os.listdir(temp_file_path) else False 

114 

115 # Add a row of zeros at the beginning of the dataframe to account for the initial condition: 

116 self.power = pd.concat([pd.DataFrame({col: 0 for col in self.power_columns}, index=[0]), self.power]).reset_index(drop=True) 

117 # Integrate the losses to obtain the cumulative power and the total power per cycle: 

118 self.cumulative_power, self.total_power_per_cycle = self.integrate_power() # Store cumulative power and total cumulative power per cycle 

119 # Store simulation time: 

120 try: 

121 with open(os.path.join(self.solution_folder, 'test_temporary', 'simulation_time.txt'), 'r') as f: 

122 self.simulation_time = float(f.readline().strip()) 

123 except: 

124 self.simulation_time = None # If the simulation time file does not exist, the simulation has not finished running. 

125 

126 # Store the rest of the post-processing data: 

127 self.time = self.power['Time'] 

128 self.instantaneous_temperature = self.load_standard_data(os.path.join(temp_file_path, 'temperature.txt'), 1, add_initial_zero=True) 

129 self.temperature = self.load_standard_data(os.path.join(temp_file_path, 'temperature.txt'), 2, add_initial_zero=True) 

130 self.I_transport = self.load_standard_data(os.path.join(temp_file_path, 'I_transport.txt'), 1) 

131 self.V_transport = self.load_standard_data(os.path.join(temp_file_path, 'V_transport.txt'), 1) 

132 self.hs_val = self.load_standard_data(os.path.join(temp_file_path, 'hs_val.txt'), 1) 

133 self.magn_fil = self.load_standard_data(os.path.join(temp_file_path, 'magn_fil.txt'), [1, 2, 3]) 

134 self.magn_matrix = self.load_standard_data(os.path.join(temp_file_path, 'magn_matrix.txt'), [1, 2, 3]) 

135 self.I = self.load_standard_data(os.path.join(temp_file_path, 'I.txt'), 1, len(self.time)) 

136 self.V = self.load_standard_data(os.path.join(temp_file_path, 'V.txt'), 1, len(self.time)) 

137 self.I_integral = self.load_standard_data(os.path.join(temp_file_path, 'I_integral.txt'), 1, len(self.time)) 

138 self.I_abs_integral = self.load_standard_data(os.path.join(temp_file_path, 'I_abs_integral.txt'), 1, len(self.time)) 

139 self.magnetic_energy_internal = self.load_standard_data(os.path.join(temp_file_path, 'magnetic_energy_internal.txt'), 1) 

140 self.Ip = self.load_standard_data(os.path.join(temp_file_path, 'Ip.txt'), 1, len(self.time)) 

141 self.Vp = self.load_standard_data(os.path.join(temp_file_path, 'Vp.txt'), 1, len(self.time)) 

142 

143 

144 def load_standard_data(self, file_path, columns, reshape = None, add_initial_zero = False): 

145 """ 

146 There are many output .txt-files with similar format. This function loads the data from one of these files and returns it as a numpy array. 

147 If the file does not exist, None is returned without raising an error. 

148 """ 

149 try: 

150 data = np.loadtxt(file_path, comments='#', usecols=columns) 

151 if reshape: 

152 data = data.reshape(-1, reshape).T 

153 except IOError: 

154 return None 

155 

156 if add_initial_zero: 

157 if len(data.shape) == 1: 

158 data = np.insert(data, 0, 0) 

159 else: 

160 zeros = np.zeros((1, data.shape[1])) 

161 data = np.vstack((zeros, data)) 

162 return data 

163 

164 def retrieve_fiqusDataModel(self): 

165 """ 

166 This function reads the YAML input-files for geometry, mesh and solve and stores them in three dictionaries which are returned. 

167 This function is to be called only once, when the object is created. 

168 """ 

169 geometry_dataModel = load_yaml(os.path.join(self.geometry_folder, 'geometry.yaml')) 

170 mesh_dataModel = load_yaml(os.path.join(self.mesh_folder, 'mesh.yaml')) 

171 solution_dataModel = load_yaml(os.path.join(self.solution_folder, 'solve.yaml')) 

172 

173 return geometry_dataModel, mesh_dataModel, solution_dataModel 

174 

175 def integrate_power(self): 

176 """  

177 This function integrates the instantaneous power over time to obtain the cumulative power.  

178 It also calculates the total cumulative power per cycle. 

179 The cumulative power is returned as a pandas dataframe and the total cumulative power per cycle is returned as a dictionary. 

180 """ 

181 find_closest_idx = lambda arr, val: np.abs(arr - val).argmin() 

182 

183 t = np.array(self.power['Time']) 

184 t_final = t[-1] 

185 t_init = find_closest_idx(t, t_final/2) 

186 

187 cumulative_power = pd.DataFrame(columns= self.power_columns) 

188 total_power_per_cycle = {} 

189 

190 cumulative_power['Time'] = self.power["Time"] 

191 for column in self.power_columns[1:]: 

192 cumulative_power[column] = np.insert(integrate.cumulative_trapezoid(self.power[column], t), 0, 0) 

193 total_power_per_cycle[column] = 2 * (cumulative_power[column].iloc[-1]-cumulative_power[column].iloc[t_init]) # / (np.pi*matrix_radius**2 * loss_factor) # Why do we divide by pi*matrix_radius**2*loss_factor? 

194 

195 return cumulative_power, total_power_per_cycle 

196 

197 def plot_instantaneous_power(self, show:bool = True, title:str = "Power", save_plot:bool = False, save_folder_path:str = None, save_file_name:str = None, overwrite:bool = False): 

198 plt.figure() 

199 plt.plot(self.power['Time'], self.power[self.power_columns[1:]] , label = self.power_columns[1:]) 

200 plt.xlabel('Time [s]') 

201 plt.ylabel('Power [W/m]') 

202 plt.legend() 

203 

204 # Configure title: 

205 # Class attributes can be accessed by using '<< ... >>' in the title string. 

206 commands = re.findall('<<(.*?)>>', title) 

207 for c in commands: 

208 title = title.replace(f"<<{c}>>", str(eval('self.'+c))) 

209 plt.title(title) 

210 

211 

212 if save_plot: # Save the plot 

213 filePath = create_non_overwriting_filepath(save_folder_path, save_file_name, '.png', overwrite) 

214 plt.savefig(filePath) 

215 

216 if show: 

217 plt.show() 

218 else: 

219 plt.close() 

220 

221class PlotPython: 

222 """ 

223 This class loads and stores the data from the simulations specified in a csv file and can apply various postprocessing operations on the data.  

224 The data from each simulation is saved as a SimulationData object which is subsequently stored in a list in this class. 

225 """ 

226 def __init__(self, fdm, csv_filename = None, lossMap_gridData_folder = None, inputs_folder_path='', outputs_folder_path='') -> None: 

227 self.fdm = fdm 

228 self.inputs_folder_path = inputs_folder_path # This is the path to the folder where the input data is stored 

229 self.model_data_output_path = outputs_folder_path # This is the path to the folder where the model output data is stored (e.g. geometries) 

230 self.outputs_folder_path = os.path.join(outputs_folder_path, fdm.magnet.postproc.batch_postproc.output_folder) # This is the path to the folder where the postprocessed data is written 

231 

232 

233 if not os.path.exists(self.outputs_folder_path): 

234 os.makedirs(self.outputs_folder_path) 

235 

236 if csv_filename is not None: 

237 try: 

238 self.input_csv = pd.read_csv(os.path.join(self.inputs_folder_path, f'{csv_filename}.csv')) # Read the csv file with the input data 

239 except: 

240 raise FileNotFoundError(f'No csv file with the name {fdm.magnet.postproc.batch_postproc.postProc_csv}.csv was found in the inputs folder.') 

241 

242 self.simulation_collection = self.retrieve_simulation_data() 

243 

244 self.avg_simulation_time = np.mean([sd.simulation_time for sd in self.simulation_collection]) 

245 self.total_simulation_time = np.sum([sd.simulation_time for sd in self.simulation_collection]) 

246 

247 print('Number of simulations considered: ', len(self.simulation_collection) ) 

248 print('Average simulation time: ', self.avg_simulation_time, 's') 

249 print('Total simulation time: ', self.total_simulation_time, 's') 

250 

251 

252 

253 elif lossMap_gridData_folder is not None: 

254 self.input_csv = None 

255 self.simulation_collection = None 

256 

257 self.totalLoss_gridData = self.load_lossMap_gridData('TotalLoss', lossMap_gridData_folder) 

258 self.filamentLoss_gridData = self.load_lossMap_gridData('FilamentLoss', lossMap_gridData_folder) 

259 self.eddyLoss_gridData = self.load_lossMap_gridData('EddyLoss', lossMap_gridData_folder) 

260 self.couplingLoss_gridData = self.load_lossMap_gridData('CouplingLoss', lossMap_gridData_folder) 

261 else: 

262 raise ValueError('No input data specified. Either a csv file or a folder with loss map grid data must be provided.') 

263 

264 

265 

266 

267 def retrieve_simulation_data(self): 

268 """ 

269 This function iterates over the input CSV-file (specifying which simulations to postprocess) and returns a list of SimulationData objects 

270 containing all the simulation data. If no CSV-file is specified, the data from the single simulation specified in the input YAML-file is returned. 

271 """ 

272 if self.input_csv is not None: 

273 simulationCollection = [] 

274 for index, row in self.input_csv.iterrows(): 

275 if pd.isna(row['input.run.geometry']) and pd.isna(row['input.run.mesh']) and pd.isna(row['input.run.solution']): 

276 continue 

277 geometry_name = 'Geometry_'+str(row['input.run.geometry']) 

278 mesh_name = 'Mesh_'+str(row['input.run.mesh']) 

279 

280 if isinstance(row['input.run.solution'], float) and row['input.run.solution'].is_integer(): 

281 solution_name = 'Solution_'+str(int(row['input.run.solution'])) 

282 else: 

283 solution_name = 'Solution_'+str(row['input.run.solution']) 

284 

285 # Check if the row refers to a valid simulation by checking if the solution folder exists: 

286 # solution_folder = os.path.join(os.getcwd(), 'tests', '_outputs', self.fdm.general.magnet_name, geometry_name, mesh_name, solution_name) 

287 solution_folder = os.path.join(self.model_data_output_path, geometry_name, mesh_name, solution_name) 

288 if os.path.exists(solution_folder): # If the solution folder exists, add the simulation to the simulationCollection 

289 sd = SimulationData(self.model_data_output_path, geometry_name, mesh_name, solution_name) 

290 if sd.simulation_time is not None: # Only add the simulation if it has finished running (and therefore has written the simulation time to a file) 

291 simulationCollection.append(sd) 

292 else: 

293 simulationCollection = [SimulationData(self.model_data_output_path, 'Geometry_'+self.fdm.run.geometry, 'Mesh_'+self.fdm.run.mesh, 'Solution_'+self.fdm.run.solution)] 

294 

295 return self.sort_simulationCollection(self.filter_simulationCollection(simulationCollection)) 

296 

297 def filter_simulationCollection(self, simulationCollection): 

298 """  

299 This function is used to filter the simulationCollection based on the filter criterion specified in the yaml input file. 

300 An example of a filter criterion is '<<solve.source_parameters.sine.frequency>> == 18', which will disregard all simulations with frequency != 18Hz. 

301 """ 

302 if self.fdm.magnet.postproc.batch_postproc.filter.apply_filter: 

303 filter_criterion = self.fdm.magnet.postproc.batch_postproc.filter.filter_criterion 

304 class_params = re.findall('<<(.*?)>>', filter_criterion) 

305 for cp in class_params: 

306 filter_criterion = filter_criterion.replace(f"<<{cp}>>", 'sd.'+cp) 

307 filtering_function = eval(f'lambda sd: {filter_criterion}') 

308 return list(filter(filtering_function, simulationCollection)) 

309 else: 

310 return simulationCollection 

311 

312 def sort_simulationCollection(self, simulationCollection): 

313 """  

314 This function is used to sort the simulationCollection based on the sort key specified in the yaml input file. 

315 An example of a sort key is 'solve.source_parameters.sine.frequency', which will sort the simulations based on frequency. 

316 """ 

317 if self.fdm.magnet.postproc.batch_postproc.sort.apply_sort: 

318 sorting_function = eval(f'lambda sd: sd.{self.fdm.magnet.postproc.batch_postproc.sort.sort_key}') 

319 return sorted(simulationCollection, key=sorting_function) 

320 else: 

321 return simulationCollection 

322 

323 def lossMap_createGridData(self, lossType = 'TotalLoss', x_val_to_include = None, y_val_to_include = None): 

324 """ 

325 This function creates the grid data needed for the loss map, based on the yaml input file. 

326 Given a collection of simulations it interpolates the loss data between the datapoints to a grid and returns the grid data. 

327 """ 

328 lm = self.fdm.magnet.postproc.batch_postproc.loss_map 

329 

330 # Extract data from simulation collection and normalize 

331 x_arr = np.array([eval('sd.'+lm.x_val)/lm.x_norm for sd in self.simulation_collection]) 

332 y_arr = np.array([eval('sd.'+lm.y_val)/lm.y_norm for sd in self.simulation_collection]) 

333 loss = np.array([sd.total_power_per_cycle[lossType]/lm.loss_norm for sd in self.simulation_collection]) 

334 

335 # Logarithmic scaling 

336 if lm.x_log: x_arr = np.log10(x_arr) 

337 if lm.y_log: y_arr = np.log10(y_arr) 

338 if lm.loss_log: loss = np.log10(loss) 

339 

340 x_arr_interpolated = np.linspace(min(x_arr), max(x_arr), lm.x_steps) 

341 y_arr_interpolated = np.linspace(min(y_arr), max(y_arr), lm.y_steps) 

342 # Insert specific values to the grid if they are not already included (useful for cross sections) 

343 if x_val_to_include is not None and x_val_to_include not in x_arr_interpolated: 

344 x_arr_interpolated = np.insert(x_arr_interpolated, np.where(x_arr_interpolated > x_val_to_include)[0][0], x_val_to_include) 

345 if y_val_to_include is not None and y_val_to_include not in y_arr_interpolated: 

346 y_arr_interpolated = np.insert(y_arr_interpolated, np.where(y_arr_interpolated > y_val_to_include)[0][0], y_val_to_include) 

347 

348 # Create grid 

349 X, Y = np.meshgrid(x_arr_interpolated, y_arr_interpolated, indexing='ij') 

350 gridPoints = np.c_[X.ravel(), Y.ravel()] 

351 dataPoints = np.c_[x_arr, y_arr] 

352 

353 # Interpolate the simulation data onto the grid 

354 V = interpolate.griddata( 

355 dataPoints, 

356 loss, 

357 gridPoints, 

358 method='linear' # Cubic produces cleaner plots. Any incentive to go back to linear? 

359 ).reshape(X.shape) 

360 

361 return X, Y, V, dataPoints 

362 

363 def save_lossMap_gridData(self, save_folder_name = 'lossMap_gridData'): 

364 """ 

365 This function calls the lossMap_createGridData function and saves the grid data. 

366 """ 

367 lm = self.fdm.magnet.postproc.batch_postproc.loss_map 

368 

369 lossTypes = ['TotalLoss', 'FilamentLoss', 'EddyLoss', 'CouplingLoss', 'CouplingLoss_dyn', 'TotalLoss_dyn'] # Only in the case dynamic correction is used, must be changed later 

370 # 1) Create a folder to store the output files 

371 gridData_folder_path = create_non_overwriting_filepath(self.outputs_folder_path, save_folder_name, '', self.fdm.run.overwrite) 

372 if not os.path.exists(gridData_folder_path): os.makedirs(gridData_folder_path) 

373 # 2) Create the grid data for each loss type and save it 

374 for lossType in lossTypes: 

375 X, Y, V, _ = self.lossMap_createGridData(lossType) 

376 if lm.x_log: X = np.power(10, X) 

377 if lm.y_log: Y = np.power(10, Y) 

378 if lm.loss_log: V = np.power(10, V) 

379 np.savetxt(os.path.join(gridData_folder_path, f'{lossType}.txt'), np.column_stack((X.ravel(), Y.ravel(), V.ravel())), delimiter=' ', header=f'{lm.x_val} {lm.y_val} {lossType}', comments='') 

380 

381 def load_lossMap_gridData(self, lossType = 'TotalLoss', save_folder_name = 'lossMap_gridData'): 

382 """ 

383 This function loads the grid data for a given loss type. 

384 """ 

385 lm = self.fdm.magnet.postproc.batch_postproc.loss_map 

386 gridData_folder_path = os.path.join(self.inputs_folder_path, save_folder_name) 

387 

388 if not os.path.exists(gridData_folder_path): 

389 raise FileNotFoundError(f'The folder {gridData_folder_path} does not exist.') 

390 

391 X, Y, V = np.loadtxt(os.path.join(gridData_folder_path, f'{lossType}.txt'), unpack=True, skiprows=1) 

392 

393 if lm.x_log: X = np.log10(X) 

394 if lm.y_log: Y = np.log10(Y) 

395 if lm.loss_log: V = np.log10(V) 

396 

397 # Get the unique counts of X and Y 

398 unique_X = np.unique(X) 

399 unique_Y = np.unique(Y) 

400 

401 # Reshape the data 

402 X = X.reshape((len(unique_X), len(unique_Y))) 

403 Y = Y.reshape((len(unique_X), len(unique_Y))) 

404 V = V.reshape((len(unique_X), len(unique_Y))) 

405 

406 return X, Y, V 

407 

408 # def add_value_to_gridData(self, X, Y, V, x_val = None, y_val = None): 

409 # """ 

410 # This function adds a value to the grid data. 

411 # Steps:  

412 # 1) Revert the grid data to 1D arrays 

413 # 2) Add x or y value or both to the arrays 

414 # 3) Reshape the arrays back to grid data, interpolating the loss to the new grid points 

415 # """ 

416 # gridPoints  

417 

418 def save_magnetization(self): 

419 """ 

420 This function saves the magnetization data for all simulations in the simulation collection. 

421 """ 

422 magnetization_folder_path = create_non_overwriting_filepath(self.outputs_folder_path, 'magnetization', '', self.fdm.run.overwrite) 

423 if not os.path.exists(magnetization_folder_path): os.makedirs(magnetization_folder_path) 

424 for sd in self.simulation_collection: 

425 magnetization = sd.magn_fil + sd.magn_matrix 

426 magnetization = np.c_[sd.time, magnetization] 

427 np.savetxt(os.path.join(magnetization_folder_path, f'magn_f{sd.solve.source_parameters.sine.frequency}_b{sd.solve.source_parameters.sine.field_amplitude}_I{sd.solve.source_parameters.sine.current_amplitude}.txt'), magnetization, delimiter=' ', header='t x y z', comments='') 

428 

429 

430 

431 

432 

433 

434 def lossMap_crossSection(self, slice_value, axis_to_cut = 'x'): 

435 """ 

436 This function returns the data corresponding to a cross section of the loss map, for all loss types. 

437 Given an axis and a value, it sweeps the other axis for the closest value and returns the data. 

438 Example: Given slice value 0 and axis x, it returns the data for the cross section at x = 0. 

439 """ 

440 

441 lm = self.fdm.magnet.postproc.batch_postproc.loss_map 

442 if axis_to_cut == 'x': 

443 x_val_to_include = slice_value 

444 y_val_to_include = None 

445 elif axis_to_cut == 'y': 

446 x_val_to_include = None 

447 y_val_to_include = slice_value 

448 X, Y, V, dataPoints = self.lossMap_createGridData('TotalLoss', x_val_to_include, y_val_to_include) 

449 _,_,FilamentLoss, _ = self.lossMap_createGridData('FilamentLoss', x_val_to_include, y_val_to_include) 

450 _,_,EddyLoss, _ = self.lossMap_createGridData('EddyLoss', x_val_to_include, y_val_to_include) 

451 _,_,CouplingLoss, _ = self.lossMap_createGridData('CouplingLoss', x_val_to_include, y_val_to_include) 

452 

453 

454 if axis_to_cut == 'x': 

455 index = np.abs(X[:, 0] - slice_value).argmin() 

456 slice_vals = Y[index, :] 

457 

458 elif axis_to_cut == 'y': 

459 index = np.abs(Y[0, :] - slice_value).argmin() 

460 slice_vals = X[:, index] 

461 

462 # Extract the loss values for the constant frequency across all applied fields 

463 totalLoss = V[index, :] if axis_to_cut == 'x' else V[:, index] 

464 filamentLoss = FilamentLoss[index, :] if axis_to_cut == 'x' else FilamentLoss[:, index] 

465 eddyLoss = EddyLoss[index, :] if axis_to_cut == 'x' else EddyLoss[:, index] 

466 couplingLoss = CouplingLoss[index, :] if axis_to_cut == 'x' else CouplingLoss[:, index] 

467 

468 return slice_vals, totalLoss, filamentLoss, eddyLoss, couplingLoss 

469 

470 def plot_lossMap_crossSection(self): 

471 """ 

472 This function calls the lossMap_crossSection function and plots the data it returns, which is the loss for all values of one axis, given a constant value of the other axis. 

473 """ 

474 

475 if is_latex_installed(): 

476 plt.rcParams['text.usetex'] = True 

477 plt.rcParams['font.family'] = 'times' 

478 plt.rcParams['font.size'] = 20 

479 

480 lm = self.fdm.magnet.postproc.batch_postproc.loss_map 

481 slice_value = lm.cross_section.cut_value 

482 axis_to_cut = lm.cross_section.axis_to_cut 

483 

484 if (lm.x_log and axis_to_cut == 'x') or (lm.y_log and axis_to_cut == 'y'): 

485 slice_value = np.log10(slice_value) 

486 

487 slice_vals, totalLoss, filamentLoss, eddyLoss, couplingLoss = self.lossMap_crossSection(slice_value, axis_to_cut = axis_to_cut) 

488 

489 def log_formatter(x, pos): 

490 """ 

491 Format the tick labels on the plot. 

492 """ 

493 return f"$10^{{{int(x)}}}$" 

494 

495 # Plot the loss with respect to applied field for the constant frequency 

496 fig, ax = plt.subplots(figsize=(8, 6)) 

497 ax.plot(slice_vals, totalLoss, label=f'Total Loss') 

498 ax.plot(slice_vals, filamentLoss, label=f'Filament Loss') 

499 ax.plot(slice_vals, eddyLoss, label=f'Eddy Loss') 

500 ax.plot(slice_vals, couplingLoss, label=f'Coupling Loss') 

501 

502 tick_formatter = FuncFormatter(log_formatter) 

503 if lm.x_log and axis_to_cut == 'y' or lm.y_log and axis_to_cut == 'x': 

504 ax.xaxis.set_major_formatter(tick_formatter) 

505 if lm.loss_log: 

506 ax.yaxis.set_major_formatter(tick_formatter) 

507 

508 

509 title = lm.cross_section.title.replace('<<cut_value>>', str(round(10**slice_value, 3))) 

510 ax.set_title(title) 

511 ax.set_xlabel(lm.ylabel if axis_to_cut == 'x' else lm.xlabel) 

512 ax.set_ylabel(lm.cross_section.ylabel) 

513 ax.legend() 

514 

515 # np.savetxt(os.path.join(self.outputs_folder_path, 'lossMaps_cut_0p2T_0A.txt'), np.column_stack((10**slice_vals, 10**totalLoss, 10**eddyLoss, 10**couplingLoss, 10**filamentLoss)), delimiter=' ', header='f total eddy coupling filament', comments='') 

516 

517 if lm.cross_section.save_plot: 

518 filePath = create_non_overwriting_filepath(self.outputs_folder_path, lm.cross_section.filename, '.png', self.fdm.run.overwrite) 

519 plt.savefig(filePath) 

520 

521 if self.fdm.run.launch_gui: plt.show() 

522 

523 def animate_lossMap_crossSection(self): 

524 """ 

525 This function is similar to the plot_lossMap_crossSection function, but instead of plotting the loss for at a constant crossection,  

526 it sweeps the crossection over a chosen axis and plots the loss for each crossection as an animation. 

527 """ 

528 lm = self.fdm.magnet.postproc.batch_postproc.loss_map 

529 axis = lm.cross_section_sweep.axis_to_sweep 

530 

531 X, Y, V, dataPoints = self.lossMap_createGridData('TotalLoss') 

532 x_vals = X[:, 0] # x-values from the loss map 

533 y_vals = Y[0, :] # y-values from the loss map 

534 

535 if axis == 'x': 

536 A = np.zeros((lm.y_steps, 4, lm.x_steps)) 

537 axis_to_sweep = x_vals 

538 constant_axis = y_vals 

539 elif axis == 'y': 

540 A = np.zeros((lm.x_steps, 4, lm.y_steps)) 

541 axis_to_sweep = y_vals 

542 constant_axis = x_vals 

543 

544 

545 for i, val in enumerate(axis_to_sweep): 

546 _, totalLoss, filamentLoss, eddyLoss, couplingLoss = self.lossMap_crossSection(val, axis_to_cut = axis) 

547 A[:, 0, i] = totalLoss 

548 A[:, 1, i] = filamentLoss 

549 A[:, 2, i] = eddyLoss 

550 A[:, 3, i] = couplingLoss 

551 

552 # Initialize the plot 

553 fig, ax = plt.subplots() 

554 lines = ax.plot(constant_axis, A[:, :, 0], lw=2, label = ['total Loss', 'filament Loss', 'eddy Loss', 'coupling Loss']) 

555 

556 # Set plot limits and labels 

557 ax.set_xlim(constant_axis[0], constant_axis[-1]) 

558 ax.set_ylim(np.min(A), np.max(A)) 

559 ax.set_xlabel(lm.ylabel if axis == 'x' else lm.xlabel) 

560 ax.set_ylabel(lm.cross_section_sweep.ylabel) 

561 

562 

563 # Define the animation update function 

564 def update(frame): 

565 for i, line in enumerate(lines): 

566 line.set_ydata(A[:, i, frame]) 

567 

568 if axis == 'x': 

569 if lm.x_log: 

570 sweep_value = 10**x_vals[frame] 

571 else: 

572 sweep_value = x_vals[frame] 

573 elif axis == 'y': 

574 if lm.y_log: 

575 sweep_value = 10**y_vals[frame] 

576 else: 

577 sweep_value = y_vals[frame] 

578 

579 title = lm.cross_section_sweep.title.replace('<<sweep_value>>', str(round(sweep_value, 3))) 

580 ax.set_title(title) 

581 return lines,# line1, line2, line3, line4 

582 

583 # Create the animation 

584 dt = 0.1 

585 ani = FuncAnimation(fig, update, frames=lm.x_steps if axis == 'x' else lm.y_steps, interval=dt*1000, blit=False) 

586 

587 # Show the animation 

588 plt.legend() 

589 plt.grid() 

590 

591 if lm.cross_section_sweep.save_plot: 

592 filepath = create_non_overwriting_filepath(folder_path=self.outputs_folder_path, base_name=lm.cross_section_sweep.filename, extension='.gif', overwrite=self.fdm.run.overwrite) 

593 ani.save(filepath, writer='imagemagick', fps=1/dt) 

594 

595 if self.fdm.run.launch_gui: plt.show() 

596 

597 def create_lossMap(self): 

598 """ 

599 This function creates a loss map based on the inputs given in the loss_map section of the input file. 

600 The loss-map can be plotted and saved as a .png file. 

601 """ 

602 lm = self.fdm.magnet.postproc.batch_postproc.loss_map 

603 

604 if self.simulation_collection: 

605 X, Y, V, dataPoints = self.lossMap_createGridData(lm.loss_type) 

606 else: 

607 X, Y, V = self.totalLoss_gridData 

608 

609 if is_latex_installed(): 

610 plt.rcParams['text.usetex'] = True 

611 plt.rcParams['font.family'] = 'times' 

612 plt.rcParams['font.size'] = 20 

613 

614 fig, ax = plt.subplots(figsize=(10,8)) 

615 

616 c = plt.pcolormesh(X, Y, V, shading='gouraud', cmap='plasma_r') 

617 c_min = min([np.ceil(np.min(V)) for V in [V]]) 

618 c_max = max([np.floor(np.max(V)) for V in [V]]) 

619 c_ticks = [int(val) for val in np.arange(c_min, c_max+1)] 

620 cont = plt.contour(X,Y,V, c_ticks, colors='k', linestyles='dashed') 

621 

622 if lm.show_datapoints: 

623 plt.scatter(dataPoints[:, 0], dataPoints[:, 1], s=50, edgecolors='k') 

624 

625 

626 if lm.show_loss_type_dominance_contour: 

627 sigmoid = lambda x: 1/(1+np.exp(-x)) 

628 if self.simulation_collection: 

629 _, _, FilamentLoss, _ = self.lossMap_createGridData(lossType='FilamentLoss') 

630 _, _, CouplingLoss, _ = self.lossMap_createGridData(lossType='CouplingLoss') 

631 _, _, EddyLoss, _ = self.lossMap_createGridData(lossType='EddyLoss') 

632 # else: 

633 # _, _, FilamentLoss = self.filamentLoss_gridData 

634 # _, _, CouplingLoss = self.couplingLoss_gridData 

635 # _, _, EddyLoss = self.eddyLoss_gridData 

636 fil_vs_coupling_loss = np.maximum(FilamentLoss, EddyLoss) - CouplingLoss 

637 fil_vs_eddy_loss = EddyLoss - np.maximum(FilamentLoss, CouplingLoss) 

638 plt.contour(X,Y,sigmoid(fil_vs_coupling_loss),[0.5], colors='k') 

639 plt.contour(X,Y,sigmoid(fil_vs_eddy_loss),[0.5], colors='k') 

640 

641 cbar = fig.colorbar(c, ticks=c_ticks)#, labels=c_labels) 

642 # cbar.ax.set_xticks([-7, -6, -5, -4, -3, -2, -1, 0, 1]) 

643 # cbar.ax.set_yticklabels([r'$10^{-7}$', r'$10^{-6}$', r'$10^{-5}$', r'$10^{-4}$', r'$10^{-3}$', r'$10^{-2}$', r'$10^{-1}$', r'$10^0$', r'$10^1$']) 

644 cbar.ax.set_yticklabels([f"$10^{{{val}}}$" for val in c_ticks]) 

645 # plt.grid(alpha=0.5) 

646 # plt.title(lm.title) 

647 # plt.xlabel(lm.xlabel) 

648 # plt.ylabel(lm.ylabel) 

649 plt.title(r'Loss per cycle (J/m)') 

650 plt.xlabel(r'Frequency $f$ (Hz)') 

651 plt.ylabel(r'Field amplitude $b$ (T)') 

652 

653 

654 # plt.annotate(r'Coupling', (np.log10(1.0), np.log10(0.007)), color='white') 

655 # plt.annotate(r'Filament', (np.log10(0.012), np.log10(0.74)), color='white') 

656 # plt.annotate(r'(uncoupled)', (np.log10(0.012), np.log10(0.55)), color='white') 

657 # plt.annotate(r'Filament', (np.log10(45), np.log10(0.38)), color='white') 

658 # plt.annotate(r'(coupled)', (np.log10(45), np.log10(0.28)), color='white') 

659 # plt.annotate(r'Eddy', (np.log10(2000), np.log10(0.03)), color='white') 

660 

661 # ax.plot(np.log10(0.03), np.log10(0.2), 'o', color='white')#, xytext=(np.log10(0.03), np.log10(0.12)), arrowprops=dict(facecolor='black', shrink=0.02)) 

662 # ax.plot(np.log10(30), np.log10(1), 'o', color='white')#, xytext=(np.log10(40), np.log10(0.8)), arrowprops=dict(facecolor='black', shrink=0.02)) 

663 # ax.plot(np.log10(3), np.log10(0.2), 'o', color='white')#, xytext=(np.log10(2), np.log10(0.2)), arrowprops=dict(facecolor='black', shrink=0.02))