Coverage for fiqus/mains/MainPancake3D.py: 40%
226 statements
« prev ^ index » next coverage.py v7.4.4, created at 2025-03-11 02:30 +0100
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
114from fiqus.data.DataFiQuS import PowerSupply
117class MainPancake3D:
118 """
119 The main class for working with simulations for high-temperature superconductor
120 pancake coil magnets.
122 Geometry can be created and saved as a BREP file. Parameters like the number of
123 turns, tape dimensions, contact layer thicknesses, and other dimensions can be
124 specified. Contact layers can be modeled as two-dimensional shells or three-dimensional
125 volumes. Moreover, multiple pancakes can be stacked on top of each other.
127 Using the BREP file created, a mesh can be generated and saved as an MSH file.
128 Winding mesh can be structured, and parameters like, azimuthal number of elements
129 per turn, axial number of elements, and radial number of elements per turn can be
130 specified for each pancake coil. The appropriate regions will be assigned to the
131 relevant volumes accordingly so that finite element simulations can be done.
133 Using the mesh files, GetDP can be used to analyze Pancake3D coils.
135 :param fdm: FiQuS data model
136 """
138 def __init__(self, fdm, verbose):
139 self.fdm: FDM = fdm
140 self.GetDP_path = None
142 self.geom_folder = None
143 self.mesh_folder = None
144 self.solution_folder = None
146 def generate_geometry(self, gui=False):
147 """
148 Generates the geometry of the magnet and save it as a BREP file. Moreover, a
149 text file with the extension VI (volume information file) is generated, which
150 stores the names of the volume tags in JSON format.
151 """
152 geometry = Geometry(
153 self.fdm,
154 self.geom_folder,
155 self.mesh_folder,
156 self.solution_folder,
157 )
158 if self.fdm.magnet.geometry.conductorWrite:
159 self.fdm.magnet.geometry.gapBetweenPancakes = 0.00024
160 innerRadius = self.fdm.magnet.geometry.winding.innerRadius
161 gap = self.fdm.magnet.geometry.winding.thickness
162 # gap = 0.01
163 turns = self.fdm.magnet.geometry.numberOfPancakes
164 # turns = 4
165 z = - self.fdm.magnet.geometry.winding.height/2
166 initialTheta = 0.0
167 # transitionNotchAngle = self.fdm.magnet.geometry.wi.t ##TODO pull from internal recalculated
168 # transitionNotchAngle = 0.7853981633974483/2
169 d = direction(0)
170 sc = spiralCurve(innerRadius,
171 gap,
172 turns,
173 z,
174 initialTheta,
175 transitionNotchAngle,
176 self.fdm.magnet.geometry,
177 direction=d, # TODO code this to understand cw direction
178 cutPlaneNormal=None)
179 else:
180 geometry.generate_geometry()
181 geometry.generate_vi_file()
183 self.model_file = geometry.brep_file
185 def load_geometry(self, gui=False):
186 """
187 Loads the previously generated geometry from the BREP file.
188 """
189 geometry = Geometry(
190 self.fdm,
191 self.geom_folder,
192 self.mesh_folder,
193 self.solution_folder,
194 )
196 geometry.load_geometry()
197 self.model_file = geometry.brep_file
199 def pre_process(self, gui=False):
200 pass
202 def mesh(self, gui=False):
203 """
204 Generates the mesh of the magnet, creates the physical regions, and saves it as
205 an MSH file. Moreover, a text file with the extension REGIONS is generated,
206 which stores the names and tags of the physical regions in YAML format.
207 """
208 mesh = Mesh(
209 self.fdm,
210 self.geom_folder,
211 self.mesh_folder,
212 self.solution_folder,
213 )
215 mesh.generate_mesh()
216 mesh.generate_regions()
217 mesh.generate_mesh_file()
219 self.model_file = mesh.mesh_file
221 return {"gamma": 0} # to be modified with mesh_parameters (see multipole)
223 def load_mesh(self, gui=False):
224 """
225 Loads the previously generated mesh from the MSH file.
226 """
227 mesh = Mesh(
228 self.fdm,
229 self.geom_folder,
230 self.mesh_folder,
231 self.solution_folder,
232 )
233 mesh.load_mesh()
235 self.model_file = mesh.mesh_file
237 def solve_and_postprocess_getdp(self, gui=False):
238 """
239 Simulates the Pancake3D magnet with GetDP using the created mesh file and post
240 processes the results.
241 """
242 # calculate inductance and time constant:
243 if self.fdm.magnet.solve.quantitiesToBeSaved is not None:
244 save_list = list(
245 map(
246 lambda quantity_object: quantity_object.quantity,
247 self.fdm.magnet.solve.quantitiesToBeSaved,
248 )
249 )
250 new_solve_object_for_inductance_and_time_constant = {
251 "type": "electromagnetic",
252 "time": {
253 "timeSteppingType": "adaptive",
254 "extrapolationOrder": 1,
255 "adaptiveSteppingSettings": {
256 "initialStep": 1,
257 "minimumStep": 0.00001,
258 "maximumStep": 20,
259 "integrationMethod": "Euler",
260 "tolerances": [
261 {
262 "quantity": "magnitudeOfMagneticField",
263 "position": {
264 "x": 0,
265 "y": 0,
266 "z": 0,
267 },
268 "relative": 0.05, # 5%
269 "absolute": 0.00002, # 0.00002 T
270 "normType": "LinfNorm",
271 }
272 ],
273 },
274 },
275 "nonlinearSolver": self.fdm.magnet.solve.nonlinearSolver,
276 "air": {
277 "permeability": 1.2566e-06,
278 },
279 "initialConditions": {
280 "temperature": 4,
281 },
282 }
283 inductance = None
284 inductance_solution_folder = os.path.join(
285 self.solution_folder,
286 "inductance",
287 )
288 inductance_file = os.path.join(
289 inductance_solution_folder, "Inductance-TimeTableFormat.csv"
290 )
291 if "inductance" in save_list:
292 # to calculate inductance and time constant, we should run the simulation
293 # with specigic power supply and material settings.
294 copy_fdm = deepcopy(self.fdm)
296 new_solve_object_for_inductance_and_time_constant["time"]["start"] = 0
297 new_solve_object_for_inductance_and_time_constant["time"]["end"] = 100
298 new_solve_object_for_inductance_and_time_constant["winding"] = {
299 "resistivity": 1e-20,
300 "thermalConductivity": 1,
301 "specificHeatCapacity": 1,
302 }
304 new_solve_object_for_inductance_and_time_constant["contactLayer"] = {
305 "resistivity": 1e-2,
306 "thermalConductivity": 1,
307 "specificHeatCapacity": 1,
308 "numberOfThinShellElements": 1,
309 }
310 new_solve_object_for_inductance_and_time_constant["terminals"] = {
311 "resistivity": 1e-10,
312 "thermalConductivity": 1,
313 "specificHeatCapacity": 1,
314 "terminalContactLayer": {
315 "resistivity": 1e-10,
316 "thermalConductivity": 1,
317 "specificHeatCapacity": 1,
318 },
319 "transitionNotch": {
320 "resistivity": 1e-10,
321 "thermalConductivity": 1,
322 "specificHeatCapacity": 1,
323 },
324 }
325 solve_object = Pancake3DSolve(
326 **new_solve_object_for_inductance_and_time_constant
327 )
328 solve_object.quantitiesToBeSaved = [Pancake3DSolveSaveQuantity(quantity="inductance")]
329 copy_fdm.magnet.solve = solve_object
331 new_power_supply_object_for_inductance = {
332 "t_control_LUT": [0, 0.1, 100],
333 "I_control_LUT": [0, 1, 1],
334 }
335 new_power_supply_object_for_inductance = PowerSupply(
336 **new_power_supply_object_for_inductance
337 )
338 copy_fdm.power_supply = new_power_supply_object_for_inductance
340 # generate the folder if it does not exist:
341 if not os.path.exists(inductance_solution_folder):
342 os.makedirs(inductance_solution_folder)
343 else:
344 shutil.rmtree(inductance_solution_folder)
345 os.makedirs(inductance_solution_folder)
347 solve = Solve(
348 copy_fdm,
349 self.GetDP_path,
350 self.geom_folder,
351 self.mesh_folder,
352 inductance_solution_folder,
353 )
354 solve.assemble_pro()
356 start_time = Time.time()
357 solve.run_getdp(solve=True, postOperation=True)
359 with open(inductance_file, "r") as f:
360 # read the last line and second column:
361 inductance = float(f.readlines()[-1].split()[1].replace(",", ""))
363 with open(inductance_file, "w") as f:
364 f.write(str(inductance))
366 if "timeConstant" in save_list:
367 # to calculate inductance and time constant, we should run the simulation
368 # with specigic power supply and material settings.
369 copy_fdm = deepcopy(self.fdm)
371 if inductance is None:
372 raise ValueError(
373 "Time constant can not be calculated without inductance. Please"
374 " add 'inductance' to the 'quantitiesToBeSaved' list in the"
375 " 'solve' section of the YAML file."
376 )
378 # equivalent radial resistance of the winding:
379 # https://doi.org/10.1109/TASC.2021.3063653
380 R_eq = 0
381 for i in range(1, int(copy_fdm.magnet.geometry.winding.numberOfTurns)):
382 r_i = (
383 copy_fdm.magnet.geometry.winding.innerRadius
384 + i * copy_fdm.magnet.geometry.winding.thickness
385 + (i - 1) * copy_fdm.magnet.geometry.contactLayer.thickness
386 + copy_fdm.magnet.geometry.contactLayer.thickness / 2
387 )
388 w_d = copy_fdm.magnet.geometry.winding.height
389 R_ct = (
390 copy_fdm.magnet.solve.contactLayer.resistivity
391 * copy_fdm.magnet.geometry.contactLayer.thickness
392 )
393 R_eq = R_eq + R_ct / (2 * math.pi * r_i * w_d)
395 estimated_time_constant = inductance / R_eq * 3
397 new_solve_object_for_inductance_and_time_constant["time"]["start"] = 0
398 new_solve_object_for_inductance_and_time_constant["time"][
399 "adaptiveSteppingSettings"
400 ]["initialStep"] = (estimated_time_constant / 18)
401 new_solve_object_for_inductance_and_time_constant["time"][
402 "adaptiveSteppingSettings"
403 ]["minimumStep"] = (estimated_time_constant / 512)
404 new_solve_object_for_inductance_and_time_constant["time"][
405 "adaptiveSteppingSettings"
406 ]["maximumStep"] = (estimated_time_constant / 18)
407 new_solve_object_for_inductance_and_time_constant["time"]["start"] = 0
408 new_solve_object_for_inductance_and_time_constant["time"]["end"] = (
409 0.1 + estimated_time_constant * 8 + 0.001
410 )
411 new_solve_object_for_inductance_and_time_constant["winding"] = (
412 copy_fdm.magnet.solve.winding.dict()
413 )
414 new_solve_object_for_inductance_and_time_constant["contactLayer"] = (
415 copy_fdm.magnet.solve.contactLayer.dict()
416 )
417 new_solve_object_for_inductance_and_time_constant["terminals"] = (
418 copy_fdm.magnet.solve.terminals.dict()
419 )
420 solve_object = Pancake3DSolve(
421 **new_solve_object_for_inductance_and_time_constant
422 )
423 solve_object.quantitiesToBeSaved = [
424 Pancake3DSolveSaveQuantity(quantity="timeConstant")
425 ]
427 solve_object.contactLayer.resistivity = copy_fdm.magnet.solve.contactLayer.resistivity
428 copy_fdm.magnet.solve = solve_object
430 new_power_supply_object_for_time_constant = {
431 "t_control_LUT": [
432 0,
433 0.1,
434 0.1 + estimated_time_constant * 6,
435 0.1 + estimated_time_constant * 6 + 0.001,
436 0.1 + estimated_time_constant * 8 + 0.001,
437 ],
438 "I_control_LUT": [0, 1, 1, 0, 0],
439 }
440 new_power_supply_object_for_time_constant = PowerSupply(
441 **new_power_supply_object_for_time_constant
442 )
443 copy_fdm.power_supply = new_power_supply_object_for_time_constant
445 time_constant_solution_folder = os.path.join(
446 self.solution_folder,
447 "time_constant",
448 )
449 # generate the folder if it does not exist:
450 if not os.path.exists(time_constant_solution_folder):
451 os.makedirs(time_constant_solution_folder)
452 solve = Solve(
453 copy_fdm,
454 self.GetDP_path,
455 self.geom_folder,
456 self.mesh_folder,
457 time_constant_solution_folder,
458 )
459 solve.assemble_pro()
461 solve.run_getdp(solve=True, postOperation=True)
463 # then compute the time constant and write it to a file:
464 # read axialComponentOfTheMagneticFieldForTimeConstant-TimeTableFormat.txt
466 time_constant_file = os.path.join(
467 time_constant_solution_folder,
468 "axialComponentOfTheMagneticFieldForTimeConstant-TimeTableFormat.csv",
469 )
470 times = []
471 Bzs = []
472 isNegative = False
473 with open(time_constant_file, "r") as f:
474 for line in f.readlines():
475 time = float(line.split()[1])
476 Bz = float(line.split()[-1])
477 times.append(time)
478 if Bz < 0:
479 isNegative = True
480 Bzs.append(-Bz)
481 else:
482 if isNegative:
483 raise ValueError(
484 "The magnetic field should not change sign."
485 )
486 Bzs.append(Bz)
488 # find the maximum value of Bz:
489 max_Bz = max(Bzs)
491 # during decay, find the time when Bz is 1/e of the maximum value (36.8%):
492 # Do linear interpolation if required:
493 percents_of_max_Bz = []
494 for Bz_i in Bzs:
495 percents_of_max_Bz.append(Bz_i / max_Bz)
497 peak_index = percents_of_max_Bz.index(max(percents_of_max_Bz))
498 for percent in percents_of_max_Bz[peak_index:]:
499 if percent < 0.368:
500 index = percents_of_max_Bz.index(percent)
501 break
503 # find the time when percents_of_max_Bz
504 x1 = times[index - 1]
505 x2 = times[index]
506 y1 = percents_of_max_Bz[index - 1]
507 y2 = percents_of_max_Bz[index]
509 time_constant = (
510 x1 + (0.368 - y1) * (x2 - x1) / (y2 - y1) - times[peak_index]
511 )
513 # write the time constant to a file:
514 time_constant_file = os.path.join(
515 time_constant_solution_folder, "time_constant.csv"
516 )
517 with open(time_constant_file, "w") as f:
518 f.write(str(time_constant))
520 # main solver:
521 if "inductance" in save_list:
522 for save_object in self.fdm.magnet.solve.quantitiesToBeSaved:
523 if save_object.quantity == "inductance":
524 self.fdm.magnet.solve.quantitiesToBeSaved.remove(save_object)
526 if "timeConstant" in save_list:
527 for save_object in self.fdm.magnet.solve.quantitiesToBeSaved:
528 if save_object.quantity == "timeConstant":
529 self.fdm.magnet.solve.quantitiesToBeSaved.remove(save_object)
531 solve = Solve(
532 self.fdm,
533 self.GetDP_path,
534 self.geom_folder,
535 self.mesh_folder,
536 self.solution_folder,
537 )
538 solve.assemble_pro()
540 start_time = Time.time()
541 solve.run_getdp(solve=True, postOperation=True)
543 return Time.time() - start_time
545 def post_process_getdp(self, gui=False):
546 solve = Solve(
547 self.fdm,
548 self.GetDP_path,
549 self.geom_folder,
550 self.mesh_folder,
551 self.solution_folder,
552 )
553 solve.assemble_pro()
555 start_time = Time.time()
556 solve.run_getdp(solve=False, postOperation=True)
558 return Time.time() - start_time
560 def post_process_python(self, gui=False):
561 """
562 To be written.
563 """
564 postprocess = Postprocess(
565 self.fdm,
566 self.geom_folder,
567 self.mesh_folder,
568 self.solution_folder,
569 )
571 if self.fdm.magnet.postproc is not None:
572 if self.fdm.magnet.postproc.timeSeriesPlots is not None:
573 postprocess.plotTimeSeriesPlots()
574 if self.fdm.magnet.postproc.magneticFieldOnCutPlane is not None:
575 postprocess.plotMagneticFieldOnCutPlane()
577 return {"overall_error": 0}
579 def plot_python(self):
580 pass