Coverage for python/gsfit/database_writers/tokamak_energy_mdsplus_new/map_results_to_database.py: 48%
138 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-14 15:25 +0000
1import typing
2from typing import TYPE_CHECKING
4import numpy as np
5import numpy.typing as npt
6from gsfit_rs.imas import equilibrium_paths as ep
7from scipy.constants import mu_0
9if TYPE_CHECKING:
10 from ...gsfit import Gsfit
11 from . import DatabaseWriterTokamakEnergyMDSplusNew
14# Pairs of (MDSplus node path, IMAS path) for everything which comes out of the `equilibrium` IDS
15# one value per time-slice. An IMAS path is a value: it holds no data, so the table below can be
16# built once, at import, and read against any IDS.
17#
18# `time_slice[:]` gathers over every time-slice, giving the `[n_time, ...]` arrays MDSplus expects.
19_TIME_SERIES_PATH_PAIRS: list[tuple[tuple[str, ...], typing.Any]] = [
20 # Plasma boundary
21 (("BOUNDARY", "GEO_AXIS", "R"), ep.time_slice[:].boundary.geometric_axis.r),
22 (("BOUNDARY", "GEO_AXIS", "Z"), ep.time_slice[:].boundary.geometric_axis.z),
23 (("BOUNDARY", "MINOR_RADIUS"), ep.time_slice[:].boundary.minor_radius),
24 (("BOUNDARY", "ELONGATION"), ep.time_slice[:].boundary.elongation),
25 (("BOUNDARY", "TRIANG"), ep.time_slice[:].boundary.triangularity),
26 (("BOUNDARY", "TRIANG_L"), ep.time_slice[:].boundary.triangularity_lower),
27 (("BOUNDARY", "TRIANG_U"), ep.time_slice[:].boundary.triangularity_upper),
28 (("BOUNDARY", "SQUARE_L_I"), ep.time_slice[:].boundary.squareness_lower_inner),
29 (("BOUNDARY", "SQUARE_L_O"), ep.time_slice[:].boundary.squareness_lower_outer),
30 (("BOUNDARY", "SQUARE_U_I"), ep.time_slice[:].boundary.squareness_upper_inner),
31 (("BOUNDARY", "SQUARE_U_O"), ep.time_slice[:].boundary.squareness_upper_outer),
32 (("BOUNDARY", "PSI"), ep.time_slice[:].boundary.psi),
33 (("BOUNDARY", "OUTLINE", "R"), ep.time_slice[:].boundary.outline.r),
34 (("BOUNDARY", "OUTLINE", "Z"), ep.time_slice[:].boundary.outline.z),
35 (("BOUNDARY", "BOUNDING", "R"), ep.time_slice[:].boundary.bounding.r),
36 (("BOUNDARY", "BOUNDING", "Z"), ep.time_slice[:].boundary.bounding.z),
37 # Convergence
38 (("CONVERGENCE", "GS_ERROR"), ep.time_slice[:].convergence.grad_shafranov_deviation_value),
39 (("CONVERGENCE", "ITERATIONS_N"), ep.time_slice[:].convergence.iterations_n),
40 # Global
41 (("GLOBAL", "CURRENT_CENT", "R"), ep.time_slice[:].global_quantities.current_centre.r),
42 (("GLOBAL", "CURRENT_CENT", "VELOCITY_Z"), ep.time_slice[:].global_quantities.current_centre.velocity_z),
43 (("GLOBAL", "CURRENT_CENT", "Z"), ep.time_slice[:].global_quantities.current_centre.z),
44 (("GLOBAL", "MAG_AXIS", "R"), ep.time_slice[:].global_quantities.magnetic_axis.r),
45 (("GLOBAL", "MAG_AXIS", "B_FIELD_PHI"), ep.time_slice[:].global_quantities.magnetic_axis.b_field_phi),
46 (("GLOBAL", "MAG_AXIS", "Z"), ep.time_slice[:].global_quantities.magnetic_axis.z),
47 (("GLOBAL", "Q_MIN", "VALUE"), ep.time_slice[:].global_quantities.q_min.value),
48 (("GLOBAL", "Q_MIN", "RHO_TOR_NORM"), ep.time_slice[:].global_quantities.q_min.rho_tor_norm),
49 (("GLOBAL", "Q_MIN", "PSI_NORM"), ep.time_slice[:].global_quantities.q_min.psi_norm),
50 (("GLOBAL", "Q_MIN", "PSI"), ep.time_slice[:].global_quantities.q_min.psi),
51 (("GLOBAL", "AREA"), ep.time_slice[:].global_quantities.area),
52 (("GLOBAL", "BETA_N"), ep.time_slice[:].global_quantities.beta_tor_norm),
53 (("GLOBAL", "BETA_P_1"), ep.time_slice[:].global_quantities.beta_pol_1),
54 (("GLOBAL", "BETA_P_2"), ep.time_slice[:].global_quantities.beta_pol_2),
55 (("GLOBAL", "BETA_P_3"), ep.time_slice[:].global_quantities.beta_pol_3),
56 (("GLOBAL", "BT_VAC_RGEO"), ep.time_slice[:].global_quantities.bt_vac_at_r_geo),
57 (("GLOBAL", "DELTA_R_SEP"), ep.time_slice[:].global_quantities.delta_r_sep),
58 (("GLOBAL", "DELTA_Z"), ep.time_slice[:].convergence.delta_z),
59 (("GLOBAL", "ENERGY_MHD"), ep.time_slice[:].global_quantities.energy_mhd),
60 (("GLOBAL", "FX"), ep.time_slice[:].global_quantities.f_x),
61 (("GLOBAL", "IP"), ep.time_slice[:].global_quantities.ip),
62 (("GLOBAL", "LI_1"), ep.time_slice[:].global_quantities.li_1),
63 (("GLOBAL", "LI_2"), ep.time_slice[:].global_quantities.li_2),
64 (("GLOBAL", "LI_3"), ep.time_slice[:].global_quantities.li_3),
65 (("GLOBAL", "LENGTH_POL"), ep.time_slice[:].global_quantities.length_pol),
66 (("GLOBAL", "PHI_DIA"), ep.time_slice[:].constraints.diamagnetic_flux.reconstructed),
67 (("GLOBAL", "PSI_MAG_AXIS"), ep.time_slice[:].global_quantities.psi_magnetic_axis),
68 (("GLOBAL", "Q_AXIS"), ep.time_slice[:].global_quantities.q_axis),
69 (("GLOBAL", "Q_95"), ep.time_slice[:].global_quantities.q_95),
70 (("GLOBAL", "SURFACE"), ep.time_slice[:].global_quantities.surface),
71 (("GLOBAL", "V_LOOP"), ep.time_slice[:].global_quantities.v_loop),
72 (("GLOBAL", "VOLUME"), ep.time_slice[:].global_quantities.volume),
73 (("GLOBAL", "XPT_DIVERTED"), ep.time_slice[:].boundary.type),
74 # Profiles_1d, on the psi_norm grid
75 (("PROFILES_1D", "PSI_NORM", "GEO_AXIS", "R"), ep.time_slice[:].profiles_1d.geometric_axis.r),
76 (("PROFILES_1D", "PSI_NORM", "GEO_AXIS", "Z"), ep.time_slice[:].profiles_1d.geometric_axis.z),
77 (("PROFILES_1D", "PSI_NORM", "AREA"), ep.time_slice[:].profiles_1d.area),
78 (("PROFILES_1D", "PSI_NORM", "DAREA_DPSI"), ep.time_slice[:].profiles_1d.darea_dpsi),
79 (("PROFILES_1D", "PSI_NORM", "DAREA_DRHO_T"), ep.time_slice[:].profiles_1d.darea_drho_tor),
80 (("PROFILES_1D", "PSI_NORM", "ELONGATION"), ep.time_slice[:].profiles_1d.elongation),
81 (("PROFILES_1D", "PSI_NORM", "F"), ep.time_slice[:].profiles_1d.f),
82 (("PROFILES_1D", "PSI_NORM", "FF_PRIME"), ep.time_slice[:].profiles_1d.f_df_dpsi),
83 (("PROFILES_1D", "PSI_NORM", "FLUX_TOR"), ep.time_slice[:].profiles_1d.phi),
84 (("PROFILES_1D", "PSI_NORM", "MAG_SHEAR"), ep.time_slice[:].profiles_1d.magnetic_shear),
85 (("PROFILES_1D", "PSI_NORM", "P_PRIME"), ep.time_slice[:].profiles_1d.dpressure_dpsi),
86 (("PROFILES_1D", "PSI_NORM", "PRESSURE"), ep.time_slice[:].profiles_1d.pressure),
87 (("PROFILES_1D", "PSI_NORM", "PSI"), ep.time_slice[:].profiles_1d.psi),
88 (("PROFILES_1D", "PSI_NORM", "PHI"), ep.time_slice[:].profiles_1d.phi),
89 (("PROFILES_1D", "PSI_NORM", "J_PHI"), ep.time_slice[:].profiles_1d.j_phi),
90 (("PROFILES_1D", "PSI_NORM", "J_PARALLEL"), ep.time_slice[:].profiles_1d.j_parallel),
91 (("PROFILES_1D", "PSI_NORM", "Q"), ep.time_slice[:].profiles_1d.q),
92 (("PROFILES_1D", "PSI_NORM", "RHO_POL"), ep.time_slice[:].profiles_1d.rho_pol),
93 (("PROFILES_1D", "PSI_NORM", "R_INBOARD"), ep.time_slice[:].profiles_1d.r_inboard),
94 (("PROFILES_1D", "PSI_NORM", "R_OUTBOARD"), ep.time_slice[:].profiles_1d.r_outboard),
95 (("PROFILES_1D", "PSI_NORM", "R_INBOARD_0"), ep.time_slice[:].profiles_1d.r_inboard_z_0),
96 (("PROFILES_1D", "PSI_NORM", "R_OUTBOARD_0"), ep.time_slice[:].profiles_1d.r_outboard_z_0),
97 (("PROFILES_1D", "PSI_NORM", "RHO_TOR"), ep.time_slice[:].profiles_1d.rho_tor),
98 (("PROFILES_1D", "PSI_NORM", "RHO_TOR_NORM"), ep.time_slice[:].profiles_1d.rho_tor_norm),
99 (("PROFILES_1D", "PSI_NORM", "SQUARE_L_I"), ep.time_slice[:].profiles_1d.squareness_lower_inner),
100 (("PROFILES_1D", "PSI_NORM", "SQUARE_L_O"), ep.time_slice[:].profiles_1d.squareness_lower_outer),
101 (("PROFILES_1D", "PSI_NORM", "SQUARE_U_I"), ep.time_slice[:].profiles_1d.squareness_upper_inner),
102 (("PROFILES_1D", "PSI_NORM", "SQUARE_U_O"), ep.time_slice[:].profiles_1d.squareness_upper_outer),
103 (("PROFILES_1D", "PSI_NORM", "TRIANG"), ep.time_slice[:].profiles_1d.triangularity),
104 (("PROFILES_1D", "PSI_NORM", "TRIANG_L"), ep.time_slice[:].profiles_1d.triangularity_lower),
105 (("PROFILES_1D", "PSI_NORM", "TRIANG_U"), ep.time_slice[:].profiles_1d.triangularity_upper),
106 (("PROFILES_1D", "PSI_NORM", "VOL"), ep.time_slice[:].profiles_1d.volume),
107 (("PROFILES_1D", "PSI_NORM", "VOL_PRIME"), ep.time_slice[:].profiles_1d.dvolume_dpsi),
108 (("PROFILES_1D", "PSI_NORM", "RHO_VOL_NORM"), ep.time_slice[:].profiles_1d.rho_volume_norm),
109 (("PROFILES_1D", "PSI_NORM", "DVOL_DPSI"), ep.time_slice[:].profiles_1d.dvolume_dpsi),
110 (("PROFILES_1D", "PSI_NORM", "DVOL_DRHO_T"), ep.time_slice[:].profiles_1d.dvolume_drho_tor),
111 (("PROFILES_1D", "PSI_NORM", "SURFACE"), ep.time_slice[:].profiles_1d.surface),
112 (("PROFILES_1D", "PSI_NORM", "GM1"), ep.time_slice[:].profiles_1d.gm1),
113 (("PROFILES_1D", "PSI_NORM", "GM2"), ep.time_slice[:].profiles_1d.gm2),
114 (("PROFILES_1D", "PSI_NORM", "GM3"), ep.time_slice[:].profiles_1d.gm3),
115 (("PROFILES_1D", "PSI_NORM", "GM4"), ep.time_slice[:].profiles_1d.gm4),
116 (("PROFILES_1D", "PSI_NORM", "GM5"), ep.time_slice[:].profiles_1d.gm5),
117 (("PROFILES_1D", "PSI_NORM", "GM6"), ep.time_slice[:].profiles_1d.gm6),
118 (("PROFILES_1D", "PSI_NORM", "GM7"), ep.time_slice[:].profiles_1d.gm7),
119 (("PROFILES_1D", "PSI_NORM", "GM8"), ep.time_slice[:].profiles_1d.gm8),
120 (("PROFILES_1D", "PSI_NORM", "GM9"), ep.time_slice[:].profiles_1d.gm9),
121 (("PROFILES_1D", "PSI_NORM", "B_FIELD_AV"), ep.time_slice[:].profiles_1d.b_field_average),
122 (("PROFILES_1D", "PSI_NORM", "B_FIELD_MIN"), ep.time_slice[:].profiles_1d.b_field_min),
123 (("PROFILES_1D", "PSI_NORM", "B_FIELD_MAX"), ep.time_slice[:].profiles_1d.b_field_max),
124 # Mid-plane profiles
125 (("PROFILES_1D", "R_MIDPLANE", "P_PRIME"), ep.time_slice[:].profiles_1d_r_midplane.dpressure_dpsi),
126 (("PROFILES_1D", "R_MIDPLANE", "F"), ep.time_slice[:].profiles_1d_r_midplane.f),
127 (("PROFILES_1D", "R_MIDPLANE", "FF_PRIME"), ep.time_slice[:].profiles_1d_r_midplane.f_df_dpsi),
128 (("PROFILES_1D", "R_MIDPLANE", "J_PHI"), ep.time_slice[:].profiles_1d_r_midplane.j_phi),
129 (("PROFILES_1D", "R_MIDPLANE", "PRESSURE"), ep.time_slice[:].profiles_1d_r_midplane.pressure),
130 (("PROFILES_1D", "R_MIDPLANE", "Q"), ep.time_slice[:].profiles_1d_r_midplane.q),
131 # SOL
132 (("PROFILES_1D", "R_MIDPLANE_H", "PSI_NORM"), ep.time_slice[:].profiles_1d_r_midplane_h.psi_norm),
133 (("PROFILES_1D", "R_MIDPLANE_H", "PSI"), ep.time_slice[:].profiles_1d_r_midplane_h.psi),
134 # Profiles_2d. `profiles_2d(0)` because GSFit solves on a single rectangular (R, Z) grid
135 (("PROFILES_2D", "R_Z", "B_FIELD_PHI"), ep.time_slice[:].profiles_2d[0].b_field_phi),
136 (("PROFILES_2D", "R_Z", "B_FIELD_R"), ep.time_slice[:].profiles_2d[0].b_field_r),
137 (("PROFILES_2D", "R_Z", "B_FIELD_Z"), ep.time_slice[:].profiles_2d[0].b_field_z),
138 (("PROFILES_2D", "R_Z", "MASK"), ep.time_slice[:].profiles_2d[0].mask),
139 (("PROFILES_2D", "R_Z", "PRESSURE"), ep.time_slice[:].profiles_2d[0].pressure),
140 (("PROFILES_2D", "R_Z", "PSI"), ep.time_slice[:].profiles_2d[0].psi),
141 (("PROFILES_2D", "R_Z", "PHI"), ep.time_slice[:].profiles_2d[0].phi),
142 (("PROFILES_2D", "R_Z", "THETA"), ep.time_slice[:].profiles_2d[0].theta),
143 (("PROFILES_2D", "R_Z", "J_PHI"), ep.time_slice[:].profiles_2d[0].j_phi),
144 (("PROFILES_2D", "R_Z", "J_PARALLEL"), ep.time_slice[:].profiles_2d[0].j_parallel),
145 # Constraints
146 (("CONSTRAINTS", "CHI_SQ_MAG"), ep.time_slice[:].constraints.chi_squared_reduced),
147 # Scrape off layer (SOL)
148 (("SOL", "HFS", "CONTOUR", "R"), ep.time_slice[:].sol.hfs.contour.r),
149 (("SOL", "HFS", "CONTOUR", "Z"), ep.time_slice[:].sol.hfs.contour.z),
150 (("SOL", "HFS", "STRIKE_POINT", "R"), ep.time_slice[:].sol.hfs.strike_point.r),
151 (("SOL", "HFS", "STRIKE_POINT", "Z"), ep.time_slice[:].sol.hfs.strike_point.z),
152 (("SOL", "LFS", "CONTOUR", "R"), ep.time_slice[:].sol.lfs.contour.r),
153 (("SOL", "LFS", "CONTOUR", "Z"), ep.time_slice[:].sol.lfs.contour.z),
154 (("SOL", "LFS", "STRIKE_POINT", "R"), ep.time_slice[:].sol.lfs.strike_point.r),
155 (("SOL", "LFS", "STRIKE_POINT", "Z"), ep.time_slice[:].sol.lfs.strike_point.z),
156]
158# The same, for quantities which are the same on every time-slice. The IDS stores them per
159# time-slice all the same, so the gather returns `[n_time, n_points]` and the first row is taken.
160_TIME_INDEPENDENT_PATH_PAIRS: list[tuple[tuple[str, ...], typing.Any]] = [
161 (("PROFILES_1D", "PSI_NORM", "PSI_NORM"), ep.time_slice[:].profiles_1d.psi_norm),
162 (("PROFILES_1D", "R_MIDPLANE", "R"), ep.time_slice[:].profiles_1d_r_midplane.r),
163 (("PROFILES_1D", "R_MIDPLANE_H", "R"), ep.time_slice[:].profiles_1d_r_midplane_h.r),
164 (("PROFILES_2D", "R_Z", "R"), ep.time_slice[:].profiles_2d[0].grid.dim1),
165 (("PROFILES_2D", "R_Z", "Z"), ep.time_slice[:].profiles_2d[0].grid.dim2),
166]
169def _assign(results: typing.Any, mdsplus_path: tuple[str, ...], value: typing.Any) -> None:
170 """Assign `value` into the nested `results` object at `mdsplus_path`."""
171 node = results
172 for key in mdsplus_path[:-1]:
173 node = node[key]
174 node[mdsplus_path[-1]] = value
177def _n_points_per_time(padded: npt.NDArray[np.float64]) -> npt.NDArray[np.int32]:
178 """Number of real points in each row of a NaN-padded `[n_time, n_points]` array.
180 The IDS stores one contour per time-slice, each its own length; gathering over time pads the
181 short rows with NaN to make a rectangle. MDSplus needs that rectangle plus a separate count of
182 how much of each row is real, which is what this recovers. Contour coordinates are always
183 finite, so a NaN can only be padding.
184 """
185 n_points: npt.NDArray[np.int32] = np.isfinite(padded).sum(axis=1).astype(np.int32)
187 return n_points
190def map_results_to_database(
191 self: "DatabaseWriterTokamakEnergyMDSplusNew",
192 gsfit_controller: "Gsfit",
193) -> None:
194 """Map the results to MDSplus structure.
195 `gsfit_controller.results` is a `NestedDict` object, which is a 1:1 mapping to the MDSplus structure.
196 This function will mutate the `gsfit_controller` object.
197 """
199 # Take class object out of the `gsfit_controller` object
200 pulseNo = gsfit_controller.pulseNo
201 settings = gsfit_controller.settings
202 plasma = gsfit_controller.plasma
203 bp_probes = gsfit_controller.bp_probes
204 flux_loops = gsfit_controller.flux_loops
205 dialoop = gsfit_controller.dialoop
206 rogowski_coils = gsfit_controller.rogowski_coils
207 passives = gsfit_controller.passives
208 coils = gsfit_controller.coils
209 pressure_sensors = gsfit_controller.pressure_sensors
210 results = gsfit_controller.results
212 # Everything below comes out of the `equilibrium` IDS through the (MDSplus path, IMAS path)
213 # tables at the top of this file. `plasma.get` reads in place, so reading quantity by quantity
214 # costs nothing extra
216 for mdsplus_path, imas_path in _TIME_SERIES_PATH_PAIRS:
217 _assign(results, mdsplus_path, plasma.get(imas_path))
219 for mdsplus_path, imas_path in _TIME_INDEPENDENT_PATH_PAIRS:
220 _assign(results, mdsplus_path, np.asarray(plasma.get(imas_path))[0])
222 # The rest need something done to them on the way out, so they are not in the tables
224 # Flux defining the LCFS. SPIDER has `psi_norm = 0.9999`
225 results["BOUNDARY"]["PSI_NORM"] = np.ones_like(results["BOUNDARY"]["PSI"])
227 # Contour lengths. The IDS stores each contour at its own length and the gather pads them into
228 # a rectangle; MDSplus wants the rectangle plus the count of real points in each row
229 results["BOUNDARY"]["OUTLINE"]["N"] = _n_points_per_time(results["BOUNDARY"]["OUTLINE"]["R"])
230 results["SOL"]["HFS"]["CONTOUR"]["N"] = _n_points_per_time(results["SOL"]["HFS"]["CONTOUR"]["R"])
231 results["SOL"]["LFS"]["CONTOUR"]["N"] = _n_points_per_time(results["SOL"]["LFS"]["CONTOUR"]["R"])
233 # The data dictionary defines `beta_tor` as a fraction, but this MDSplus node has always held a
234 # percentage, so the factor of 100 is put back here rather than changing what consumers read
235 results["GLOBAL"]["BETA_T"] = 100.0 * np.asarray(plasma.get(ep.time_slice[:].global_quantities.beta_tor))
237 # The rod current is not a data dictionary node. It is recovered from the vacuum toroidal
238 # field, `f_vac = r0 * b0 = mu_0 * i_rod / (2 * pi)`
239 r0 = plasma.get(ep.vacuum_toroidal_field.r0)
240 b0 = np.asarray(plasma.get(ep.vacuum_toroidal_field.b0))
241 results["GLOBAL"]["I_ROD"] = 2.0 * np.pi * r0 * b0 / mu_0
244 for sensor_name in bp_probes.keys():
245 # results["CONSTRAINTS"]["BP_PROBE"][pf_name]["EXACT"]
246 results["CONSTRAINTS"]["BP_PROBE"][sensor_name]["INCLUDE"] = np.int32(bp_probes.get_bool([sensor_name, "fit_settings", "include"]))
247 results["CONSTRAINTS"]["BP_PROBE"][sensor_name]["MEASURED"] = bp_probes.get_array1([sensor_name, "b", "measured", "value"])
248 results["CONSTRAINTS"]["BP_PROBE"][sensor_name]["RECONSTRUCT"] = bp_probes.get_array1([sensor_name, "b", "calculated", "value"])
249 results["CONSTRAINTS"]["BP_PROBE"][sensor_name]["WEIGHT"] = bp_probes.get_f64([sensor_name, "fit_settings", "weight"])
251 for sensor_name in flux_loops.keys():
252 # results["CONSTRAINTS"]["FLUX_LOOP"][pf_name]["EXACT"]
253 results["CONSTRAINTS"]["FLUX_LOOP"][sensor_name]["INCLUDE"] = np.int32(flux_loops.get_bool([sensor_name, "fit_settings", "include"]))
254 results["CONSTRAINTS"]["FLUX_LOOP"][sensor_name]["MEASURED"] = flux_loops.get_array1([sensor_name, "psi", "measured", "value"])
255 results["CONSTRAINTS"]["FLUX_LOOP"][sensor_name]["RECONSTRUCT"] = flux_loops.get_array1([sensor_name, "psi", "calculated", "value"])
256 results["CONSTRAINTS"]["FLUX_LOOP"][sensor_name]["WEIGHT"] = flux_loops.get_f64([sensor_name, "fit_settings", "weight"])
258 for sensor_name in rogowski_coils.keys():
259 # results["CONSTRAINTS"]["ROGOWSKI"][pf_name]["EXACT"]
260 results["CONSTRAINTS"]["ROGOWSKI"][sensor_name]["INCLUDE"] = np.int32(rogowski_coils.get_bool([sensor_name, "fit_settings", "include"]))
261 results["CONSTRAINTS"]["ROGOWSKI"][sensor_name]["MEASURED"] = rogowski_coils.get_array1([sensor_name, "i", "measured", "value"])
262 results["CONSTRAINTS"]["ROGOWSKI"][sensor_name]["RECONSTRUCT"] = rogowski_coils.get_array1([sensor_name, "i", "calculated", "value"])
263 results["CONSTRAINTS"]["ROGOWSKI"][sensor_name]["WEIGHT"] = rogowski_coils.get_f64([sensor_name, "fit_settings", "weight"])
265 # Diamagnetic flux (single diamagnetic flux loop "DIALOOP")
266 for sensor_name in dialoop.keys():
267 results["CONSTRAINTS"]["DIAMAG_FLUX"]["INCLUDE"] = np.int32(dialoop.get_bool([sensor_name, "fit_settings", "include"]))
268 results["CONSTRAINTS"]["DIAMAG_FLUX"]["MEASURED"] = dialoop.get_array1([sensor_name, "b", "measured", "value"])
269 results["CONSTRAINTS"]["DIAMAG_FLUX"]["RECONSTRUCT"] = dialoop.get_array1([sensor_name, "b", "calculated", "value"])
270 results["CONSTRAINTS"]["DIAMAG_FLUX"]["WEIGHT"] = dialoop.get_f64([sensor_name, "fit_settings", "weight"])
272 for pf_name in coils.keys(["pf"]):
273 # results["CONSTRAINTS"]["PF_CURRENT"][pf_name]["EXACT"]
274 # results["CONSTRAINTS"]["PF_CURRENT"][pf_name]["INCLUDE"] = np.int32(coils.get_bool(["pf", pf_name, "fit_settings", "include"]))
275 results["CONSTRAINTS"]["PF_CURRENT"][pf_name]["MEASURED"] = coils.get_array1(["pf", pf_name, "i", "measured", "value"])
276 # results["CONSTRAINTS"]["PF_CURRENT"][pf_name]["RECONSTRUCT"] = coils.get_array1(["pf", pf_name, "i", "calculated", "value"])
277 # results["CONSTRAINTS"]["PF_CURRENT"][pf_name]["WEIGHT"] = coils.get_f64(["pf", pf_name, "fit_settings", "weight"])
278 # TODO: need to handle circuits vs coils better
279 results["CONSTRAINTS"]["PF_CURRENT"]["BVL"]["MEASURED"] = coils.get_array1(["pf", "BVLT", "i", "measured", "value"])
280 results["CONSTRAINTS"]["PF_CURRENT"]["DIV"]["MEASURED"] = coils.get_array1(["pf", "DIVT", "i", "measured", "value"])
281 results["CONSTRAINTS"]["PF_CURRENT"]["PSH"]["MEASURED"] = coils.get_array1(["pf", "PSHT", "i", "measured", "value"])
283 # # X-points
284 # # TODO!!!!!!!!!!!
285 # results["XPOINTS"]["UPPER"]["R"] = plasma.get_array1(["xpoints", "upper", "r"])
286 # results["XPOINTS"]["UPPER"]["Z"] = plasma.get_array1(["xpoints", "upper", "z"])
287 # results["XPOINTS"]["LOWER"]["R"] = plasma.get_array1(["xpoints", "lower", "r"])
288 # results["XPOINTS"]["LOWER"]["Z"] = plasma.get_array1(["xpoints", "lower", "z"])
290 # Passives
291 for passive_name in passives.keys():
292 if passive_name == "IVC":
293 results["CONSTRAINTS"]["PF_PASSIVE"]["IVC"]["DOF"]["EIG_01"]["RECONSTRUCT"] = passives.get_array1(["IVC", "dof", "eig_01", "calculated"])
294 results["CONSTRAINTS"]["PF_PASSIVE"]["IVC"]["DOF"]["EIG_01"]["I_DIST"] = passives.get_array1(["IVC", "dof", "eig_01", "current_distribution"])
295 results["CONSTRAINTS"]["PF_PASSIVE"]["IVC"]["DOF"]["EIG_02"]["RECONSTRUCT"] = passives.get_array1(["IVC", "dof", "eig_02", "calculated"])
296 results["CONSTRAINTS"]["PF_PASSIVE"]["IVC"]["DOF"]["EIG_02"]["I_DIST"] = passives.get_array1(["IVC", "dof", "eig_02", "current_distribution"])
297 results["CONSTRAINTS"]["PF_PASSIVE"]["IVC"]["DOF"]["EIG_03"]["RECONSTRUCT"] = passives.get_array1(["IVC", "dof", "eig_03", "calculated"])
298 results["CONSTRAINTS"]["PF_PASSIVE"]["IVC"]["DOF"]["EIG_03"]["I_DIST"] = passives.get_array1(["IVC", "dof", "eig_03", "current_distribution"])
299 results["CONSTRAINTS"]["PF_PASSIVE"]["IVC"]["DOF"]["EIG_04"]["RECONSTRUCT"] = passives.get_array1(["IVC", "dof", "eig_04", "calculated"])
300 results["CONSTRAINTS"]["PF_PASSIVE"]["IVC"]["DOF"]["EIG_04"]["I_DIST"] = passives.get_array1(["IVC", "dof", "eig_04", "current_distribution"])
301 results["CONSTRAINTS"]["PF_PASSIVE"]["IVC"]["DOF"]["EIG_05"]["RECONSTRUCT"] = passives.get_array1(["IVC", "dof", "eig_05", "calculated"])
302 results["CONSTRAINTS"]["PF_PASSIVE"]["IVC"]["DOF"]["EIG_05"]["I_DIST"] = passives.get_array1(["IVC", "dof", "eig_05", "current_distribution"])
303 results["CONSTRAINTS"]["PF_PASSIVE"]["IVC"]["DOF"]["EIG_06"]["RECONSTRUCT"] = passives.get_array1(["IVC", "dof", "eig_06", "calculated"])
304 results["CONSTRAINTS"]["PF_PASSIVE"]["IVC"]["DOF"]["EIG_06"]["I_DIST"] = passives.get_array1(["IVC", "dof", "eig_06", "current_distribution"])
305 results["CONSTRAINTS"]["PF_PASSIVE"]["IVC"]["DOF"]["EIG_07"]["RECONSTRUCT"] = passives.get_array1(["IVC", "dof", "eig_07", "calculated"])
306 results["CONSTRAINTS"]["PF_PASSIVE"]["IVC"]["DOF"]["EIG_07"]["I_DIST"] = passives.get_array1(["IVC", "dof", "eig_07", "current_distribution"])
307 results["CONSTRAINTS"]["PF_PASSIVE"]["IVC"]["DOF"]["EIG_08"]["RECONSTRUCT"] = passives.get_array1(["IVC", "dof", "eig_08", "calculated"])
308 results["CONSTRAINTS"]["PF_PASSIVE"]["IVC"]["DOF"]["EIG_08"]["I_DIST"] = passives.get_array1(["IVC", "dof", "eig_08", "current_distribution"])
309 results["CONSTRAINTS"]["PF_PASSIVE"]["IVC"]["DOF"]["EIG_09"]["RECONSTRUCT"] = passives.get_array1(["IVC", "dof", "eig_09", "calculated"])
310 results["CONSTRAINTS"]["PF_PASSIVE"]["IVC"]["DOF"]["EIG_09"]["I_DIST"] = passives.get_array1(["IVC", "dof", "eig_09", "current_distribution"])
311 results["CONSTRAINTS"]["PF_PASSIVE"]["IVC"]["DOF"]["EIG_10"]["RECONSTRUCT"] = passives.get_array1(["IVC", "dof", "eig_10", "calculated"])
312 results["CONSTRAINTS"]["PF_PASSIVE"]["IVC"]["DOF"]["EIG_10"]["I_DIST"] = passives.get_array1(["IVC", "dof", "eig_10", "current_distribution"])
313 results["CONSTRAINTS"]["PF_PASSIVE"]["IVC"]["DOF"]["EIG_11"]["RECONSTRUCT"] = passives.get_array1(["IVC", "dof", "eig_11", "calculated"])
314 results["CONSTRAINTS"]["PF_PASSIVE"]["IVC"]["DOF"]["EIG_11"]["I_DIST"] = passives.get_array1(["IVC", "dof", "eig_11", "current_distribution"])
315 results["CONSTRAINTS"]["PF_PASSIVE"]["IVC"]["DOF"]["EIG_12"]["RECONSTRUCT"] = passives.get_array1(["IVC", "dof", "eig_12", "calculated"])
316 results["CONSTRAINTS"]["PF_PASSIVE"]["IVC"]["DOF"]["EIG_12"]["I_DIST"] = passives.get_array1(["IVC", "dof", "eig_12", "current_distribution"])
317 results["CONSTRAINTS"]["PF_PASSIVE"]["IVC"]["DOF"]["EIG_13"]["RECONSTRUCT"] = passives.get_array1(["IVC", "dof", "eig_13", "calculated"])
318 results["CONSTRAINTS"]["PF_PASSIVE"]["IVC"]["DOF"]["EIG_13"]["I_DIST"] = passives.get_array1(["IVC", "dof", "eig_13", "current_distribution"])
319 results["CONSTRAINTS"]["PF_PASSIVE"]["IVC"]["DOF"]["EIG_14"]["RECONSTRUCT"] = passives.get_array1(["IVC", "dof", "eig_14", "calculated"])
320 results["CONSTRAINTS"]["PF_PASSIVE"]["IVC"]["DOF"]["EIG_14"]["I_DIST"] = passives.get_array1(["IVC", "dof", "eig_14", "current_distribution"])
321 results["CONSTRAINTS"]["PF_PASSIVE"]["IVC"]["DOF"]["EIG_15"]["RECONSTRUCT"] = passives.get_array1(["IVC", "dof", "eig_15", "calculated"])
322 results["CONSTRAINTS"]["PF_PASSIVE"]["IVC"]["DOF"]["EIG_15"]["I_DIST"] = passives.get_array1(["IVC", "dof", "eig_15", "current_distribution"])
323 results["CONSTRAINTS"]["PF_PASSIVE"]["IVC"]["GEOMETRY"]["ANGLE_1"] = passives.get_array1(["IVC", "geometry", "angle_1"])
324 results["CONSTRAINTS"]["PF_PASSIVE"]["IVC"]["GEOMETRY"]["ANGLE_2"] = passives.get_array1(["IVC", "geometry", "angle_2"])
325 results["CONSTRAINTS"]["PF_PASSIVE"]["IVC"]["GEOMETRY"]["D_R"] = passives.get_array1(["IVC", "geometry", "d_r"])
326 results["CONSTRAINTS"]["PF_PASSIVE"]["IVC"]["GEOMETRY"]["D_Z"] = passives.get_array1(["IVC", "geometry", "d_z"])
327 results["CONSTRAINTS"]["PF_PASSIVE"]["IVC"]["GEOMETRY"]["R"] = passives.get_array1(["IVC", "geometry", "r"])
328 results["CONSTRAINTS"]["PF_PASSIVE"]["IVC"]["GEOMETRY"]["Z"] = passives.get_array1(["IVC", "geometry", "z"])
329 else:
330 results["CONSTRAINTS"]["PF_PASSIVE"][passive_name]["DOF"]["CONSTANT_J"]["RECONSTRUCT"] = passives.get_array1(
331 [passive_name, "dof", "constant_current_density", "calculated"]
332 )
333 results["CONSTRAINTS"]["PF_PASSIVE"][passive_name]["DOF"]["CONSTANT_J"]["I_DIST"] = passives.get_array1(
334 [passive_name, "dof", "constant_current_density", "current_distribution"]
335 )
336 results["CONSTRAINTS"]["PF_PASSIVE"][passive_name]["GEOMETRY"]["ANGLE_1"] = passives.get_array1([passive_name, "geometry", "angle_1"])
337 results["CONSTRAINTS"]["PF_PASSIVE"][passive_name]["GEOMETRY"]["ANGLE_2"] = passives.get_array1([passive_name, "geometry", "angle_2"])
338 results["CONSTRAINTS"]["PF_PASSIVE"][passive_name]["GEOMETRY"]["D_R"] = passives.get_array1([passive_name, "geometry", "d_r"])
339 results["CONSTRAINTS"]["PF_PASSIVE"][passive_name]["GEOMETRY"]["D_Z"] = passives.get_array1([passive_name, "geometry", "d_z"])
340 results["CONSTRAINTS"]["PF_PASSIVE"][passive_name]["GEOMETRY"]["R"] = passives.get_array1([passive_name, "geometry", "r"])
341 results["CONSTRAINTS"]["PF_PASSIVE"][passive_name]["GEOMETRY"]["Z"] = passives.get_array1([passive_name, "geometry", "z"])
344 if len(pressure_sensors.keys()) > 0:
345 sensor_names = list(pressure_sensors.keys())
347 # Per-sensor nodes (keyed by sensor name)
348 for sensor_name in sensor_names:
349 results["CONSTRAINTS"]["PRESSURE"][sensor_name]["MEASURED"] = pressure_sensors.get_array1(
350 [sensor_name, "pressure", "measured", "value"]
351 ) # shape = [n_time]
352 results["CONSTRAINTS"]["PRESSURE"][sensor_name]["RECONSTRUCT"] = pressure_sensors.get_array1(
353 [sensor_name, "pressure", "calculated", "value"]
354 ) # shape = [n_time]
355 results["CONSTRAINTS"]["PRESSURE"][sensor_name]["WEIGHT"] = pressure_sensors.get_f64([sensor_name, "fit_settings", "weight"]) # scalar
356 results["CONSTRAINTS"]["PRESSURE"][sensor_name]["POSITION"]["R"] = pressure_sensors.get_f64([sensor_name, "geometry", "r"]) # scalar
357 results["CONSTRAINTS"]["PRESSURE"][sensor_name]["POSITION"]["Z"] = pressure_sensors.get_f64([sensor_name, "geometry", "z"]) # scalar
358 results["CONSTRAINTS"]["PRESSURE"][sensor_name]["POSITION"]["PSI"] = pressure_sensors.get_array1(
359 [sensor_name, "pressure", "calculated", "psi"]
360 ) # shape = [n_time]
362 # ALL aggregate node
363 results["CONSTRAINTS"]["PRESSURE"]["ALL"]["MEASURED"] = pressure_sensors.get_array2(
364 ["*", "pressure", "measured", "value"]
365 ) # shape = [n_time, n_points]
366 results["CONSTRAINTS"]["PRESSURE"]["ALL"]["RECONSTRUCT"] = pressure_sensors.get_array2(
367 ["*", "pressure", "calculated", "value"]
368 ) # shape = [n_time, n_points]
369 results["CONSTRAINTS"]["PRESSURE"]["ALL"]["WEIGHT"] = pressure_sensors.get_array1(["*", "fit_settings", "weight"]) # shape = [n_points]
370 results["CONSTRAINTS"]["PRESSURE"]["ALL"]["POSITION"]["R"] = pressure_sensors.get_array1(["*", "geometry", "r"]) # shape = [n_points]
371 results["CONSTRAINTS"]["PRESSURE"]["ALL"]["POSITION"]["Z"] = pressure_sensors.get_array1(["*", "geometry", "z"]) # shape = [n_points]
372 results["CONSTRAINTS"]["PRESSURE"]["ALL"]["POSITION"]["PSI"] = pressure_sensors.get_array2(
373 ["*", "pressure", "calculated", "psi"]
374 ) # shape = [n_time, n_points]
375 results["CONSTRAINTS"]["PRESSURE"]["ALL"]["NAMES"] = np.array(sensor_names)
377 # Store "WORKFLOW"
378 database_reader_method = settings["GSFIT_code_settings.json"]["database_reader"]["method"]
380 code_names = settings["GSFIT_code_settings.json"]["database_reader"][database_reader_method]["workflow"].keys()
382 for code_name in code_names:
383 pulseNo_json = settings["GSFIT_code_settings.json"]["database_reader"][database_reader_method]["workflow"][code_name]["pulseNo"]
384 if pulseNo_json is not None:
385 results["INPUT"]["WORKFLOW"][code_name]["PULSE"] = pulseNo_json
386 else:
387 results["INPUT"]["WORKFLOW"][code_name]["PULSE"] = pulseNo
389 run_name = settings["GSFIT_code_settings.json"]["database_reader"][database_reader_method]["workflow"][code_name]["run_name"]
390 results["INPUT"]["WORKFLOW"][code_name]["RUN"] = run_name
392 usage = settings["GSFIT_code_settings.json"]["database_reader"][database_reader_method]["workflow"][code_name]["usage"]
393 results["INPUT"]["WORKFLOW"][code_name]["USAGE"] = usage