Coverage for fiqus/mains/MainMultipole.py: 59%

1import logging 

2import os 

3import gmsh 

4import time 

5 

6from fiqus.plotters.PlotPythonMultipole import PlotPythonMultipole 

7from fiqus.utils.Utils import GmshUtils 

8from fiqus.utils.Utils import FilesAndFolders as Util 

9from fiqus.data import DataFiQuS as dF 

10from fiqus.data.DataRoxieParser import FiQuSGeometry 

11from fiqus.geom_generators.GeometryMultipole import Geometry 

12from fiqus.mesh_generators.MeshMultipole import Mesh 

13from fiqus.getdp_runners.RunGetdpMultipole import RunGetdpMultipole 

14from fiqus.getdp_runners.RunGetdpMultipole import AssignNaming 

15from fiqus.post_processors.PostProcessMultipole import PostProcess 

16 

17 

18class MainMultipole: 

19 def __init__(self, fdm: dF.FDM = None, rgd_path: str = None, verbose: bool = None, inputs_folder_path = None): 

20 """ 

21 Main class for working with simulations for multipole type magnets 

22 :param fdm: FiQuS data model 

23 :param rgd_path: ROXIE geometry data path 

24 :param verbose: if True, more info is printed in the console 

25 """ 

26 self.fdm = fdm 

27 self.rgd = rgd_path 

28 self.verbose = verbose 

29 

30 self.GetDP_path = None 

31 self.geom_folder = None 

32 self.mesh_folder = None 

33 self.solution_folder = None 

34 self.inputs_folder_path = inputs_folder_path 

35 

36 def force_symmetry(self): 

37 fdm = self.fdm.__deepcopy__() 

38 fdm.magnet.geometry.electromagnetics.symmetry = 'x' 

39 return fdm 

40 

41 def generate_geometry(self, gui: bool = False): 

42 geom = Util.read_data_from_yaml(self.rgd, FiQuSGeometry) 

43 fdm = self.force_symmetry() if 'solenoid' in geom.Roxie_Data.coil.coils[1].type else self.fdm # todo: this should be handled by pydantic 

44 if self.fdm.magnet.geometry.plot_preview: 

45 plotter = PlotPythonMultipole(geom, self.fdm) 

46 plotter.plot_coil_wedges() 

47 gg = Geometry(data=fdm, geom=geom, geom_folder=self.geom_folder, verbose=self.verbose) 

48 gg.saveHalfTurnCornerPositions() 

49 

50 geometry_settings = {'EM': fdm.magnet.geometry.electromagnetics, 'TH': self.fdm.magnet.geometry.thermal} 

51 geometry_type_list = [] 

52 

53 if geometry_settings['EM'].create: geometry_type_list.append('EM') 

54 if geometry_settings['TH'].create: geometry_type_list.append('TH') 

55 for geometry_type in geometry_type_list: 

56 gg.saveStrandPositions(geometry_type) 

57 if any(geometry_settings[geometry_type].areas): 

58 gg.constructIronGeometry(geometry_settings[geometry_type].symmetry if geometry_type == 'EM' else 'none', geometry_settings[geometry_type], run_type = geometry_type) 

59 gg.constructCoilGeometry(geometry_type) 

60 if geometry_settings[geometry_type].with_wedges: 

61 gg.constructWedgeGeometry(geometry_settings[geometry_type].use_TSA if geometry_type == 'TH' else False) 

62 gmsh.model.occ.synchronize() 

63 if geometry_type == 'TH': 

64 if geometry_settings[geometry_type].use_TSA: 

65 gg.constructThinShells(geometry_settings[geometry_type].with_wedges) 

66 else: 

67 gg.constructInsulationGeometry() 

68 if geometry_settings[geometry_type].use_TSA_new: 

69 gg.constructAdditionalThinShells() 

70 

71 gg.buildDomains(geometry_type, geometry_settings[geometry_type].symmetry if geometry_type == 'EM' else 'none') 

72 if geometry_type == 'EM': 

73 gg.fragment() 

74 if geometry_type =='TH' and 'poles' in geometry_settings[geometry_type].areas: 

75 # make sure geometry is connected 

