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