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