Coverage for python/gsfit/database_writers/rtgsfit_mdsplus/map_results_to_database.py: 0%
303 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-14 15:25 +0000
1import copy
2from typing import TYPE_CHECKING
4import numpy as np
5import numpy.typing as npt
6import shapely.geometry
7from gsfit_rs.imas import equilibrium_paths as ep
8from gsfit_rs.imas import wall_paths as wp
9from scipy.constants import mu_0
11from .greens_with_boundary_points import greens_with_boundary_points
12from .poisson_matrix import compute_lup_bands
14if TYPE_CHECKING:
15 from ...gsfit import Gsfit
16 from . import DatabaseWriterRTGSFitMDSplus
19def map_results_to_database(self: "DatabaseWriterRTGSFitMDSplus", gsfit_controller: "Gsfit") -> None:
20 """
21 Map the results to MDSplus structure.
23 `gsfit_controller.results` is a `NestedDict` object, which has a 1:1 mapping to the MDSplus structure.
25 This function mutates the `gsfit_controller` object.
26 """
28 print("rtgsfit_mdsplus running")
30 # TODO: move this to a *.json file. This format hasn't yet settled down...
31 # TODO: alphabetical order
32 rtgsfit_psus = [
33 {"power_supply_name": "SOL", "coils": ["SOL"]},
34 {"power_supply_name": "MCT", "coils": ["MCT"]},
35 {"power_supply_name": "MCB", "coils": ["MCB"]},
36 {"power_supply_name": "DIV", "coils": ["DIVT", "DIVB"]},
37 {"power_supply_name": "BVL", "coils": ["BVLT", "BVLB"]},
38 {"power_supply_name": "BVUT", "coils": ["BVUT"]},
39 {"power_supply_name": "BVUB", "coils": ["BVUB"]},
40 {"power_supply_name": "PSH", "coils": ["PSHT", "PSHB"]},
41 ]
43 # Get objects out of `gsfit_controller`
44 plasma = gsfit_controller.plasma
45 passives = gsfit_controller.passives
46 flux_loops = gsfit_controller.flux_loops
47 bp_probes = gsfit_controller.bp_probes
48 rogowski_coils = gsfit_controller.rogowski_coils
49 results = gsfit_controller.results
51 wall = gsfit_controller.wall
53 # Geometry. `profiles_2d(0)` because GSFit solves on a single rectangular (R, Z) grid
54 r = plasma.get(ep.time_slice[0].profiles_2d[0].grid.dim1)
55 z = plasma.get(ep.time_slice[0].profiles_2d[0].grid.dim2)
56 n_z = len(z)
57 n_r = len(r)
58 d_r = np.mean(r[1:] - r[0:-1])
59 d_z = np.mean(z[1:] - z[0:-1])
61 # Number of power supplies
62 n_psu = len(rtgsfit_psus)
64 # Store geomery
65 results["PRESHOT"]["DR"] = d_r
66 results["PRESHOT"]["DZ"] = d_z
67 results["PRESHOT"]["R_VEC"] = r
68 results["PRESHOT"]["Z_VEC"] = z
69 results["PRESHOT"]["N_R"] = np.int32(n_r)
70 results["PRESHOT"]["N_Z"] = np.int32(n_z)
72 r_grid, z_grid = np.meshgrid(r, z)
73 r_flat = r_grid.flatten()
74 z_flat = z_grid.flatten()
75 results["PRESHOT"]["R_GRID"] = r_flat
76 results["PRESHOT"]["Z_GRID"] = z_flat
77 inv_r_mu0 = 1.0 / (r_flat * mu_0)
78 results["PRESHOT"]["INV_R_MU0"] = inv_r_mu0
79 r_mu0_dz2 = 2 * np.pi * r_flat * mu_0 * d_z**2
80 results["PRESHOT"]["R_MU0_DZ2"] = r_mu0_dz2
81 results["PRESHOT"]["N_COIL"] = np.int32(n_psu)
82 n_grid = n_r * n_z
83 results["PRESHOT"]["N_GRID"] = np.int32(n_grid)
84 n_ltrb = 2 * n_r + 2 * n_z - 4 # Number of points on the boundary, removing the double counting at the 4 corners
85 results["PRESHOT"]["N_LTRB"] = np.int32(n_ltrb)
87 # Collect the greens for "grid-coils"
88 g_grid_coil = np.zeros((n_z, n_r, n_psu))
89 # The Green's tables are an array of structures, so the coil name is looked up once to give an
90 # index to read by
91 pf_active_names: list[str] = plasma.get(ep.greens.pf_active[:].name).tolist()
92 psu_names = []
93 for i_psu, power_supply in enumerate(rtgsfit_psus):
94 psu_names.append(power_supply["power_supply_name"])
95 coil_names = power_supply["coils"]
96 for coil_name in coil_names:
97 i_coil = pf_active_names.index(coil_name)
98 g_grid_coil[:, :, i_psu] += plasma.get(ep.greens.pf_active[i_coil].psi)
100 # Store in MDSplus
101 results["PRESHOT"]["GREENS"]["GRID_COIL"] = g_grid_coil.flatten()
102 results["PRESHOT"]["COIL_NAMES"] = np.array(psu_names)
104 # Get the "included" sensors
105 flux_loops_to_include = flux_loops.get_vec_bool(["*", "fit_settings", "include"])
106 n_flux_loops_to_include: int = sum(flux_loops_to_include)
107 bp_probes_to_include = bp_probes.get_vec_bool(["*", "fit_settings", "include"])
108 n_bp_probes_to_include: int = sum(bp_probes_to_include)
109 rogowski_coils_to_include = rogowski_coils.get_vec_bool(["*", "fit_settings", "include"])
110 n_rogowski_coils_to_include: int = sum(rogowski_coils_to_include)
112 # Count the number of passive degrees of freedom, and regularisations
113 n_passive_dofs = 0
114 n_regularisations = 0
115 passive_names = passives.keys() # Note: this includes the IVC
116 print(passive_names)
117 # passive_names = ["IVC", "OVC", "BVLTCASE", "BVLBCASE", "DIVPSRT", "DIVPSRB", "HFSPSRT", "HFSPSRB"] # TODO: TEMPORARY while debugging
118 for passive_name in passive_names:
119 # The data structure looks like this:
120 # passives(["BVLBCASE", "dof", "constant_current_density"])
121 # passives(["BVLTCASE", "dof", "constant_current_density"])
122 # ...
123 # passives(["IVC", "dof", "eig_01"])
124 # passives(["IVC", "dof", "eig_02"])
125 # ...
126 n_passive_dofs += len(passives.keys([passive_name, "dof"]))
128 # The data structure looks like this:
129 # passives(["BVLBCASE", "regularisations"]) # shape = [n_regularisations, n_dof]
130 # passives(["BVLTCASE", "regularisations"])
131 # ...
132 # passives(["IVC", "regularisations"])
133 # ...
134 [n_regularisations_local, _] = passives.get_array2([passive_name, "regularisations"]).shape
135 n_regularisations += n_regularisations_local
137 # Total number of constraints
138 n_constraints = n_flux_loops_to_include + n_bp_probes_to_include + n_rogowski_coils_to_include + n_regularisations
140 # Store the number of constraints and regularisations
141 results["PRESHOT"]["N_F_LOOPS"] = np.int32(n_flux_loops_to_include)
142 results["PRESHOT"]["N_BP_PROBES"] = np.int32(n_bp_probes_to_include)
143 results["PRESHOT"]["N_ROG_COILS"] = np.int32(n_rogowski_coils_to_include)
144 results["PRESHOT"]["N_MEAS"] = np.int32(n_constraints)
145 results["PRESHOT"]["N_REG"] = np.int32(n_regularisations)
147 # Store number of plasma degrees of freedom
148 n_p_prime = len(plasma.get(ep.time_slice[0].source_functions.p_prime.coefficients))
149 n_ff_prime = len(plasma.get(ep.time_slice[0].source_functions.ff_prime.coefficients))
150 n_delta_z = 1
151 n_plasma_dof = n_p_prime + n_ff_prime + n_delta_z
152 results["PRESHOT"]["N_PLS"] = np.int32(n_plasma_dof)
153 # Total number of degrees of freedom
154 n_coef = n_passive_dofs + n_plasma_dof
155 results["PRESHOT"]["N_COEF"] = np.int32(n_coef)
157 # Greens for "sensors-coils"
158 g_measured_coil = np.zeros((n_constraints, n_psu))
160 # Weights for the constraints
161 constraints_weight = np.zeros(n_constraints)
163 # Lis of constraint names
164 constraint_names = []
166 # Greens between the measurements and the degrees of freedom
167 # Note: the plasma's dof's are calculated during real-time and are set to zero.
168 # So g_dof_meas[:, 0 : n_plasma_dof] will be zero.
169 g_dof_meas = np.zeros((n_constraints, n_coef))
171 # Collect the greens for "grid-measurements"
172 g_grid_meas = np.zeros((n_constraints, n_r * n_z))
174 # flux_loops are the first set of constraints
175 i_constraint = 0
176 for i_flux_loop, floop_name in enumerate(flux_loops.keys()):
177 if flux_loops_to_include[i_flux_loop]:
178 # Add flux loop name in PCS (Plasma Control System) format
179 floop_name_pcs = floop_name.replace("L", "PSI_FLOOP_")
180 constraint_names.append(floop_name_pcs)
181 # Add the Greens between measurements and coils
182 for i_psu, power_supply in enumerate(rtgsfit_psus):
183 coil_names = power_supply["coils"]
184 for coil_name in coil_names:
185 g_measured_coil[i_constraint, i_psu] += flux_loops.get_f64([floop_name, "greens", "pf", coil_name])
186 # Add the weight
187 constraints_weight[i_constraint] = (
188 2 * np.pi * flux_loops.get_f64([floop_name, "fit_settings", "weight"]) / flux_loops.get_f64([floop_name, "fit_settings", "expected_value"])
189 )
190 # Add the Greens between measurements and degrees of freedom
191 i_vessel_dof = 0
192 for passive_name in passive_names:
193 dof_names = passives.keys([passive_name, "dof"])
194 for dof_name in dof_names:
195 g_dof_meas[i_constraint, n_plasma_dof + i_vessel_dof] = flux_loops.get_f64([floop_name, "greens", "passives", passive_name, dof_name])
196 i_vessel_dof += 1
197 # Greens between the sensors and the plasma grid
198 g_grid_meas[i_constraint, :] = flux_loops.get_array1([floop_name, "greens", "plasma"])
199 # Set-up index for next sensor
200 i_constraint += 1
202 # bp_probes are the second set of constraints
203 for i_bp_probe, bp_name in enumerate(bp_probes.keys()):
204 if bp_probes_to_include[i_bp_probe]:
205 # Add flux loop name in PCS (Plasma Control System) format
206 bp_name_pcs = bp_name.replace("P", "B_BPPROBE_")
207 constraint_names.append(bp_name_pcs)
208 # Add the Greens between measurements and coils
209 for i_psu, power_supply in enumerate(rtgsfit_psus):
210 coil_names = power_supply["coils"]
211 for coil_name in coil_names:
212 g_measured_coil[i_constraint, i_psu] += bp_probes.get_f64([bp_name, "greens", "pf", coil_name])
213 # Add the weight
214 constraints_weight[i_constraint] = bp_probes.get_f64([bp_name, "fit_settings", "weight"]) / bp_probes.get_f64(
215 [bp_name, "fit_settings", "expected_value"]
216 )
217 # Add the Greens between measurements and degrees of freedom
218 i_vessel_dof = 0
219 for passive_name in passive_names:
220 dof_names = passives.keys([passive_name, "dof"])
221 for dof_name in dof_names:
222 g_dof_meas[i_constraint, n_plasma_dof + i_vessel_dof] = bp_probes.get_f64([bp_name, "greens", "passives", passive_name, dof_name])
223 i_vessel_dof += 1
224 # Greens between the sensors and the plasma grid
225 g_grid_meas[i_constraint, :] = bp_probes.get_array1([bp_name, "greens", "plasma"])
226 # Set-up index for next sensor
227 i_constraint += 1
229 rogowski_coils_names_rtgsfit_order = [
230 "INIVC000",
231 "BVLT",
232 "BVLB",
233 "GASBFLT",
234 "GASBFLB",
235 "HFSPSRT",
236 "HFSPSRB",
237 "DIVPSRT",
238 "DIVPSRB",
239 ]
241 # rogowski_coils are the third set of constraints
242 # for i_rogowski_coil, rogowski_coil_name in enumerate(rogowski_coils.keys()):
243 for i_rogowski_coil, rogowski_coil_name in enumerate(rogowski_coils_names_rtgsfit_order):
244 # if rogowski_coils_to_include[i_rogowski_coil]:
245 if rogowski_coils.get_vec_bool([rogowski_coil_name, "fit_settings", "include"]):
246 # Add Rogowski coil name in PCS (Plasma Control System) format
247 rogowski_coil_name_pcs = f"I_ROG_{rogowski_coil_name}"
248 constraint_names.append(rogowski_coil_name_pcs)
249 # Add the Greens between measurements and coils
250 for i_psu, power_supply in enumerate(rtgsfit_psus):
251 coil_names = power_supply["coils"]
252 for coil_name in coil_names:
253 g_measured_coil[i_constraint, i_psu] += rogowski_coils.get_f64([rogowski_coil_name, "greens", "pf", coil_name])
254 # Add the weight
255 constraints_weight[i_constraint] = rogowski_coils.get_f64([rogowski_coil_name, "fit_settings", "weight"]) / rogowski_coils.get_f64(
256 [rogowski_coil_name, "fit_settings", "expected_value"]
257 )
258 # Add the Greens between measurements and degrees of freedom
259 i_vessel_dof = 0
260 for passive_name in passive_names:
261 dof_names = passives.keys([passive_name, "dof"])
262 for dof_name in dof_names:
263 g_dof_meas[i_constraint, n_plasma_dof + i_vessel_dof] = rogowski_coils.get_f64(
264 [rogowski_coil_name, "greens", "passives", passive_name, dof_name]
265 )
266 i_vessel_dof += 1
267 # Greens between the sensors and the plasma grid
268 g_grid_meas[i_constraint, :] = rogowski_coils.get_array1([rogowski_coil_name, "greens", "plasma"])
269 # Set-up index for next sensor
270 i_constraint += 1
272 # "passive regularisations" are the fourth set of constraints
273 # Loop over all passives and add regularisations if they exist
274 i_dof_start = n_plasma_dof
275 for passive_name in passive_names:
276 # Find the number of degrees of freedom for this passive
277 n_dof_local = len(passives.keys([passive_name, "dof"]))
279 passive_regularisation_local = passives.get_array2([passive_name, "regularisations"])
280 [n_reg_local, _] = passive_regularisation_local.shape
282 i_dof_end = i_dof_start + n_dof_local
284 for i_reg in range(n_reg_local):
285 g_dof_meas[i_constraint, i_dof_start:i_dof_end] = passive_regularisation_local[i_reg, :]
286 regularisation_scaling = 0.001 * np.pi
287 constraints_weight[i_constraint] = regularisation_scaling * passives.get_array1([passive_name, "regularisations_weight"])[i_reg]
289 # Set-up index for next sensor / constraint
290 i_constraint += 1
292 # Set-up index for next passive
293 i_dof_start += n_dof_local
295 # Store in MDSplus
296 results["PRESHOT"]["GREENS"]["MEAS_COIL"] = g_measured_coil.flatten()
297 results["PRESHOT"]["WEIGHT"] = constraints_weight
298 g_dof_meas_weight = np.dot(np.diag(constraints_weight), g_dof_meas).T
299 results["PRESHOT"]["GREENS"]["COEF_MEAS_W"] = g_dof_meas_weight.flatten() # .reshape((n_coef, n_constraints))
300 g_grid_meas_weight = np.dot(np.diag(constraints_weight), g_grid_meas).T * d_r * d_z
301 results["PRESHOT"]["GREENS"]["GRID_MEAS_W"] = g_grid_meas_weight.flatten()
303 # Collect the greens for "grid-vessel"
304 g_grid_vessel = np.zeros((n_z * n_r, n_passive_dofs))
305 # Loop over all passives
306 # `greens/pf_passive` is in the same order as `passives.keys()`, and each passive's `dof` in the
307 # same order as its own dof keys, because they were built in that order
308 i_dof = 0
309 for i_passive, passive_name in enumerate(passive_names):
310 current_distribution_dof_names = passives.keys([passive_name, "dof"])
311 # `current_distribution_dof_names` can be "constant_current_density", "eig_01", "eig_02", etc.
312 for i_passive_dof, _current_distribution_dof_name in enumerate(current_distribution_dof_names):
313 g_grid_vessel[:, i_dof] = plasma.get(ep.greens.pf_passive[i_passive].dof[i_passive_dof].psi)
314 i_dof += 1
315 # Store in MDSplus
316 results["PRESHOT"]["GREENS"]["GRID_VESSEL"] = g_grid_vessel.flatten() # .reshape((n_z * n_r, n_passive_dofs))
317 results["PRESHOT"]["N_VESS"] = np.int32(n_passive_dofs)
319 # Store some settings
320 rtgsfit_code_settings = gsfit_controller.settings["RTGSFIT_code_settings.json"]
321 results["PRESHOT"]["N_XPT_MAX"] = np.int32(rtgsfit_code_settings["n_xpt_max"])
322 results["PRESHOT"]["N_LCFS_MAX"] = np.int32(rtgsfit_code_settings["n_lcfs_max"])
323 results["PRESHOT"]["N_INTRP"] = np.int32(rtgsfit_code_settings["n_intrp"])
324 results["PRESHOT"]["THRESH"] = rtgsfit_code_settings["thresh"]
325 results["PRESHOT"]["FRAC"] = rtgsfit_code_settings["frac"]
327 # Add initial conditions
328 flux_norm = gsfit_controller.settings["rtgsfit_initial_conditions.json"]["flux_norm"]
329 mask = gsfit_controller.settings["rtgsfit_initial_conditions.json"]["mask"]
330 psi_total = gsfit_controller.settings["rtgsfit_initial_conditions.json"]["psi_total"]
331 results["PRESHOT"]["INITIAL_COND"]["FLUX_NORM"] = np.array(flux_norm).astype(np.float64)
332 results["PRESHOT"]["INITIAL_COND"]["MASK"] = np.array(mask).astype(np.int32)
333 results["PRESHOT"]["INITIAL_COND"]["PSI_TOTAL"] = np.array(psi_total).astype(np.float64)
335 # Vessel
336 # `unit(0)` is the vacuum vessel contour, which is the region the plasma is allowed to occupy
337 vessel_r = wall.get(wp.description_2d[0].limiter.unit[0].outline.r)
338 vessel_z = wall.get(wp.description_2d[0].limiter.unit[0].outline.z)
339 vessel_polygon = shapely.geometry.Polygon(np.column_stack((vessel_r, vessel_z)))
341 # Test if grid-points are inside the vessel polygon
342 grid_points = []
343 for i_z in range(n_z):
344 for i_r in range(n_r):
345 grid_points.append(shapely.geometry.Point(r[i_r], z[i_z]))
346 mask_lim = vessel_polygon.contains(grid_points)
348 # # Ensure that the (R,Z) grid cell nearest to the vessel is included in the mask
349 # for r_v, z_v in zip(vessel_r, vessel_z):
350 # # Find the nearest grid-point to the vessel point
351 # i_r_nearest = np.argmin(np.abs(r - r_v))
352 # i_z_nearest = np.argmin(np.abs(z - z_v))
353 # mask_lim[i_z_nearest * n_r + i_r_nearest] = True
355 # results["PRESHOT"]["MASK_LIM"] = mask_lim.astype(np.int32)
357 def compute_limit_idx_and_weights(
358 r: npt.NDArray[np.float64],
359 z: npt.NDArray[np.float64],
360 lim_r: npt.NDArray[np.float64],
361 lim_z: npt.NDArray[np.float64],
362 n_intrp: int,
363 ) -> tuple[npt.NDArray[np.int32], npt.NDArray[np.float64]]:
364 n_r: int = len(r)
365 n_lim: int = len(lim_r)
366 limit_idx = np.zeros(n_lim * n_intrp, dtype=int)
367 limit_w = np.zeros(n_lim * n_intrp, dtype=float)
369 for i, (lr, lz) in enumerate(zip(lim_r, lim_z)):
370 if lr < r[0] or lr > r[-1] or lz < z[0] or lz > z[-1]:
371 raise ValueError(f"Limiter point ({lr}, {lz}) is out of bounds of the grid.")
372 r_idx = np.searchsorted(r, lr) - 1
373 z_idx = np.searchsorted(z, lz) - 1
375 limit_idx[n_intrp * i + 0] = n_r * z_idx + r_idx # (r_idx, z_idx)
376 limit_idx[n_intrp * i + 1] = n_r * z_idx + r_idx + 1 # (r_idx + 1, z_idx)
377 limit_idx[n_intrp * i + 2] = n_r * (z_idx + 1) + r_idx # (r_idx, z_idx + 1)
378 limit_idx[n_intrp * i + 3] = n_r * (z_idx + 1) + r_idx + 1 # (r_idx + 1, z_idx + 1)
380 r0, r1 = r[r_idx], r[r_idx + 1]
381 z0, z1 = z[z_idx], z[z_idx + 1]
382 dr = r1 - r0
383 dz = z1 - z0
384 limit_w[n_intrp * i + 0] = (r1 - lr) * (z1 - lz) / (dr * dz) # (r_idx, z_idx)
385 limit_w[n_intrp * i + 1] = (lr - r0) * (z1 - lz) / (dr * dz) # (r_idx + 1, z_idx)
386 limit_w[n_intrp * i + 2] = (r1 - lr) * (lz - z0) / (dr * dz) # (r_idx, z_idx + 1)
387 limit_w[n_intrp * i + 3] = (lr - r0) * (lz - z0) / (dr * dz) # (r_idx + 1, z_idx + 1)
389 return limit_idx, limit_w
391 n_intrp: int = int(rtgsfit_code_settings["n_intrp"])
392 # Every limiter unit contributes candidate limit points, so they are concatenated
393 n_limiter_unit = len(wall.get(wp.description_2d[0].limiter.unit[:].name))
394 lim_r = np.concatenate([wall.get(wp.description_2d[0].limiter.unit[i].outline.r) for i in range(n_limiter_unit)])
395 lim_z = np.concatenate([wall.get(wp.description_2d[0].limiter.unit[i].outline.z) for i in range(n_limiter_unit)])
396 # Remove indices where |lim_z| > 0.7 m
397 # lim_r = lim_r[np.abs(lim_z) < 0.7]
398 # lim_z = lim_z[np.abs(lim_z) < 0.7]
399 n_lim = len(lim_r)
400 limit_idx, limit_w = compute_limit_idx_and_weights(r, z, lim_r, lim_z, n_intrp)
401 results["PRESHOT"]["N_LIMIT"] = np.int32(n_lim)
402 results["PRESHOT"]["LIMIT_IDX"] = limit_idx.astype(np.int32)
403 results["PRESHOT"]["LIMIT_W"] = limit_w.astype(np.float64)
404 results["PRESHOT"]["LIMIT_R"] = lim_r.astype(np.float64)
405 results["PRESHOT"]["LIMIT_Z"] = lim_z.astype(np.float64)
407 # If grid point is within sqrt(d_r^2 + d_z^2) of
408 # a limiter point then include it in the mask_lim
409 d_rz = np.sqrt(d_r**2 + d_z**2) / 2
410 for i, (lr, lz) in enumerate(zip(lim_r, lim_z)):
411 points_near_lim = np.sqrt((lr - r_flat) ** 2 + (lz - z_flat) ** 2) < d_rz
412 mask_lim = np.logical_or(mask_lim, points_near_lim)
414 results["PRESHOT"]["MASK_LIM"] = mask_lim.astype(np.int32)
416 r_ltrb = np.concatenate(
417 (
418 [r[0]], # (bottom, left)
419 r[0] * np.ones(len(z[1:-1])), # traverse (bottom, left) to (top, left) (excluding corners)
420 [r[0]], # (top, left)
421 r[1:-1], # traverse (top, left) to (top, right) (excluding corners)
422 [r[-1]], # (top, right)
423 r[-1] * np.ones(len(z[1:-1])), # traverse (top, right) to (bottom, right) (excluding corners)
424 [r[-1]], # (bottom, right)
425 np.flip(r[1:-1]), # traverse (bottom, right) to (bottom, left) (excluding corners)
426 )
427 )
428 inv_r_ltrb_mu0 = 1.0 / (r_ltrb * mu_0)
429 results["PRESHOT"]["INV_R_L_MU0"] = inv_r_ltrb_mu0.astype(np.float64)
431 lower_band, upper_band, perm_idx = compute_lup_bands(r, z)
432 results["PRESHOT"]["LOWER_BAND"] = lower_band.astype(np.float64)
433 results["PRESHOT"]["UPPER_BAND"] = upper_band.astype(np.float64)
434 results["PRESHOT"]["PERM_IDX"] = perm_idx.astype(np.int32)
436 # Greens with the boundary points
437 g_ltrb = greens_with_boundary_points(plasma)
438 results["PRESHOT"]["GREENS"]["LTRB"] = g_ltrb
440 # Sometimes important sensors are broken, mainly Rogowski coils which measure current in
441 # passive plates. As a work around we can use the nearby flux loop voltage (V = I * R) as
442 # a replacement. This section replaces the damaged "bad" sensors with working "good" sensors.
443 sensors_pcs_should_read = copy.deepcopy(constraint_names)
444 if "sensor_replacement.json" in gsfit_controller.settings:
445 sensor_replacements = gsfit_controller.settings["sensor_replacement.json"]
446 # sensor_replacement is a dictonary containing the sensors which should be replaced, the data structure is:
447 # {
448 # "sensor_name_to_replace": {
449 # "replacements": ["sensor_name_to_use_01", "sensor_name_to_use_02", ...],
450 # "coefficients": [1.0, 0.5, ...]
451 # },
452 # ...
453 # }
455 # Remove the sensors that are being replaced
456 for sensor_replacement_name in sensor_replacements.keys():
457 if sensor_replacement_name in sensors_pcs_should_read:
458 sensors_pcs_should_read.remove(sensor_replacement_name)
460 # Add the sensors that are replacing the removed sensors
461 for sensor_replacement_item in sensor_replacements.values():
462 sensor_replacement_names = sensor_replacement_item["replacements"]
463 for sensor_replacement_name in sensor_replacement_names:
464 # Don't double add sensor names
465 if sensor_replacement_name not in sensors_pcs_should_read:
466 sensors_pcs_should_read.append(sensor_replacement_name)
468 # Construct a matrix where:
469 # `sensors_rtgsfit_wants = sensor_replacement_matrix * sensors_pcs_should_read`
470 sensor_replacement_matrix = np.zeros((len(constraint_names), len(sensors_pcs_should_read)), dtype=np.float64)
471 for i_constraint, constraint_name in enumerate(constraint_names):
472 # If the sensor is being replaced, only use the replacement coefficients
473 if constraint_name in sensor_replacements:
474 sensor_replacement = sensor_replacements[constraint_name]
475 sensor_replacement_names = sensor_replacement["replacements"]
476 sensor_replacement_coefficients = sensor_replacement["coefficients"]
477 n_replacements = len(sensor_replacement_names)
478 for i_replacement in range(n_replacements):
479 replacement_name = sensor_replacement_names[i_replacement]
480 i_pcs_sensor = sensors_pcs_should_read.index(replacement_name)
481 sensor_replacement_matrix[i_constraint, i_pcs_sensor] = sensor_replacement_coefficients[i_replacement]
482 elif constraint_name in sensors_pcs_should_read:
483 # Sensor is not being replaced, set identity
484 i_pcs_sensor = sensors_pcs_should_read.index(constraint_name)
485 sensor_replacement_matrix[i_constraint, i_pcs_sensor] = 1.0
487 else:
488 # No sensor replacement, so the matrix is just the identity matrix
489 sensor_replacement_matrix = np.eye(len(constraint_names), len(sensors_pcs_should_read), dtype=np.float64)
491 results["PRESHOT"]["SENS_NAMES"] = np.array(sensors_pcs_should_read)
492 results["PRESHOT"]["SENS_REP_MAT"] = sensor_replacement_matrix.flatten()
493 results["PRESHOT"]["N_SENS_PCS"] = np.int32(len(sensors_pcs_should_read))
495 # Save IVC geometry data
496 results["PASSIVES"]["IVC"]["GEOMETRY"]["R"] = passives.get_array1(["IVC", "geometry", "r"])
497 results["PASSIVES"]["IVC"]["GEOMETRY"]["Z"] = passives.get_array1(["IVC", "geometry", "z"])
498 results["PASSIVES"]["IVC"]["GEOMETRY"]["D_R"] = passives.get_array1(["IVC", "geometry", "d_r"])
499 results["PASSIVES"]["IVC"]["GEOMETRY"]["D_Z"] = passives.get_array1(["IVC", "geometry", "d_z"])
500 results["PASSIVES"]["IVC"]["GEOMETRY"]["AREA"] = results["PASSIVES"]["IVC"]["GEOMETRY"]["D_R"] * results["PASSIVES"]["IVC"]["GEOMETRY"]["D_Z"]
501 n_eigs = gsfit_controller.settings["passive_dof_regularisation.json"]["IVC"]["n_dof"]
502 n_segs = len(passives.get_array1(["IVC", "dof", f"eig_01", "current_distribution"]))
503 current_dofs = np.zeros((n_eigs, n_segs))
504 for eig_num in range(n_eigs):
505 current_dofs[eig_num, :] = passives.get_array1(["IVC", "dof", f"eig_{eig_num + 1:02d}", "current_distribution"])
506 results["PASSIVES"]["IVC"]["GEOMETRY"]["CURRENT_DOFS"] = current_dofs
508 # Create coef names list and save it
509 coef_names = ["pls0", "pls1", "pls2"]
510 for passive_name in passives.keys():
511 dof_names = passives.keys([passive_name, "dof"])
512 for dof_name in dof_names:
513 if dof_name == "constant_current_density":
514 coef_names.append(passive_name)
515 elif dof_name.startswith("eig_"):
516 coef_names.append(dof_name)
517 else:
518 raise ValueError(f"Unknown DoF name: {dof_name}")
519 results["PRESHOT"]["COEF_NAMES"] = np.array(coef_names)
521 # create meas_names list and save it
522 meas_names = []
523 for i_flux_loop, floop_name in enumerate(flux_loops.keys()):
524 if flux_loops_to_include[i_flux_loop]:
525 meas_names.append(floop_name)
526 for i_bp_probe, bp_name in enumerate(bp_probes.keys()):
527 if bp_probes_to_include[i_bp_probe]:
528 meas_names.append(bp_name)
529 for rog_name in rogowski_coils_names_rtgsfit_order:
530 meas_names.append(rog_name)
531 for passive_name in passive_names:
532 passive_regularisation_local = passives.get_array2([passive_name, "regularisations"])
533 [n_reg_local, _] = passive_regularisation_local.shape
534 for i_reg in range(n_reg_local):
535 meas_names.append(f"{passive_name}_reg_{i_reg}")
536 results["PRESHOT"]["MEAS_NAMES"] = np.array(meas_names)