76 gmsh.model.occ.removeAllDuplicates() 

77 gmsh.model.occ.synchronize() 

78 

79 gg.saveBoundaryRepresentationFile(geometry_type) 

80 gg.loadBoundaryRepresentationFile(geometry_type) 

81 gg.updateTags(geometry_type, geometry_settings[geometry_type].symmetry if geometry_type == 'EM' else 'none') 

82 gg.saveAuxiliaryFile(geometry_type) 

83 gg.clear() 

84 gg.ending_step(gui) 

85 

86 def load_geometry(self, gui: bool = False): 

87 pass 

88 # gu = GmshUtils(self.geom_folder, self.verbose) 

89 # gu.initialize(verbosity_Gmsh=self.fdm.run.verbosity_Gmsh) 

90 # model_file = os.path.join(self.geom_folder, self.fdm.general.magnet_name) 

91 # gmsh.option.setString(name='Geometry.OCCTargetUnit', value='M') # set units to meters 

92 # gmsh.open(model_file + '_EM.brep') 

93 # gmsh.open(model_file + '_TH.brep') 

94 # if gui: gu.launch_interactive_GUI() 

95 

96 def pre_process(self, gui: bool = False): 

97 pass 

98 

99 def load_geometry_for_mesh(self, run_type): 

100 gu = GmshUtils(self.geom_folder, self.verbose) 

101 gu.initialize(verbosity_Gmsh=self.fdm.run.verbosity_Gmsh) 

102 model_file = os.path.join(self.geom_folder, self.fdm.general.magnet_name) 

103 gmsh.option.setString(name='Geometry.OCCTargetUnit', value='M') # set units to meters 

104 gmsh.open(model_file + f'_{run_type}.brep') 

105 

106 def mesh(self, gui: bool = False): 

107 def _create_physical_group_for_reference(self): 

108 """ 

109 This code generates the reference models for the FALCOND_C 

110 """ 

111 ### hardcoded, we need only one reference 

112 if self.fdm.general.magnet_name == 'TEST_MULTIPOLE_FALCOND_C_TSA_COLLAR_POLE': 

113 CUT_REFERENCE = True # self specify loop and surf 

114 

115 L_col = [56, 57, 58, 35, 16, 49, 48, 47] 

116 l1 = gmsh.model.occ.addLine(746,48) 

117 l2 = gmsh.model.occ.addLine(41,691) 

118 L_inner = list(range(854, 754-1, -1)) 

119 loop1 = gmsh.model.occ.addCurveLoop([l1] + L_col + [l2] + L_inner) 

120 

121 R_col = [52, 53, 54, 26, 6, 45, 44, 43] 

122 r1 = gmsh.model.occ.addLine(38,902) 

123 r2 = gmsh.model.occ.addLine(847,45) 

124 R_inner = list(range(910, 1010+1)) 

125 loop2 = gmsh.model.occ.addCurveLoop([r2]+ R_col + [r1] + R_inner) 

126 

127 surf1 = gmsh.model.occ.addPlaneSurface([loop1]) 

128 surf2 = gmsh.model.occ.addPlaneSurface([loop2]) 

129 surf = [surf1, surf2] # tag 

130 

131 else: raise Exception("Reference meshing is not implemented for this magnet.") 

132 gmsh.model.occ.synchronize() 

133 

134 #add to physical groups / domains -> write to aux file 

135 file = os.path.join(self.geom_folder, self.fdm.general.magnet_name + '_TH.aux') 

136 with open(file) as f: 

137 lines = f.readlines() 

138 updated_lines = [] 

139 flag = 0 

140 for line in lines: 

141 updated_lines.append(line) 

142 if flag == 0 and 'groups_entities:' in line: 

143 flag = 1 

144 if flag==1 and ('ref_mesh: {}' in line): 

145 if type(surf) == list: 

146 updated_lines[-1] = f' ref_mesh:\n {self.fdm.magnet.solve.thermal.insulation_TSA.between_collar.material}: {surf}\n' 

147 else: 

148 updated_lines[-1] = f' ref_mesh:\n {self.fdm.magnet.solve.thermal.insulation_TSA.between_collar.material}: [{surf}]\n' 

