Coverage for python/gsfit/database_readers/st40_astra_mdsplus/setup_coils.py: 0%
134 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-14 15:25 +0000
1# mypy: ignore-errors
2# TODO: need to fix mypy errors
4import typing
5from typing import TYPE_CHECKING
7import mdsthin
8import numpy as np
9import numpy.typing as npt
10from gsfit_rs import Coils
12from .astra_coils_reader import astra_coils_reader
14if TYPE_CHECKING:
15 from . import DatabaseReader
17# Read coil geometry from coil.dat
18astra_coil_description = astra_coils_reader()
21def dividing_parallograms(coil_dictionary) -> tuple[int, int]:
22 """
23 The PF coils are described as parallograms, and discretized in `ASTRA/SRC/Trecur_2.f`, in the `equil_astra_tepm` GitLab repository.
25 ```
26 SW = SQRT( NDIVA*WSIZE/HSIZE )
27 SH = SQRT( NDIVA*HSIZE/WSIZE )
28 SW = SW + 0.5
29 SH = SH + 0.5
30 C
31 NDIVW = IDINT(SW)
32 NDIVH = IDINT(SH)
33 ```
34 """
35 ndi = coil_dictionary["ndi"]
36 wc = coil_dictionary["wc"]
37 hc = coil_dictionary["hc"]
39 sw = np.sqrt(ndi * wc / hc)
40 sh = np.sqrt(ndi * hc / wc)
41 sw = sw + 0.5
42 sh = sh + 0.5
43 ndivw = int(sw)
44 ndivh = int(sh)
45 return ndivw, ndivh
48def calculate_coil_filament_positions(coil_dictionary) -> tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]:
49 """
50 Calculate the coil filament positions from the coil description dictionary.
52 :param coil_dictionary: Dictionary containing the coil description
53 :return: Tuple of numpy arrays containing the r and z positions of the coil filaments
55 ```
56 C**************************************************************
57 C SUBROUTINE FOR DIVIDING OF PARALLELOGRAM
58 C ( FOR PFC CROSS-SECTION )
59 C**************************************************************
60 C INPUT DATE:
61 C ----------
62 C RC,ZC - CILINDER COORDINATES OF CENTER OF PARALLELOGRAM
63 C WC - PROJECTION OF THE FIRST SIDE OF PARALLELOGRAM
64 C ON AXIS "R" (IN METER)
65 C HC - PROJECTION OF THE SECOND SIDE OF PARALLELOGRAM
66 C ON AXIS "Z" (IN METER)
67 C AWC - ANGLE BETWEEN THE FIRST SIDE OF PARAL. AND
68 C AXIS "R" (IN DEGREES, IT MUST NOT BE EQUAL 90 ,
69 C IT MUST BE : -90 < AWC < 90 )
70 C AHC - ANGLE BETWEEN THE SECOND SIDE OF PARAL. AND
71 C AXIS "R" (IN DEGREES, IT MUST NOT BE EQUAL 0 ,
72 C IT MUST BE : 0 < AHC < 180 )
73 C CURC - CURRENT OF PARALLELOGRAM (IN MA)
74 C NDIV - APPROXIMATE NUMBER OF CELLS OF DIVIDING:
75 C NDIV=0 - A SPECIAL CASE: AUTOMATICALLY NDIVRE=1,
76 C RS(1)=RC, ZS(1)=ZC, PS=CURC
77 C IF NDIV > 0 THEN WE HAVE THE MOST TOTAL ALGORITHM OF
78 C DIVIDING
79 C IF NDIV < 0 THEN NDIVW AND NDIVH ARE CUT OFF BY ABS(NDIV)
80 C
81 C OUTPUT DATE:
82 C ----------
83 C NDIVRE - REAL NUMBER OF CELLS OF DIVIDING ( = NDIVW*NDIVH )
84 C NDIVW - NUMBER OF DIVIDING OF THE FIRST SIDE OF PARAL.
85 C NDIVH - NUMBER OF DIVIDING OF THE SECOND SIDE OF PARAL.
86 C RS(L),ZS(L) - CILINDER COORDINATES OF CENTERS OF CELLS OF DIVIDING
87 C L = 1,2,...,NDIVRE ! (IN METER)
88 C PS - CURRENT OF EVERY CELL OF DIVIDING (IN MA)
89 C VERS - VERTICAL (OR LINEAR, OR RADIUS) SIZE OF CELL CROSS-S.
90 C HORS - HORIZONTAL SIZE OF CELL CROSS-SECTION
91 C
92 C**************************************************************
93 C
94 SUBROUTINE DIVPAR( RC, ZC, WC, HC, AWC, AHC, CURC, NDIV,
95 * NDIVRE, NDIVW, NDIVH, RS, ZS, PS,
96 * VERS, HORS )
97 C
98 C
99 include 'double.inc'
100 C
101 DIMENSION RS(1), ZS(1)
102 C
103 C**************************************************************
104 SIN(X) = DSIN(X)
105 COS(X) = DCOS(X)
106 ATAN(X) = DATAN(X)
107 SQRT(X) = DSQRT(X)
108 C**************************************************************
109 C
110 IF(NDIV.EQ.0) THEN
111 NDIVW = 1
112 NDIVH = 1
113 NDIVRE = 1
114 RS(1) = RC
115 ZS(1) = ZC
116 PS = CURC
117 VERS = HC
118 HORS = WC
119 RETURN
120 END IF
121 C***************************************
122 C
123 IF(AHC.LT.0) AHC = AHC + 180.
124 C
125 IF(NDIV.GE.0) THEN
126 NDIVA = NDIV
127 ELSE
128 NDIVA = -NDIV
129 END IF
130 C***************************************
131 C PARAMETERS OF PARALLELOGRAM
132 C
133 XX = 1.
134 PI = 4.*ATAN(XX)
135 C
136 AWCR = AWC * PI /180.
137 AHCR = AHC * PI /180.
138 C
139 R0 = RC - 0.5*( WC + HC * COS(AHCR)/SIN(AHCR) )
140 Z0 = ZC - 0.5*( HC + WC * SIN(AWCR)/COS(AWCR) )
141 C
142 WSIZE = WC / COS(AWCR)
143 HSIZE = HC / SIN(AHCR)
144 C
145 WR = WC
146 WZ = WC * SIN(AWCR) / COS(AWCR)
147 HR = HC * COS(AHCR) / SIN(AHCR)
148 HZ = HC
149 C***************************************
150 C CALCULATION NDIVW, NDIVH, NDIVRE, PS
151 C
152 SW = SQRT( NDIVA*WSIZE/HSIZE )
153 SH = SQRT( NDIVA*HSIZE/WSIZE )
154 SW = SW + 0.5
155 SH = SH + 0.5
156 C
157 NDIVW = IDINT(SW)
158 NDIVH = IDINT(SH)
159 C
160 IF(NDIVW.EQ.0) NDIVW = 1
161 IF(NDIVH.EQ.0) NDIVH = 1
162 C
163 IF((NDIV.LT.0).AND.(NDIVW.GT.NDIVA)) NDIVW = NDIVA
164 IF((NDIV.LT.0).AND.(NDIVH.GT.NDIVA)) NDIVH = NDIVA
165 C
166 NDIVRE = NDIVW * NDIVH
167 PS = CURC / NDIVRE
168 C***************************************
169 C CALCULATION RS(L), ZS(L) : L = 1,2,...,NDIVRE
170 C
171 WR = WR / NDIVW
172 WZ = WZ / NDIVW
173 HR = HR / NDIVH
174 HZ = HZ / NDIVH
175 C
176 HORS = WR
177 VERS = HZ
178 C
179 RS(1) = R0 + 0.5*(WR + HR)
180 ZS(1) = Z0 + 0.5*(WZ + HZ)
181 C
182 DO 1 I=1,NDIVW
183 DO 1 J=1,NDIVH
184 L = (I-1)*NDIVH + J
185 RS(L) = RS(1) + (I-1)*WR + (J-1)*HR
186 ZS(L) = ZS(1) + (I-1)*WZ + (J-1)*HZ
187 1 CONTINUE
188 C***************************************
189 C
190 RETURN
191 END
192 ```
193 """
195 # Extract required parameters
196 rc = float(coil_dictionary["rc"]) # center R
197 zc = float(coil_dictionary["zc"]) # center Z
198 wc = float(coil_dictionary["wc"]) # projection of first side on R
199 hc = float(coil_dictionary["hc"]) # projection of second side on Z
200 # Optional angles (degrees). Default to AWC=0, AHC=90 if not provided,
201 # matching the most common cases in the coil table.
202 awc = float(coil_dictionary.get("awc", 0.0))
203 ahc = float(coil_dictionary.get("ahc", 90.0))
204 # Approximate number of cells (NDIV/NDIVA)
205 ndia = int(coil_dictionary.get("ndi", 1))
207 # Handle the special case NDIV == 0
208 if ndia == 0:
209 rs = np.array([rc], dtype=np.float64)
210 zs = np.array([zc], dtype=np.float64)
211 return rs, zs
213 # Convert angles to radians, applying the Fortran adjustment for AHC < 0
214 if ahc < 0.0:
215 ahc = ahc + 180.0
217 pi = np.pi
218 awcr = awc * pi / 180.0
219 ahcr = ahc * pi / 180.0
221 # Guard against invalid angles (avoid division by zero)
222 # In Fortran code: it assumes -90 < AWC < 90 and 0 < AHC < 180
223 # If cos(awcr) or sin(ahcr) are too close to zero, nudge slightly.
224 eps = 1e-12
225 c_aw = np.cos(awcr)
226 s_aw = np.sin(awcr)
227 s_ah = np.sin(ahcr)
228 c_ah = np.cos(ahcr)
229 if abs(c_aw) < eps:
230 c_aw = eps if c_aw >= 0 else -eps
231 if abs(s_ah) < eps:
232 s_ah = eps if s_ah >= 0 else -eps
234 # Fortran mapping
235 r0 = rc - 0.5 * (wc + hc * c_ah / s_ah)
236 z0 = zc - 0.5 * (hc + wc * s_aw / c_aw)
238 wsize = wc / c_aw
239 hsize = hc / s_ah
241 wr = wc
242 wz = wc * s_aw / c_aw
243 hr = hc * c_ah / s_ah
244 hz = hc
246 # Calculate NDIVW, NDIVH per Fortran's DIVPAR (NDIV >= 0 path)
247 sw = np.sqrt(ndia * wsize / hsize) + 0.5
248 sh = np.sqrt(ndia * hsize / wsize) + 0.5
249 ndivw = int(sw)
250 ndivh = int(sh)
251 if ndivw == 0:
252 ndivw = 1
253 if ndivh == 0:
254 ndivh = 1
256 # Real number of divided cells
257 ndivre = ndivw * ndivh
259 # Cell sizes
260 wr_cell = wr / ndivw
261 wz_cell = wz / ndivw
262 hr_cell = hr / ndivh
263 hz_cell = hz / ndivh
265 # First cell center
266 r1 = r0 + 0.5 * (wr_cell + hr_cell)
267 z1 = z0 + 0.5 * (wz_cell + hz_cell)
269 # Generate centers for all cells (Fortran order L = (I-1)*NDIVH + J)
270 rs = np.empty(ndivre, dtype=np.float64)
271 zs = np.empty(ndivre, dtype=np.float64)
272 idx = 0
273 for i in range(1, ndivw + 1):
274 for j in range(1, ndivh + 1):
275 r = r1 + (i - 1) * wr_cell + (j - 1) * hr_cell
276 z = z1 + (i - 1) * wz_cell + (j - 1) * hz_cell
277 rs[idx] = r
278 zs[idx] = z
279 idx += 1
281 return rs, zs
284def setup_coils(
285 self: "DatabaseReader",
286 settings: dict[str, typing.Any],
287 pulseNo: int,
288) -> Coils:
289 """
290 This method initialises the Rust `Coils` class.
292 :param pulseNo: Pulse number, used to read from the database
293 :param settings: Dictionary containing the JSON settings read from the `settings` directory
295 **This method is specific to ST40's ASTRA stored on MDSplus.**
297 See `python/gsfit/database_readers/interface.py` for more details on how a new database_reader should be implemented.
298 """
300 # Initialise the Coils Rust class
301 coils = Coils()
303 # Extract the astra_run_name from settings
304 astra_run_name = settings["GSFIT_code_settings.json"]["database_reader"]["st40_astra_mdsplus"]["workflow"]["astra"]["run_name"]
306 # Connect to MDSplus
307 conn = mdsthin.Connection("smaug")
308 conn.openTree("ASTRA", pulseNo)
310 # ASTRA time
311 time = conn.get(f"\\ASTRA::TOP.{astra_run_name}:TIME").data().astype(np.float64)
313 # PF coil currents
314 # ["IPL", "MC", "PSH", "DIV", "BVL", "BVUT", "BVUB", "CS", "MCVC"]
315 # [ 0, 1, 2, 3, 4, 5, 6, 7 8]
316 currents = conn.get(f"\\ASTRA::TOP.{astra_run_name}.COILS.PSU2PF:I").data().astype(np.float64) * 1.0e6
318 # Add BVLB PF coil
319 coil_name = "BVLB"
320 coil_r, coil_z = calculate_coil_filament_positions(astra_coil_description[coil_name])
321 # Read in BVL PSU current, and normalise by the number of filaments
322 currents_local = currents[:, 4] * 4.0 * 4.0 / len(coil_r)
323 coils.add_pf_coil(
324 coil_name,
325 coil_r,
326 coil_z,
327 d_r=0.0 * coil_r,
328 d_z=0.0 * coil_z,
329 time=time,
330 measured=currents_local,
331 )
332 # Add BVLT PF coil
333 coil_name = "BVLT"
334 coil_r, coil_z = calculate_coil_filament_positions(astra_coil_description[coil_name])
335 coils.add_pf_coil(
336 coil_name,
337 coil_r,
338 coil_z,
339 d_r=0.0 * coil_r,
340 d_z=0.0 * coil_z,
341 time=time,
342 measured=currents_local,
343 )
345 # Read in BVUB PSU current, and normalise by the number of filaments
346 coil_name = "BVUB"
347 coil_r, coil_z = calculate_coil_filament_positions(astra_coil_description[coil_name])
348 currents_local = currents[:, 6] * 4.0 * 6.0 / len(coil_r)
349 coils.add_pf_coil(
350 coil_name,
351 coil_r,
352 coil_z,
353 d_r=0.0 * coil_r,
354 d_z=0.0 * coil_z,
355 time=time,
356 measured=currents_local,
357 )
359 # Read in BVUT PSU current, and normalise by the number of filaments
360 coil_name = "BVUT"
361 coil_r, coil_z = calculate_coil_filament_positions(astra_coil_description[coil_name])
362 currents_local = currents[:, 5] * 4.0 * 6.0 / len(coil_r)
363 coils.add_pf_coil(
364 coil_name,
365 coil_r,
366 coil_z,
367 d_r=0.0 * coil_r,
368 d_z=0.0 * coil_z,
369 time=time,
370 measured=currents_local,
371 )
373 # Read in SOL PSU current, and normalise by the number of filaments
374 coil_name = "SOL"
375 coil_r, coil_z = calculate_coil_filament_positions(astra_coil_description[coil_name])
376 currents_local = currents[:, 7] * 2.0 * 95.0 / len(coil_r)
377 coils.add_pf_coil(
378 coil_name,
379 coil_r,
380 coil_z,
381 d_r=0.0 * coil_r,
382 d_z=0.0 * coil_z,
383 time=time,
384 measured=currents_local,
385 )
387 # Add DIVT PF coil, and normalise by the number of filaments
388 coil_name = "DIVT"
389 coil_r, coil_z = calculate_coil_filament_positions(astra_coil_description[coil_name])
390 currents_local = currents[:, 3] * 4.0 * 7.0 / len(coil_r)
391 coils.add_pf_coil(
392 coil_name,
393 coil_r,
394 coil_z,
395 d_r=0.0 * coil_r,
396 d_z=0.0 * coil_z,
397 time=time,
398 measured=currents_local,
399 )
400 # Add DIVB PF coil
401 coil_name = "DIVB"
402 coil_r, coil_z = calculate_coil_filament_positions(astra_coil_description[coil_name])
403 coils.add_pf_coil(
404 coil_name,
405 coil_r,
406 coil_z,
407 d_r=0.0 * coil_r,
408 d_z=0.0 * coil_z,
409 time=time,
410 measured=currents_local,
411 )
413 # Add MCT PF coil, and normalise by the number of filaments
414 coil_name = "MCT"
415 coil_r, coil_z = calculate_coil_filament_positions(astra_coil_description[coil_name])
416 currents_local = -currents[:, 8] * 11.0 / len(coil_r) / 2.0 + currents[:, 1] * 11.0 / len(coil_r)
417 coils.add_pf_coil(
418 coil_name,
419 coil_r,
420 coil_z,
421 d_r=0.0 * coil_r,
422 d_z=0.0 * coil_z,
423 time=time,
424 measured=currents_local,
425 )
426 # Add MCB PF coil
427 coil_name = "MCB"
428 coil_r, coil_z = calculate_coil_filament_positions(astra_coil_description[coil_name])
429 currents_local = currents[:, 8] * 11.0 / len(coil_r) / 2.0 + currents[:, 1] * 11.0 / len(coil_r)
430 coils.add_pf_coil(
431 coil_name,
432 coil_r,
433 coil_z,
434 d_r=0.0 * coil_r,
435 d_z=0.0 * coil_z,
436 time=time,
437 measured=currents_local,
438 )
440 # Add PSHT PF coil
441 coil_name = "PSHT1"
442 coil_r, coil_z = calculate_coil_filament_positions(astra_coil_description[coil_name])
443 currents_local = currents[:, 2] * 4.0 / len(coil_r)
444 coils.add_pf_coil(
445 coil_name,
446 coil_r,
447 coil_z,
448 d_r=0.0 * coil_r,
449 d_z=0.0 * coil_z,
450 time=time,
451 measured=currents_local,
452 )
453 # Add PSHB PF coil
454 coil_name = "PSHB1"
455 coil_r, coil_z = calculate_coil_filament_positions(astra_coil_description[coil_name])
456 coils.add_pf_coil(
457 coil_name,
458 coil_r,
459 coil_z,
460 d_r=0.0 * coil_r,
461 d_z=0.0 * coil_z,
462 time=time,
463 measured=currents_local,
464 )
465 # Add PSHT2 PF coil
466 coil_name = "PSHT2"
467 coil_r, coil_z = calculate_coil_filament_positions(astra_coil_description[coil_name])
468 coils.add_pf_coil(
469 coil_name,
470 coil_r,
471 coil_z,
472 d_r=0.0 * coil_r,
473 d_z=0.0 * coil_z,
474 time=time,
475 measured=currents_local,
476 )
477 # Add PSHB2 PF coil
478 coil_name = "PSHB2"
479 coil_r, coil_z = calculate_coil_filament_positions(astra_coil_description[coil_name])
480 coils.add_pf_coil(
481 coil_name,
482 coil_r,
483 coil_z,
484 d_r=0.0 * coil_r,
485 d_z=0.0 * coil_z,
486 time=time,
487 measured=currents_local,
488 )
490 return coils