Coverage for fiqus/mains/MainPancake3D.py: 41%
204 statements
« prev ^ index » next coverage.py v6.4.4, created at 2024-05-20 03:24 +0200
1import os
2import sys
3import math
4import time as Time
5import shutil
6from copy import deepcopy
8from fiqus.data.DataFiQuSPancake3D import (
9 Pancake3DGeometry,
10 Pancake3DMesh,
11 Pancake3DSolve,
12 Pancake3DPostprocess,
13 Pancake3DSolveSaveQuantity,
14)
16if len(sys.argv) == 3:
17 sys.path.insert(0, os.path.join(os.getcwd(), "steam-fiqus-dev"))
20class Base:
21 """
22 Base class for geometry, mesh, and solution classes. It is created to avoid code
23 duplication and to make the code more readable. Moreover, it guarantees that
24 all the classes have the same fundamental methods and attributes.
25 """
27 def __init__(
28 self,
29 fdm,
30 geom_folder=None,
31 mesh_folder=None,
32 solution_folder=None,
33 ) -> None:
34 """
35 Define the fundamental attributes of the class.
37 :param fdm: fiqus data model
38 :param geom_folder: folder where the geometry related files are stored
39 :type geom_folder: str
40 :param mesh_folder: folder where the mesh related files are stored
41 :type mesh_folder: str
42 :param solution_folder: folder where the solution related files are stored
43 :type solution_folder: str
44 """
46 self.magnet_name = fdm.general.magnet_name
47 self.geom_folder = geom_folder
48 self.mesh_folder = mesh_folder
49 self.solution_folder = solution_folder
51 self.dm = fdm # Data model
52 self.geo: Pancake3DGeometry = fdm.magnet.geometry # Geometry data
53 self.mesh: Pancake3DMesh = fdm.magnet.mesh # Mesh data
54 self.solve: Pancake3DSolve = fdm.magnet.solve # Solve data
55 self.pp: Pancake3DPostprocess = fdm.magnet.postproc # Postprocess data
57 self.geo_gui = False
58 self.mesh_gui = False
59 self.solve_gui = False
60 self.python_postprocess_gui = False
62 if fdm.run.launch_gui:
63 if fdm.run.type == "start_from_yaml":
64 self.python_postprocess_gui = True
65 elif fdm.run.type == "geometry_only":
66 self.geo_gui = True
67 elif fdm.run.type == "mesh_only":
68 self.mesh_gui = True
69 elif fdm.run.type == "geometry_and_mesh":
70 self.mesh_gui = True
71 elif fdm.run.type == "solve_only":
72 self.solve_gui = True
73 elif fdm.run.type == "post_process_getdp_only":
74 self.solve_gui = True
75 elif fdm.run.type == "solve_with_post_process_python":
76 self.python_postprocess_gui = True
77 elif fdm.run.type == "post_process_python_only":
78 self.python_postprocess_gui = True
80 # Geometry related files:
81 if self.geom_folder is not None:
82 self.brep_file = os.path.join(self.geom_folder, f"{self.magnet_name}.brep")
83 self.vi_file = os.path.join(self.geom_folder, f"{self.magnet_name}.vi")
84 self.geometry_data_file = os.path.join(self.geom_folder, "geometry.yaml")
86 # Mesh related files:
87 if self.mesh_folder is not None:
88 self.mesh_file = os.path.join(self.mesh_folder, f"{self.magnet_name}.msh")
89 self.regions_file = os.path.join(
90 self.mesh_folder, f"{self.magnet_name}.regions"
91 )
92 self.mesh_data_file = os.path.join(self.mesh_folder, "mesh.yaml")
93 # Solution related files:
94 if self.solution_folder is not None:
95 self.pro_file = os.path.join(
96 self.solution_folder, f"{self.magnet_name}.pro"
97 )
100from fiqus.geom_generators.GeometryPancake3D import Geometry
101from fiqus.mesh_generators.MeshPancake3D import Mesh
102from fiqus.getdp_runners.RunGetdpPancake3D import Solve
103from fiqus.post_processors.PostProcessPancake3D import Postprocess
104from fiqus.data.DataFiQuS import FDM
105from fiqus.data.DataFiQuS import PowerSupply
108class MainPancake3D:
109 """
110 The main class for working with simulations for high-temperature superconductor
111 pancake coil magnets.
113 Geometry can be created and saved as a BREP file. Parameters like the number of
114 turns, tape dimensions, contact layer thicknesses, and other dimensions can be
115 specified. Contact layers can be modeled as two-dimensional shells or three-dimensional
116 volumes. Moreover, multiple pancakes can be stacked on top of each other.
118 Using the BREP file created, a mesh can be generated and saved as an MSH file.
119 Winding mesh can be structured, and parameters like, azimuthal number of elements
120 per turn, axial number of elements, and radial number of elements per turn can be
121 specified for each pancake coil. The appropriate regions will be assigned to the
122 relevant volumes accordingly so that finite element simulations can be done.
124 Using the mesh files, GetDP can be used to analyze Pancake3D coils.
126 :param fdm: FiQuS data model
127 """
129 def __init__(self, fdm, verbose):
130 self.fdm: FDM = fdm
131 self.GetDP_path = None
133 self.geom_folder = None
134 self.mesh_folder = None
135 self.solution_folder = None
137 def generate_geometry(self, gui=False):
138 """
139 Generates the geometry of the magnet and save it as a BREP file. Moreover, a
140 text file with the extension VI (volume information file) is generated, which
141 stores the names of the volume tags in JSON format.
142 """
143 geometry = Geometry(
144 self.fdm,
145 self.geom_folder,
146 self.mesh_folder,
147 self.solution_folder,
148 )
150 geometry.generate_geometry()
151 geometry.generate_vi_file()
153 self.model_file = geometry.brep_file
155 def load_geometry(self, gui=False):
156 """
157 Loads the previously generated geometry from the BREP file.
158 """
159 geometry = Geometry(
160 self.fdm,
161 self.geom_folder,
162 self.mesh_folder,
163 self.solution_folder,
164 )
166 geometry.load_geometry()
167 self.model_file = geometry.brep_file
169 def pre_process(self, gui=False):
170 pass
172 def mesh(self, gui=False):
173 """
174 Generates the mesh of the magnet, creates the physical regions, and saves it as
175 an MSH file. Moreover, a text file with the extension REGIONS is generated,
176 which stores the names and tags of the physical regions in YAML format.
177 """
178 mesh = Mesh(
179 self.fdm,
180 self.geom_folder,
181 self.mesh_folder,
182 self.solution_folder,
183 )
185 mesh.generate_mesh()
186 mesh.generate_regions()
187 mesh.generate_mesh_file()
189 self.model_file = mesh.mesh_file
191 return {"gamma": 0} # to be modified with mesh_parameters (see multipole)
193 def load_mesh(self, gui=False):
194 """
195 Loads the previously generated mesh from the MSH file.
196 """
197 mesh = Mesh(
198 self.fdm,
199 self.geom_folder,
200 self.mesh_folder,
201 self.solution_folder,
202 )
203 mesh.load_mesh()
205 self.model_file = mesh.mesh_file
207 def solve_and_postprocess_getdp(self, gui=False):
208 """
209 Simulates the Pancake3D magnet with GetDP using the created mesh file and post
210 processes the results.
211 """
212 # calculate inductance and time constant:
213 if self.fdm.magnet.solve.save is not None:
214 save_list = list(
215 map(
216 lambda quantity_object: quantity_object.quantity,
217 self.fdm.magnet.solve.save,
218 )
219 )
220 new_solve_object_for_inductance_and_time_constant = {
221 "type": "electromagnetic",
222 "time": {
223 "timeSteppingType": "adaptive",
224 "extrapolationOrder": 1,
225 "adaptiveSteppingSettings": {
226 "initialStep": 1,
227 "minimumStep": 0.00001,
228 "maximumStep": 20,
229 "integrationMethod": "Euler",
230 "tolerances": [
231 {
232 "quantity": "magnitudeOfMagneticField",
233 "position": {
234 "x": 0,
235 "y": 0,
236 "z": 0,
237 },
238 "relative": 0.05, # 5%
239 "absolute": 0.00002, # 0.00002 T
240 "normType": "LinfNorm",
241 }
242 ],
243 },
244 },
245 "nonlinearSolver": self.fdm.magnet.solve.nls,
246 "air": {
247 "permeability": 1.2566e-06,
248 },
249 "initialConditions": {
250 "temperature": 4,
251 },
252 }
253 inductance = None
254 inductance_solution_folder = os.path.join(
255 self.solution_folder,
256 "inductance",
257 )
258 inductance_file = os.path.join(
259 inductance_solution_folder, "Inductance-TimeTableFormat.txt"
260 )
261 if "inductance" in save_list:
262 # to calculate inductance and time constant, we should run the simulation
263 # with specigic power supply and material settings.
264 copy_fdm = deepcopy(self.fdm)
266 new_solve_object_for_inductance_and_time_constant["time"]["start"] = 0
267 new_solve_object_for_inductance_and_time_constant["time"]["end"] = 100
268 new_solve_object_for_inductance_and_time_constant["winding"] = {
269 "resistivity": 1e-20,
270 "thermalConductivity": 1,
271 "specificHeatCapacity": 1,
272 }
274 new_solve_object_for_inductance_and_time_constant["contactLayer"] = {
275 "resistivity": 1e-2,
276 "thermalConductivity": 1,
277 "specificHeatCapacity": 1,
278 "numberOfThinShellElements": 1,
279 }
280 new_solve_object_for_inductance_and_time_constant["terminals"] = {
281 "resistivity": 1e-10,
282 "thermalConductivity": 1,
283 "specificHeatCapacity": 1,
284 "terminalContactLayer": {
285 "resistivity": 1e-10,
286 "thermalConductivity": 1,
287 "specificHeatCapacity": 1,
288 },
289 "transitionNotch": {
290 "resistivity": 1e-10,
291 "thermalConductivity": 1,
292 "specificHeatCapacity": 1,
293 },
294 }
295 solve_object = Pancake3DSolve(
296 **new_solve_object_for_inductance_and_time_constant
297 )
298 solve_object.save = [Pancake3DSolveSaveQuantity(quantity="inductance")]
299 copy_fdm.magnet.solve = solve_object
301 new_power_supply_object_for_inductance = {
302 "t_control_LUT": [0, 0.1, 100],
303 "I_control_LUT": [0, 1, 1],
304 }
305 new_power_supply_object_for_inductance = PowerSupply(
306 **new_power_supply_object_for_inductance
307 )
308 copy_fdm.power_supply = new_power_supply_object_for_inductance
310 # generate the folder if it does not exist:
311 if not os.path.exists(inductance_solution_folder):
312 os.makedirs(inductance_solution_folder)
313 else:
314 shutil.rmtree(inductance_solution_folder)
315 os.makedirs(inductance_solution_folder)
317 solve = Solve(
318 copy_fdm,
319 self.GetDP_path,
320 self.geom_folder,
321 self.mesh_folder,
322 inductance_solution_folder,
323 )
324 solve.assemble_pro()
326 start_time = Time.time()
327 solve.run_getdp(solve=True, postOperation=True)
329 with open(inductance_file, "r") as f:
330 # read the last line and second column:
331 inductance = float(f.readlines()[-1].split()[1].replace(",", ""))
333 with open(inductance_file, "w") as f:
334 f.write(str(inductance))
336 if "timeConstant" in save_list:
337 # to calculate inductance and time constant, we should run the simulation
338 # with specigic power supply and material settings.
339 copy_fdm = deepcopy(self.fdm)
341 if inductance is None:
342 raise ValueError(
343 "Time constant can not be calculated without inductance. Please"
344 " add 'inductance' to the 'quantitiesToBeSaved' list in the"
345 " 'solve' section of the YAML file."
346 )
348 # equivalent radial resistance of the winding:
349 # https://doi.org/10.1109/TASC.2021.3063653
350 R_eq = 0
351 for i in range(1, int(copy_fdm.magnet.geometry.wi.N)):
352 r_i = (
353 copy_fdm.magnet.geometry.wi.r_i
354 + i * copy_fdm.magnet.geometry.wi.t
355 + (i - 1) * copy_fdm.magnet.geometry.ii.t
356 + copy_fdm.magnet.geometry.ii.t / 2
357 )
358 w_d = copy_fdm.magnet.geometry.wi.h
359 R_ct = (
360 copy_fdm.magnet.solve.ii.resistivity
361 * copy_fdm.magnet.geometry.ii.t
362 )
363 R_eq = R_eq + R_ct / (2 * math.pi * r_i * w_d)
365 estimated_time_constant = inductance / R_eq * 3
367 new_solve_object_for_inductance_and_time_constant["time"]["start"] = 0
368 new_solve_object_for_inductance_and_time_constant["time"][
369 "adaptiveSteppingSettings"
370 ]["initialStep"] = (estimated_time_constant / 18)
371 new_solve_object_for_inductance_and_time_constant["time"][
372 "adaptiveSteppingSettings"
373 ]["minimumStep"] = (estimated_time_constant / 512)
374 new_solve_object_for_inductance_and_time_constant["time"][
375 "adaptiveSteppingSettings"
376 ]["maximumStep"] = (estimated_time_constant / 18)
377 new_solve_object_for_inductance_and_time_constant["time"]["start"] = 0
378 new_solve_object_for_inductance_and_time_constant["time"]["end"] = (
379 0.1 + estimated_time_constant * 8 + 0.001
380 )
381 new_solve_object_for_inductance_and_time_constant["winding"] = (
382 copy_fdm.magnet.solve.wi.dict()
383 )
384 new_solve_object_for_inductance_and_time_constant["contactLayer"] = (
385 copy_fdm.magnet.solve.ii.dict()
386 )
387 new_solve_object_for_inductance_and_time_constant["terminals"] = (
388 copy_fdm.magnet.solve.ti.dict()
389 )
390 solve_object = Pancake3DSolve(
391 **new_solve_object_for_inductance_and_time_constant
392 )
393 solve_object.save = [
394 Pancake3DSolveSaveQuantity(quantity="timeConstant")
395 ]
397 solve_object.ii.resistivity = copy_fdm.magnet.solve.ii.resistivity
398 copy_fdm.magnet.solve = solve_object
400 new_power_supply_object_for_time_constant = {
401 "t_control_LUT": [
402 0,
403 0.1,
404 0.1 + estimated_time_constant * 6,
405 0.1 + estimated_time_constant * 6 + 0.001,
406 0.1 + estimated_time_constant * 8 + 0.001,
407 ],
408 "I_control_LUT": [0, 1, 1, 0, 0],
409 }
410 new_power_supply_object_for_time_constant = PowerSupply(
411 **new_power_supply_object_for_time_constant
412 )
413 copy_fdm.power_supply = new_power_supply_object_for_time_constant
415 time_constant_solution_folder = os.path.join(
416 self.solution_folder,
417 "time_constant",
418 )
419 # generate the folder if it does not exist:
420 if not os.path.exists(time_constant_solution_folder):
421 os.makedirs(time_constant_solution_folder)
422 solve = Solve(
423 copy_fdm,
424 self.GetDP_path,
425 self.geom_folder,
426 self.mesh_folder,
427 time_constant_solution_folder,
428 )
429 solve.assemble_pro()
431 solve.run_getdp(solve=True, postOperation=True)
433 # then compute the time constant and write it to a file:
434 # read axialComponentOfTheMagneticFieldForTimeConstant-TimeTableFormat.txt
436 time_constant_file = os.path.join(
437 time_constant_solution_folder,
438 "axialComponentOfTheMagneticFieldForTimeConstant-TimeTableFormat.txt",
439 )
440 times = []
441 Bzs = []
442 isNegative = False
443 with open(time_constant_file, "r") as f:
444 for line in f.readlines():
445 time = float(line.split()[1])
446 Bz = float(line.split()[-1])
447 times.append(time)
448 if Bz < 0:
449 isNegative = True
450 Bzs.append(-Bz)
451 else:
452 if isNegative:
453 raise ValueError(
454 "The magnetic field should not change sign."
455 )
456 Bzs.append(Bz)
458 # find the maximum value of Bz:
459 max_Bz = max(Bzs)
461 # during decay, find the time when Bz is 1/e of the maximum value (36.8%):
462 # Do linear interpolation if required:
463 percents_of_max_Bz = []
464 for Bz_i in Bzs:
465 percents_of_max_Bz.append(Bz_i / max_Bz)
467 peak_index = percents_of_max_Bz.index(max(percents_of_max_Bz))
468 for percent in percents_of_max_Bz[peak_index:]:
469 if percent < 0.368:
470 index = percents_of_max_Bz.index(percent)
471 break
473 # find the time when percents_of_max_Bz
474 x1 = times[index - 1]
475 x2 = times[index]
476 y1 = percents_of_max_Bz[index - 1]
477 y2 = percents_of_max_Bz[index]
479 time_constant = (
480 x1 + (0.368 - y1) * (x2 - x1) / (y2 - y1) - times[peak_index]
481 )
483 # write the time constant to a file:
484 time_constant_file = os.path.join(
485 time_constant_solution_folder, "time_constant.txt"
486 )
487 with open(time_constant_file, "w") as f:
488 f.write(str(time_constant))
490 # main solver:
491 solve = Solve(
492 self.fdm,
493 self.GetDP_path,
494 self.geom_folder,
495 self.mesh_folder,
496 self.solution_folder,
497 )
498 solve.assemble_pro()
500 start_time = Time.time()
501 solve.run_getdp(solve=True, postOperation=True)
503 return Time.time() - start_time
505 def post_process_getdp(self, gui=False):
506 solve = Solve(
507 self.fdm,
508 self.GetDP_path,
509 self.geom_folder,
510 self.mesh_folder,
511 self.solution_folder,
512 )
513 solve.assemble_pro()
515 start_time = Time.time()
516 solve.run_getdp(solve=False, postOperation=True)
518 return Time.time() - start_time
520 def post_process_python(self, gui=False):
521 """
522 To be written.
523 """
524 postprocess = Postprocess(
525 self.fdm,
526 self.geom_folder,
527 self.mesh_folder,
528 self.solution_folder,
529 )
531 if self.fdm.magnet.postproc is not None:
532 if self.fdm.magnet.postproc.timeSeriesPlots is not None:
533 postprocess.plotTimeSeriesPlots()
534 if self.fdm.magnet.postproc.magneticFieldOnCutPlane is not None:
535 postprocess.plotMagneticFieldOnCutPlane()
537 return {"overall_error": 0}
539 def plot_python(self):
540 pass