149 #if not CUT_REFERENCE else f' ref_mesh:\n {self.fdm.magnet.solve.thermal.insulation_TSA.between_collar.material}: '+ str(surf) +'\n' 

150 flag = -1 

151 

152 with open(file, 'w') as f: 

153 f.writelines(updated_lines) 

154 logger = logging.getLogger('FiQuS') 

155 logger.warning("Overwrite the .aux file to include the new domain") 

156 

157 mm = Mesh(data=self.fdm, mesh_folder=self.mesh_folder, verbose=self.verbose) ## same mesh object is used for both thermal and EM 

158 geom = Util.read_data_from_yaml(self.rgd, FiQuSGeometry) 

159 fdm = self.force_symmetry() if 'solenoid' in geom.Roxie_Data.coil.coils[1].type else self.fdm 

160 geometry_settings = {'EM': fdm.magnet.geometry.electromagnetics, 'TH': self.fdm.magnet.geometry.thermal} 

161 mesh_settings = {'EM': fdm.magnet.mesh.electromagnetics, 'TH': fdm.magnet.mesh.thermal} 

162 

163 mesh_type_list = [] 

164 if mesh_settings['EM'].create: mesh_type_list.append('EM') 

165 if mesh_settings['TH'].create: mesh_type_list.append('TH') 

166 for physics_solved in mesh_type_list: 

167 self.load_geometry_for_mesh(physics_solved) 

168 if physics_solved == 'TH' and mesh_settings['TH'].reference.enabled and 'collar' in geometry_settings['TH'].areas: 

169 if self.fdm.magnet.geometry.thermal.use_TSA_new or self.fdm.magnet.geometry.thermal.use_TSA: 

170 raise Exception('Reference solution is not implemented for collar with TSA') 

171 if 'iron_yoke' in geometry_settings['TH'].areas: 

172 raise Exception('Reference solution is intended (read: hardcoded) without iron yoke') 

173 _create_physical_group_for_reference(self) 

174 mm.loadAuxiliaryFile(physics_solved) 

175 if any(geometry_settings[physics_solved].areas): 

176 mm.getIronCurvesTags(physics_solved) 

177 mm.defineMesh(geometry_settings[physics_solved], mesh_settings[physics_solved], physics_solved) 

178 mm.createPhysicalGroups(geometry_settings[physics_solved]) 

179 mm.updateAuxiliaryFile(physics_solved) 

180 if geometry_settings[physics_solved].model_dump().get('use_TSA', False): 

181 mm.rearrangeThinShellsData() # rearrange data for the pro file, technically optional 

182 mm.assignRegionsTags(geometry_settings[physics_solved], mesh_settings[physics_solved]) 

183 mm.saveRegionFile(physics_solved) 

184 mm.setMeshOptions() 

185 mm.generateMesh() 

186 mm.checkMeshQuality() 

187 mm.saveMeshFile(physics_solved) 

188 if geometry_settings[physics_solved].model_dump().get('use_TSA', False): 

189 mm.saveClosestNeighboursList() 

190 if self.fdm.magnet.mesh.thermal.isothermal_conductors: mm.selectMeshNodes(elements='conductors') 

191 if self.fdm.magnet.geometry.thermal.with_wedges and self.fdm.magnet.mesh.thermal.isothermal_wedges: mm.selectMeshNodes(elements='wedges') 

192 if geometry_settings[physics_solved].model_dump().get('use_TSA_new', False): 

193 mm.saveClosestNeighboursList_new_TSA() 

194 mm.saveHalfTurnCornerPositions() 

195 mm.saveRegionCoordinateFile(physics_solved) 

196 mm.clear() 

197 mm.ending_step(gui) 

198 return mm.mesh_parameters 

199 

200 def load_mesh(self, gui: bool = False): 

201 gu = GmshUtils(self.geom_folder, self.verbose) 

202 gu.initialize(verbosity_Gmsh=self.fdm.run.verbosity_Gmsh) 

203 gmsh.open(f"{os.path.join(self.mesh_folder, self.fdm.general.magnet_name)}.msh") 

