Coverage for fiqus/post_processors/PostProcessMultipole.py: 97%
198 statements
« prev ^ index » next coverage.py v6.4.4, created at 2024-05-20 03:24 +0200
1import os
2from pathlib import Path
3import gmsh
4import json
5import numpy as np
6import matplotlib.pyplot as plt
7import matplotlib.lines as lines
9from fiqus.utils.Utils import GmshUtils
10from fiqus.utils.Utils import GeometricFunctions as Func
11from fiqus.utils.Utils import RoxieParsers as Pars
12from fiqus.data import DataFiQuS as dF
15class PostProcess:
16 def __init__(self, data: dF.FDM() = None, sett: dF.FiQuSSettings() = None, solution_folder: str = None,
17 verbose: bool = False):
18 """
19 Class to post process results
20 :param data: FiQuS data model
21 :param sett: settings data model
22 :param verbose: If True more information is printed in python console.
23 """
24 self.data: dF.FDM() = data
25 self.set: dF.FiQuSSettings() = sett
26 self.solution_folder = solution_folder
27 self.verbose: bool = verbose
29 self.gu = GmshUtils(self.solution_folder, self.verbose)
30 self.gu.initialize()
32 self.brep_iron_curves = {1: set(), 2: set(), 3: set(), 4: set()}
33 self.strands = None
34 self.crns = None
35 self.BB_err_mean = []
36 self.BB_err_min = []
37 self.BB_err_max = []
38 self.postprocess_parameters = dict.fromkeys(['overall_error', 'minimum_diff', 'maximum_diff'])
39 self.geom_folder = os.path.dirname(os.path.dirname(self.solution_folder))
40 self.model_file = f"{os.path.join(self.solution_folder, self.data.general.magnet_name)}.map2d"
42 self.II = (self.set.Model_Data_GS.general_parameters.I_ref[0] if self.data.magnet.postproc.compare_to_ROXIE
43 else self.data.magnet.solve.I_initial[0])
45 if self.data.magnet.postproc.plot_all != 'False':
46 self.fiqus = None
47 self.roxie = None
48 plt.figure(1)
49 self.ax = plt.axes()
50 self.ax.set_xlabel('x [m]')
51 self.ax.set_ylabel('y [m]')
52 self.ax.set_xlim(0, 0.09)
53 self.ax.set_ylim(0, 0.04)
55 if self.data.magnet.postproc.compare_to_ROXIE:
56 fig2 = plt.figure(2)
57 self.ax2 = fig2.add_subplot(projection='3d')
58 self.ax2.set_xlabel('x [m]')
59 self.ax2.set_ylabel('y [m]')
60 self.ax2.set_zlabel('Absolute Error [T]')
61 self.fig4 = plt.figure(4)
62 self.ax4 = plt.axes()
63 self.ax4.set_xlabel('x [cm]')
64 self.ax4.set_ylabel('y [cm]')
65 self.ax4.set_aspect('equal', 'box')
66 fig3 = plt.figure(3)
67 self.ax3 = fig3.add_subplot(projection='3d')
68 self.ax3.set_xlabel('x [m]')
69 self.ax3.set_ylabel('y [m]')
70 self.ax3.set_zlabel('norm(B) [T]')
72 name = 'b_Omega_p'
73 gmsh.open(os.path.join(self.solution_folder, f"{name}.pos"))
75 def ending_step(self, gui: bool = False):
76 if gui:
77 self.gu.launch_interactive_GUI()
78 else:
79 gmsh.clear()
80 gmsh.finalize()
82 def loadStrandPositions(self):
83 self.strands = json.load(open(f"{os.path.join(self.geom_folder, self.data.general.magnet_name)}.strs"))
85 def loadHalfTurnCornerPositions(self):
86 self.crns = json.load(open(f"{os.path.join(self.geom_folder, self.data.general.magnet_name)}.crns"))
88 def plotHalfTurnGeometry(self):
89 for i in range(len(self.crns['iL'])):
90 self.ax.add_line(lines.Line2D([self.crns['iL'][i][0], self.crns['iR'][i][0]],
91 [self.crns['iL'][i][1], self.crns['iR'][i][1]], color='green'))
92 self.ax.add_line(lines.Line2D([self.crns['oL'][i][0], self.crns['oR'][i][0]],
93 [self.crns['oL'][i][1], self.crns['oR'][i][1]], color='green'))
94 self.ax.add_line(lines.Line2D([self.crns['oR'][i][0], self.crns['iR'][i][0]],
95 [self.crns['oR'][i][1], self.crns['iR'][i][1]], color='green'))
96 self.ax.add_line(lines.Line2D([self.crns['iL'][i][0], self.crns['oL'][i][0]],
97 [self.crns['iL'][i][1], self.crns['oL'][i][1]], color='green'))
98 cc_fiqus = Func.centroid([self.crns['iL'][i][0], self.crns['iR'][i][0],
99 self.crns['oR'][i][0], self.crns['oL'][i][0]],
100 [self.crns['iL'][i][1], self.crns['iR'][i][1],
101 self.crns['oR'][i][1], self.crns['oL'][i][1]])
103 if self.data.magnet.postproc.compare_to_ROXIE:
104 self.ax.add_line(lines.Line2D([self.crns['iLr'][i][0], self.crns['iRr'][i][0]],
105 [self.crns['iLr'][i][1], self.crns['iRr'][i][1]],
106 color='red', linestyle='dashed'))
107 self.ax.add_line(lines.Line2D([self.crns['oLr'][i][0], self.crns['oRr'][i][0]],
108 [self.crns['oLr'][i][1], self.crns['oRr'][i][1]],
109 color='red', linestyle='dashed'))
110 self.ax.add_line(lines.Line2D([self.crns['oRr'][i][0], self.crns['iRr'][i][0]],
111 [self.crns['oRr'][i][1], self.crns['iRr'][i][1]],
112 color='red', linestyle='dashed'))
113 self.ax.add_line(lines.Line2D([self.crns['iLr'][i][0], self.crns['oLr'][i][0]],
114 [self.crns['iLr'][i][1], self.crns['oLr'][i][1]],
115 color='red', linestyle='dashed'))
116 cc_roxie = Func.centroid(
117 [self.crns['iLr'][i][0], self.crns['iRr'][i][0], self.crns['oRr'][i][0], self.crns['oLr'][i][0]],
118 [self.crns['iLr'][i][1], self.crns['iRr'][i][1], self.crns['oRr'][i][1], self.crns['oLr'][i][1]])
120 self.fiqus = self.ax.scatter(cc_fiqus[0], cc_fiqus[1], c="green")
121 if self.data.magnet.postproc.compare_to_ROXIE:
122 self.roxie = self.ax.scatter(cc_roxie[0], cc_roxie[1], edgecolor='r', facecolor='none')
124 def postProcess(self):
125 def _printState(bb):
126 perc = round(bb / blocks_nr * 100)
127 print("Info : [" + f"{' ' if perc < 10 else ''}" + f"{' ' if perc < 100 else ''}" + f"{perc}" +
128 "%] Interpolating within block" + f"{str(bb)}")
130 def _fetchRoxieData():
131 roxie_x.append(matrix[row, 3])
132 roxie_y.append(matrix[row, 4])
133 strands_area.append(matrix[row, 7])
134 BB_roxie_x.append(matrix[row, 5])
135 BB_roxie_y.append(matrix[row, 6])
136 BB_roxie.append(np.linalg.norm(np.array([matrix[row, 5], matrix[row, 6]])))
138 def _probeFromView():
139 BB_x = []
140 BB_y = []
141 BB = []
142 b_list = [b for b, field in enumerate(self.data.magnet.postproc.variables) if field == 'b']
143 if not b_list:
144 raise Exception("The interpolation of the field at the strands locations can not be executed: "
145 "the field 'b' is not listed in 'post_processors' -> 'variables'")
146 b_field_index = 0
147 while self.data.magnet.postproc.volumes[b_list[b_field_index]] != 'Omega_p':
148 b_field_index += 1
149 # self.data.magnet.post_proc.variables.index('b')
150 for ss in range(len(strands_x)):
151 probe_data = gmsh.view.probe(gmsh.view.getTags()[0] if len(b_list) == 1 else b_list[b_field_index],
152 strands_x[ss], strands_y[ss], 0)[0]
153 BB_x.append(probe_data[0])
154 BB_y.append(probe_data[1])
155 BB.append(np.linalg.norm(np.array([BB_x[-1], BB_y[-1]])))
156 return BB_x, BB_y, BB
158 def _updateMap2dFile(blk_nr, ht_nr):
159 with open(self.model_file, 'a') as f:
160 content = []
161 s = row - strands
162 for x, y, Bx, By in zip(strands_x, strands_y, BB_strands_x, BB_strands_y):
163 s += 1
164 content.append(f" {blk_nr} {ht_nr} {s} {x * 1e3:.4f} "
165 f"{y * 1e3:.4f} {Bx:.4f} {By:.4f} {1.0} "
166 f"{self.II / strands:.2f} {1.0}\n")
167 f.writelines(content)
169 def _computeError():
170 BB_err = BB_strands - np.array(BB_roxie)
171 self.BB_err_mean.append(np.mean(abs(BB_err)))
172 self.BB_err_min.append(np.min(abs(BB_err)))
173 self.BB_err_min.append(np.max(abs(BB_err)))
174 return BB_err
176 def _plotData():
177 map2d = None
178 pos_roxie = None
179 if self.data.magnet.postproc.compare_to_ROXIE:
180 map2d = self.ax.scatter(roxie_x, roxie_y, edgecolor='black', facecolor='black', s=10)
181 corners = conductorPositionsList[c + cond_nr].xyCorner
182 for corner in range(len(corners)):
183 self.ax.scatter(corners[corner][0] / 1e3, corners[corner][1] / 1e3, edgecolor='black',
184 facecolor='black', s=10)
185 if flag_contraction:
186 self.ax.scatter([x * (1 - 0.002) for x in parser_x],
187 [y * (1 - 0.002) for y in parser_y], c='r', s=10)
188 else:
189 self.ax.scatter(parser_x, parser_y, c='r', s=10)
190 self.ax2.scatter3D(roxie_x, roxie_y, BB_abs_err, c=BB_abs_err, cmap='viridis') # , vmin=-0.2, vmax=0.2)
191 self.ax4.scatter(np.array(roxie_x) * 1e2, np.array(roxie_y) * 1e2, s=1, c=np.array(BB_abs_err) * 1e3, cmap='viridis')
192 pos_roxie = self.ax3.scatter3D(roxie_x, roxie_y, BB_roxie, c=BB_roxie, cmap='Reds', vmin=0, vmax=10)
193 pos = self.ax3.scatter3D(strands_x, strands_y, BB_strands, c=BB_strands, cmap='Greens', vmin=0, vmax=10)
194 return map2d, pos_roxie, pos
196 print(f"Info : {self.data.general.magnet_name} - I n t e r p o l a t i n g . . .")
197 print("Info : Interpolating magnetic flux density ...")
199 if self.data.magnet.postproc.compare_to_ROXIE:
200 roxie_path = self.data.magnet.postproc.compare_to_ROXIE # f"C:\\Users\\avitrano\\PycharmProjects\\steam_sdk\\tests\\builders\\model_library\\magnets\\{self.data.general.magnet_name}\\input"
201 flag_contraction = False
202 # flag_self_field = False
203 path_map2d = Path(roxie_path)
204 # path_map2d = Path(roxie_path, "MQXA_All_" +
205 # f"{'WithIron_' if self.data.magnet.geometry.with_iron_yoke else 'NoIron_'}" +
206 # f"{'WithSelfField' if flag_self_field else 'NoSelfField'}" +
207 # f"{'' if flag_contraction else '_no_contraction'}" + ".map2d")
208 matrix = Pars.parseMap2d(map2dFile=path_map2d)
210 if self.data.magnet.postproc.plot_all != 'False':
211 path_cond2d = Path(os.path.join(os.path.dirname(roxie_path), self.data.general.magnet_name + ".cond2d"))
212 # path_cond2d = Path(os.path.dirname(roxie_path), "MQXA_All_NoIron_NoSelfField" +
213 # f"{'' if flag_contraction else '_no_contraction'}" + ".cond2d")
214 conductorPositionsList = Pars.parseCond2d(path_cond2d)
216 with open(self.model_file, 'w') as file:
217 file.write(" BL. COND. NO. X-POS/MM Y-POS/MM BX/T BY/T"
218 " AREA/MM**2 CURRENT FILL FAC.\n\n")
220 cond_nr = 0
221 row = 0
222 c_nr_list = [] # list of conductors per block
223 blocks_nr = self.strands['block'][-1]
224 blocks, hts = np.array(self.strands['block']), np.array(self.strands['ht'])
225 for blk in range(blocks_nr):
226 _printState(blk + 1)
227 c_nr_list.append(hts[np.where(blocks == blk + 1)[0][-1]] - sum(c_nr_list))
228 for c in range(c_nr_list[blk]):
229 parser_x, parser_y, roxie_x, roxie_y, strands_area, BB_roxie, BB_roxie_x, BB_roxie_y = \
230 [], [], [], [], [], [], [], []
231 strands = 0
232 ht = hts[row]
233 while ht == c + 1 + cond_nr and row < len(self.strands['x']):
234 parser_x.append(self.strands['x'][row])
235 parser_y.append(self.strands['y'][row])
236 if self.data.magnet.postproc.compare_to_ROXIE:
237 _fetchRoxieData()
238 strands += 1
239 row += 1
240 if row < len(self.strands['x']) - 1:
241 ht = hts[row]
243 strands_x = roxie_x if self.data.magnet.postproc.compare_to_ROXIE else parser_x
244 strands_y = roxie_y if self.data.magnet.postproc.compare_to_ROXIE else parser_y
246 BB_strands_x, BB_strands_y, BB_strands = _probeFromView()
248 _updateMap2dFile(blk + 1, c + 1 + cond_nr)
250 if self.data.magnet.postproc.compare_to_ROXIE:
251 BB_abs_err = _computeError()
253 if self.data.magnet.postproc.plot_all != 'False': # and (blk == 0 or blk == 1):
254 if self.data.magnet.postproc.compare_to_ROXIE:
255 map2d_strands, scatter3D_pos_roxie, scatter3D_pos = _plotData()
256 else:
257 scatter3D_pos = _plotData()[2]
259 cond_nr += c_nr_list[blk]
261 if self.data.magnet.postproc.compare_to_ROXIE:
262 self.postprocess_parameters['overall_error'] = np.linalg.norm(self.BB_err_mean)
263 self.postprocess_parameters['minimum_diff'] = np.min(self.BB_err_mean)
264 self.postprocess_parameters['maximum_diff'] = np.max(self.BB_err_mean)
266 print(f"Info : {self.data.general.magnet_name} - E n d I n t e r p o l a t i n g")
268 if self.data.magnet.postproc.plot_all != 'False' and self.data.magnet.postproc.plot_all != False:
269 self.plotHalfTurnGeometry()
270 if self.data.magnet.postproc.compare_to_ROXIE:
271 cax4 = self.fig4.add_axes([self.ax4.get_position().x1 + 0.02, self.ax4.get_position().y0,
272 0.02, self.ax4.get_position().height])
273 cbar = plt.colorbar(self.ax4.get_children()[2], cax=cax4)
274 cbar.ax.set_ylabel('Absolute error [mT]', rotation=270)
275 self.ax3.legend([scatter3D_pos, scatter3D_pos_roxie], ['FiQuS', 'ROXIE'], numpoints=1)
276 # pickle.dump(fig, open('Bfield.fig.pickle', 'wb'))
277 self.ax.legend([self.fiqus, self.roxie, map2d_strands], ['FiQuS', 'ROXIEparser', 'ROXIE'], numpoints=1)
279 self.fig4.savefig(f"{os.path.join(self.solution_folder, self.data.general.magnet_name)}.svg", bbox_inches='tight')
281 if self.data.magnet.postproc.plot_all == 'True':
282 plt.show()
284 # os.remove(os.path.join(self.solution_folder, 'b_Omega_p.pos'))
285 # os.remove(f"{os.path.join(self.solution_folder, self.data.general.magnet_name)}.pre")
286 # os.remove(f"{os.path.join(os.path.dirname(self.solution_folder), self.data.general.magnet_name)}.msh")