Coverage for src/hallmd/data/__init__.py: 94%
142 statements
« prev ^ index » next coverage.py v7.6.10, created at 2025-02-20 20:43 +0000
1"""The `hallmd.data` package contains a folder for each unique thruster. The experimental data for each thruster
2is further divided by folders for each individual paper or reference. The raw experimental data is contained within
3these folders in any supported format (currently only .csv). Any additional documentation
4for the datasets is encouraged (e.g. citations, descriptions, summaries, etc.) and can be included in the data folders.
6## Thrusters
8### SPT-100
9Currently the only thruster with available data. Data for the SPT-100 comes from four sources:
111. [Diamant et al. 2014](https://arc.aiaa.org/doi/10.2514/6.2014-3710) - provides thrust and ion current density data as a function of chamber background pressure.
122. [Macdonald et al. 2019](https://arc.aiaa.org/doi/10.2514/1.B37133) - provides ion velocity profiles for varying chamber pressures.
133. [Sankovic et al. 1993](https://www.semanticscholar.org/paper/Performance-evaluation-of-the-Russian-SPT-100-at-Sankovic-Hamley/81b7d985669b21aa1a8419277c52e7a879bf3b46) - provides thrust at varying operating conditions.
144. [Jorns and Byrne. 2021](https://pepl.engin.umich.edu/pdf/2021_PSST_Jorns.pdf) - provides cathode coupling voltages at same conditions as Diamant et al. 2014.
16Citations:
17``` title="SPT-100.bib"
18--8<-- "hallmd/data/spt100/spt100.bib:citations"
19```
20""" # noqa: E501
22from dataclasses import dataclass, fields
23from pathlib import Path
24from typing import Any, Generic, Optional, TypeAlias, TypeVar
26import numpy as np
27import numpy.typing as npt
29Array: TypeAlias = npt.NDArray[np.floating[Any]]
30PathLike: TypeAlias = str | Path
33@dataclass(frozen=True)
34class OperatingCondition:
35 """Operating conditions for a Hall thruster. Currently includes background pressure (Torr),
36 discharge voltage (V), and anode mass flow rate (kg/s).
37 """
38 background_pressure_Torr: float
39 discharge_voltage_V: float
40 anode_mass_flow_rate_kg_s: float
43T = TypeVar("T", np.float64, Array)
46@dataclass(frozen=True)
47class Measurement(Generic[T]):
48 """A measurement object that includes a mean and standard deviation. The mean is the best estimate of the
49 quantity being measured, and the standard deviation is the uncertainty in the measurement. Can be used to specify
50 a scalar measurement quantity or a field quantity (e.g. a profile) in the form of a `numpy` array.
51 """
52 mean: T
53 std: T
55 def __str__(self):
56 return f"(μ = {self.mean}, σ = {self.std})"
59def _gauss_logpdf(mean: T, std: T, observation: T) -> np.float64:
60 var = std**2
61 term1 = np.mean(np.log(2 * np.pi * var)) / 2
62 term2 = np.mean((mean - observation) ** 2 / (2 * var))
63 return -term1 - term2
66def _measurement_gauss_logpdf(data: Measurement[T] | None, observation: Measurement[T] | None) -> np.float64:
67 if data is None or observation is None:
68 return np.float64(0.0)
70 return _gauss_logpdf(data.mean, data.std, observation.mean)
73def _interp_gauss_logpdf(
74 coords: Array | None,
75 data: Measurement[Array] | None,
76 obs_coords: Array | None,
77 observation: Measurement[Array] | None,
78) -> np.float64:
79 if coords is None or data is None or obs_coords is None or observation is None:
80 return np.float64(0.0)
82 obs_interp_mean = np.interp(coords, obs_coords, observation.mean)
83 obs_interp = Measurement(obs_interp_mean, np.zeros(1))
84 return _measurement_gauss_logpdf(data, obs_interp)
87@dataclass
88class ThrusterData:
89 """Class for Hall thruster data. Contains fields for all relevant performance metrics and quantities of interest."""
90 # Cathode
91 cathode_coupling_voltage_V: Optional[Measurement[np.float64]] = None
92 # Thruster
93 thrust_N: Optional[Measurement[np.float64]] = None
94 discharge_current_A: Optional[Measurement[np.float64]] = None
95 ion_current_A: Optional[Measurement[np.float64]] = None
96 efficiency_current: Optional[Measurement[np.float64]] = None
97 efficiency_mass: Optional[Measurement[np.float64]] = None
98 efficiency_voltage: Optional[Measurement[np.float64]] = None
99 efficiency_anode: Optional[Measurement[np.float64]] = None
100 ion_velocity_coords_m: Optional[Array] = None
101 ion_velocity_m_s: Optional[Measurement[Array]] = None
102 # Plume
103 divergence_angle_rad: Optional[Measurement[np.float64]] = None
104 ion_current_density_radius_m: Optional[float] = None
105 ion_current_density_coords_m: Optional[Array] = None
106 ion_current_density_A_m2: Optional[Measurement[Array]] = None
108 def __str__(self) -> str:
109 fields_str = ",\n".join([f"\t{field.name} = {val}" for field in fields(ThrusterData)
110 if (val := getattr(self, field.name)) is not None])
111 return f"ThrusterData(\n{fields_str}\n)\n"
114def pem_to_thrusterdata(
115 operating_conditions: list[OperatingCondition], outputs
116) -> dict[OperatingCondition, ThrusterData]:
117 # Assemble output dict from operating conditions -> results
118 # Note, we assume that the pem outputs are ordered based on the input operating conditions
119 NaN = np.float64(np.nan)
120 output_dict = {
121 opcond: ThrusterData(
122 cathode_coupling_voltage_V=Measurement(outputs["V_cc"][i], NaN),
123 thrust_N=Measurement(outputs["T"][i], NaN),
124 discharge_current_A=Measurement(outputs["I_d"][i], NaN),
125 ion_current_A=Measurement(outputs["I_B0"][i], NaN),
126 ion_velocity_coords_m=outputs["u_ion_coords"][i],
127 ion_velocity_m_s=Measurement(outputs["u_ion"][i], np.full_like(outputs["u_ion"][i], NaN)),
128 ion_current_density_coords_m=outputs["j_ion_coords"][i],
129 ion_current_density_A_m2=Measurement(outputs["j_ion"][i], np.full_like(outputs["j_ion"][i], NaN)),
130 ion_current_density_radius_m=1.0,
131 efficiency_mass=Measurement(outputs["eta_m"][i], NaN),
132 efficiency_current=Measurement(outputs["eta_c"][i], NaN),
133 efficiency_voltage=Measurement(outputs["eta_v"][i], NaN),
134 efficiency_anode=Measurement(outputs["eta_a"][i], NaN),
135 )
136 for (i, opcond) in enumerate(operating_conditions)
137 }
139 return output_dict
142def log_likelihood(data: ThrusterData, observation: ThrusterData) -> np.float64:
143 log_likelihood = (
144 # Add contributions from global performance metrics
145 _measurement_gauss_logpdf(data.cathode_coupling_voltage_V, observation.cathode_coupling_voltage_V)
146 + _measurement_gauss_logpdf(data.thrust_N, observation.thrust_N)
147 + _measurement_gauss_logpdf(data.discharge_current_A, observation.discharge_current_A)
148 + _measurement_gauss_logpdf(data.ion_current_A, observation.ion_current_A)
149 + _measurement_gauss_logpdf(data.efficiency_current, observation.efficiency_current)
150 + _measurement_gauss_logpdf(data.efficiency_mass, observation.efficiency_mass)
151 + _measurement_gauss_logpdf(data.efficiency_voltage, observation.efficiency_voltage)
152 + _measurement_gauss_logpdf(data.efficiency_anode, observation.efficiency_anode)
153 + _measurement_gauss_logpdf(data.divergence_angle_rad, observation.divergence_angle_rad)
154 # interpolated average pointwise error from ion velocity and ion current density
155 + _interp_gauss_logpdf(
156 data.ion_velocity_coords_m,
157 data.ion_velocity_m_s,
158 observation.ion_velocity_coords_m,
159 observation.ion_velocity_m_s,
160 )
161 + _interp_gauss_logpdf(
162 data.ion_current_density_coords_m,
163 data.ion_current_density_A_m2,
164 observation.ion_current_density_coords_m,
165 observation.ion_current_density_A_m2,
166 )
167 )
168 return log_likelihood
171def load(files: list[PathLike] | PathLike) -> dict[OperatingCondition, ThrusterData]:
172 """Load all data from the given files into a dict map of `OperatingCondition` -> `ThrusterData`.
173 Each thruster operating condition corresponds to one set of thruster measurements or quantities of interest (QoIs).
175 :param files: A list of file paths or a single file path to load data from (only .csv supported).
176 :return: A dict map of `OperatingCondition` -> `ThrusterData` objects.
177 """
178 data: dict[OperatingCondition, ThrusterData] = {}
179 if isinstance(files, list):
180 # Recursively load resources in this list (possibly list of lists)
181 for file in files:
182 data.update(load(file))
183 else:
184 data.update(_load_single(files))
186 return data
189def _load_single(file: PathLike) -> dict[OperatingCondition, ThrusterData]:
190 """Load data from a single file into a dict map of `OperatingCondition` -> `ThrusterData`."""
191 if not Path(file).suffix == '.csv':
192 raise ValueError(f"Unsupported file format: {Path(file).suffix}. Only .csv files are supported.")
194 table = _table_from_file(file, delimiter=",", comments="#")
195 data: dict[OperatingCondition, ThrusterData] = {}
197 # Compute anode mass flow rate, if not present
198 mdot_a_key = "anode flow rate (mg/s)"
199 mdot_t_key = "total flow rate (mg/s)"
200 flow_ratio_key = "anode-cathode flow ratio"
201 keys = list(table.keys())
203 if mdot_a_key in keys:
204 mdot_a = table[mdot_a_key] * 1e-6 # convert to kg/s
205 else:
206 if mdot_t_key in keys and flow_ratio_key in keys:
207 flow_ratio = table[flow_ratio_key]
208 anode_flow_fraction = flow_ratio / (flow_ratio + 1)
209 mdot_a = table[mdot_t_key] * anode_flow_fraction * 1e-6 # convert to kg/s
210 else:
211 raise KeyError(
212 f"{file}: No mass flow rate provided."
213 + " Need either `anode flow rate (mg/s)` or [`total flow rate (mg/s)` and `anode-cathode flow ratio`]"
214 )
216 # Get background pressure and discharge voltage
217 P_B = table["background pressure (torr)"]
218 V_a = table["anode voltage (v)"]
220 num_rows = len(table[keys[0]])
221 row_num = 0
222 opcond_start_row = 0
223 opcond = OperatingCondition(P_B[0], V_a[0], mdot_a[0])
225 while True:
226 next_opcond = opcond
227 if row_num < num_rows:
228 next_opcond = OperatingCondition(P_B[row_num], V_a[row_num], mdot_a[row_num])
229 if next_opcond == opcond:
230 row_num += 1
231 continue
233 # either at end of operating condition or end of table
234 # fill up thruster data object for this row
235 data[opcond] = ThrusterData()
237 # Fill up data
238 # We assume all errors (expressed as value +/- error) correspond to two standard deviations
239 for key, val in table.items():
240 if key == "thrust (mn)":
241 # Load thrust data
242 T = val[opcond_start_row] * 1e-3 # convert to Newtons
243 T_std = table["thrust relative uncertainty"][opcond_start_row] * T / 2
244 data[opcond].thrust_N = Measurement(mean=T, std=T_std)
246 elif key == "anode current (a)":
247 # Load discharge current data
248 # assume a 0.1-A std deviation for discharge current
249 data[opcond].discharge_current_A = Measurement(mean=val[opcond_start_row], std=np.float64(0.1))
251 elif key == "cathode coupling voltage (v)":
252 # Load cathode coupling data
253 V_cc = val[opcond_start_row]
254 V_cc_std = table["cathode coupling voltage absolute uncertainty (v)"][opcond_start_row] / 2
255 data[opcond].cathode_coupling_voltage_V = Measurement(mean=V_cc, std=V_cc_std)
257 elif key == "ion velocity (m/s)":
258 # Load ion velocity data
259 uion: Array = val[opcond_start_row:row_num]
260 uion_std: Array = table["ion velocity absolute uncertainty (m/s)"][opcond_start_row:row_num] / 2
261 data[opcond].ion_velocity_m_s = Measurement(mean=uion, std=uion_std)
262 data[opcond].ion_velocity_coords_m = table["axial position from anode (m)"][opcond_start_row:row_num]
264 elif key == "ion current density (ma/cm^2)":
265 # Load ion current density data
266 jion: Array = val[opcond_start_row:row_num] * 10 # Convert to A / m^2
267 jion_std: Array = table["ion current density relative uncertainty"][opcond_start_row:row_num] * jion / 2
268 r = table["radial position from thruster exit (m)"][0]
269 jion_coords: Array = table["angular position from thruster centerline (deg)"][opcond_start_row:row_num]
271 # Keep only measurements at angles less than 90 degrees
272 keep_inds = jion_coords < 90
273 data[opcond].ion_current_density_coords_m = jion_coords[keep_inds] * np.pi / 180
274 data[opcond].ion_current_density_radius_m = r
275 data[opcond].ion_current_density_A_m2 = Measurement(mean=jion[keep_inds], std=jion_std[keep_inds])
277 # Advance to next operating condition or break out of loop if we're at the end of the table
278 if row_num == num_rows:
279 break
281 opcond = next_opcond
282 opcond_start_row = row_num
284 return data
287def _table_from_file(file: PathLike, delimiter=",", comments="#") -> dict[str, Array]:
288 """Return a `dict` of `numpy` arrays from a CSV file. The keys of the dict are the column names in the CSV."""
289 # Read header of file to get column names
290 # We skip comments (lines starting with the string in the `comments` arg)
291 header_start = 0
292 header = ""
293 with open(file, "r") as f:
294 for i, line in enumerate(f):
295 if not line.startswith(comments):
296 header = line.rstrip()
297 header_start = i
298 break
300 if header == "":
301 return {}
303 column_names = header.split(delimiter)
304 data = np.genfromtxt(file, delimiter=delimiter, comments=comments, skip_header=header_start + 1)
306 table: dict[str, Array] = {column.casefold(): data[:, i] for (i, column) in enumerate(column_names)}
307 return table