204 if gui: gu.launch_interactive_GUI() 

205 

206 def solve_and_postprocess_getdp(self, gui: bool = False): 

207 an = AssignNaming(data=self.fdm) 

208 rg = RunGetdpMultipole(data=an, solution_folder=self.solution_folder, GetDP_path=self.GetDP_path, 

209 verbose=self.verbose) 

210 rg.loadRegionFiles() 

211 if self.fdm.magnet.solve.thermal.solve_type and self.fdm.magnet.geometry.thermal.use_TSA: 

212 rg.loadRegionCoordinateFile() 

213 rg.assemblePro() 

214 start_time = time.time() 

215 rg.solve_and_postprocess() 

216 rg.ending_step(gui) 

217 return time.time() - start_time 

218 

219 def solve_and_postprocess_getdp(self, gui: bool = False): 

220 an = AssignNaming(data=self.fdm) 

221 rg = RunGetdpMultipole(data=an, solution_folder=self.solution_folder, GetDP_path=self.GetDP_path, 

222 verbose=self.verbose) 

223 rg.loadRegionFiles() 

224 if self.fdm.magnet.solve.thermal.solve_type and self.fdm.magnet.geometry.thermal.use_TSA: 

225 rg.loadRegionCoordinateFile() 

226 rg.read_aux_file(os.path.join(self.mesh_folder, f"{self.fdm.general.magnet_name}_EM.aux")) 

227 rg.extract_half_turn_blocks() 

228 if self.fdm.magnet.geometry.thermal.use_TSA_new: 

229 rg.read_aux_file(os.path.join(self.mesh_folder, 

230 f"{self.fdm.general.magnet_name}_TH.aux")) # now load the thermal aux file 

231 rg.extract_specific_TSA_lines() 

232 rg.assemblePro() 

233 start_time = time.time() 

234 rg.solve_and_postprocess() 

235 rg.ending_step(gui) 

236 return time.time() - start_time 

237 

238 

239 def post_process_getdp(self, gui: bool = False): 

240 an = AssignNaming(data=self.fdm) 

241 rg = RunGetdpMultipole(data=an, solution_folder=self.solution_folder, GetDP_path=self.GetDP_path, verbose=self.verbose) 

242 rg.loadRegionFiles() 

243 if self.fdm.magnet.solve.thermal.solve_type and self.fdm.magnet.geometry.thermal.use_TSA: 

244 rg.loadRegionCoordinateFile() 

245 rg.assemblePro() 

246 rg.postprocess() 

247 rg.ending_step(gui) 

248 

249 def post_process_python(self, gui: bool = False): 

250 if self.fdm.run.type == 'post_process_python_only': 

251 an = AssignNaming(data=self.fdm) 

252 data = an.data 

253 else: data = self.fdm 

254 

255 run_types = [] 

256 if self.fdm.magnet.solve.electromagnetics.solve_type: run_types.append('EM') 

257 if self.fdm.magnet.solve.thermal.solve_type: run_types.append('TH') 

258 pp_settings = {'EM': self.fdm.magnet.postproc.electromagnetics, 'TH': self.fdm.magnet.postproc.thermal} 

259 pp = PostProcess(data=data, solution_folder=self.solution_folder, verbose=self.verbose) 

260 for run_type in run_types: 

261 pp.prepare_settings(pp_settings[run_type]) 

262 pp.loadStrandPositions(run_type) 

263 pp.loadAuxiliaryFile(run_type) 

264 if pp_settings[run_type].plot_all != 'False': pp.loadHalfTurnCornerPositions() 

265 if pp_settings[run_type].model_dump().get('take_average_conductor_temperature', False): pp.loadRegionFile() 

266 pp.postProcess(pp_settings[run_type]) 

267 if run_type == 'EM' and self.fdm.magnet.geometry.electromagnetics.symmetry != 'none': pp.completeMap2d() 

268 pp.clear() 

269 pp.ending_step(gui) 

270 return pp.postprocess_parameters 

271 

272 def plot_python(self): 

273 os.chdir(self.solution_folder) 

274 p = PlotPythonMultipole(self.fdm, self.fdm) 

275 p.plot_coil_wedges()