From d4056a0853a37528e82c998d834e3ef74fac6269 Mon Sep 17 00:00:00 2001 From: andytorrestb Date: Mon, 27 Jul 2026 22:53:38 -0500 Subject: [PATCH 01/14] rpod: make n_firings mean exactly N JFH entries print_jfh_1d_approach_n_fire accepted an n_firings argument and then ignored it: the emitted Jet Firing History held whatever number of entries the deceleration simulation happened to produce. - add validate_n_firings() to approach_maneuvers as the single definition of the rule (positive integer; zero, negative, fractional and non-numeric values are rejected rather than coerced); - give compute_1d_approach an optional n_firings that stops the simulation at exactly that many entries and raises when the approach completes in fewer firings than requested, instead of silently returning a shorter history; - validate the argument in print_jfh_1d_approach_n_fire, pass it through, and assert the written file's entry count. Omitting n_firings keeps the previous behavior exactly, so no existing caller changes. Co-Authored-By: Claude Opus 5 --- pyrpod/rpod/PlumeStrikeEstimationStudy.py | 15 ++++++ pyrpod/rpod/approach_maneuvers.py | 56 +++++++++++++++++++++++ 2 files changed, 71 insertions(+) diff --git a/pyrpod/rpod/PlumeStrikeEstimationStudy.py b/pyrpod/rpod/PlumeStrikeEstimationStudy.py index edb7998..2659897 100644 --- a/pyrpod/rpod/PlumeStrikeEstimationStudy.py +++ b/pyrpod/rpod/PlumeStrikeEstimationStudy.py @@ -37,6 +37,7 @@ from pyrpod.rpod.approach_maneuvers import ( ApproachInputs, compute_1d_approach, + validate_n_firings, ) from pyrpod.rpod.io import ensure_results_dirs, write_jfh from pyrpod.rpod.PlumeStudyExport import PlumeStudyExport @@ -1411,6 +1412,15 @@ def calc_time_multiplier(self, v_ida: float, v_o: float, def print_jfh_1d_approach_n_fire(self, v_ida: float, v_o: float, r_o: float, n_firings: int, trade_study: bool = False) -> None: + """Write a 1D-approach JFH holding EXACTLY n_firings entries. + + ``n_firings`` is the number of entries written to the Jet Firing + History. It is validated up front and enforced by + ``compute_1d_approach``; previously it was accepted and then ignored, + so the file's length was whatever the deceleration simulation + happened to produce. + """ + n_firings = validate_n_firings(n_firings) # Delegate to approach_maneuvers.compute_1d_approach and rpod.io.write_jfh tm = self.calc_time_multiplier(v_ida, v_o, r_o) inputs = ApproachInputs(v_ida=float(v_ida), v_o=float(v_o), r_o=float(r_o), group='neg_x') @@ -1430,6 +1440,7 @@ def calc_v_e(self, group: str) -> float: grouping=_GroupingAdapter(self), cant_rad=self.vv.decel_cant, dt_strategy={"multiplier": tm}, + n_firings=n_firings, ) r = [results["x"], results["y"], results["z"]] @@ -1448,6 +1459,10 @@ def calc_v_e(self, group: str) -> float: gen_start = time.perf_counter() os.makedirs(os.path.dirname(jfh_path), exist_ok=True) write_jfh(t_values, r, rot, jfh_path, mode="1d") + if len(t_values) != n_firings: + raise ValueError( + f"JFH {jfh_path} was written with {len(t_values)} entries but " + f"{n_firings} firings were requested") _log_jfh_generation_complete(jfh_path, len(t_values), gen_start) diff --git a/pyrpod/rpod/approach_maneuvers.py b/pyrpod/rpod/approach_maneuvers.py index 42acaf9..b567b7c 100644 --- a/pyrpod/rpod/approach_maneuvers.py +++ b/pyrpod/rpod/approach_maneuvers.py @@ -21,6 +21,38 @@ class ApproachInputs: group: str = "neg_x" # thruster group for decel +def validate_n_firings(value: Any) -> int: + """Validate a requested firing count. + + ``n_firings`` is the EXACT number of entries a Jet Firing History will + contain, so it must be a positive integer: zero, negative, fractional and + non-numeric values are rejected rather than coerced or silently ignored. + + Parameters + ---------- + value : Any + The requested firing count. + + Returns + ------- + int + The validated count. + + Raises + ------ + ValueError + If the value is not a positive integer. + """ + if isinstance(value, bool) or not isinstance(value, (int, float)): + raise ValueError(f"n_firings must be a positive integer, got {value!r}") + if isinstance(value, float) and not float(value).is_integer(): + raise ValueError(f"n_firings must be a positive integer, got {value!r}") + n_firings = int(value) + if n_firings < 1: + raise ValueError(f"n_firings must be a positive integer, got {value!r}") + return n_firings + + def choose_dt_from_mib(vv: Any, group: str, time_multiplier: float = 1.0) -> float: """Default time step based on min impulse bit and a multiplier. @@ -46,15 +78,29 @@ def compute_1d_approach( grouping: Any, cant_rad: float, dt_strategy: Dict[str, Any] | None = None, + n_firings: int | None = None, ) -> Dict[str, np.ndarray]: """Discrete 1D deceleration under constant-thrust firings. Returns dict with arrays: t, x, y, z, v, dv, mass, dm_total, rot No file I/O; callers can serialize via rpod.io.write_jfh. + + Parameters + ---------- + n_firings : int, optional + When given, the returned arrays hold EXACTLY this many entries: the + deceleration is simulated until the requested number of firings has + been produced. If the approach reaches its target velocity in fewer + firings than requested, the mismatch is raised rather than silently + ignored -- ``n_firings`` is a hard count, not a hint. When omitted + (the default) the simulation runs to completion exactly as before. """ v_ida, v_o, r_o = inputs.v_ida, inputs.v_o, inputs.r_o group = inputs.group + if n_firings is not None: + n_firings = validate_n_firings(n_firings) + # Pre-compute thrust/mass-flow characteristics m_dot_sum = grouping.calc_m_dot_sum(group) v_e = grouping.calc_v_e(group) @@ -97,6 +143,10 @@ def compute_1d_approach( rot.append(rot_mat) while dv_req > 0: + # A requested firing count is exact: stop as soon as the arrays hold + # n_firings entries (one JFH entry per entry). + if n_firings is not None and len(t) >= n_firings: + break m_o = mass[-1] mass.append(m_o - dm_firing) @@ -119,6 +169,12 @@ def compute_1d_approach( # guard against runaway loops break + if n_firings is not None and len(t) != n_firings: + raise ValueError( + f"requested n_firings={n_firings} but the approach reached its " + f"target velocity after {len(t)} firings; reduce n_firings or " + "adjust the approach conditions (v_o, v_ida, time step)") + return { "t": np.array(t), "x": np.array(x), From f58a0815b7fe7e0edcb6965e45f8c84015393340 Mon Sep 17 00:00:00 2001 From: andytorrestb Date: Mon, 27 Jul 2026 22:53:38 -0500 Subject: [PATCH 02/14] mdao: add the YAML study configuration A study configuration is a declarative layer ON TOP OF an existing PyRPOD case: the case's config.ini keeps owning the vehicle, thruster, plume-model and target assets, so every existing case and public API is untouched, while the YAML adds only what a trade study needs and an INI cannot express -- swept plate angles and source distances, the sweep decomposition mode, the exact firing count, explicitly prescribed firings, the moment reference point, coefficient normalization, output/plot settings and an optional external reference-data location. Validation is strict and specific: a missing case directory, a case without a config.ini, an unsupported plume model, an unknown sweep mode, an empty or non-positive sweep axis, non-orthogonal target axes, duplicate component names, a non-positive normalization value or a firing list whose length disagrees with the mode's implied total each raise StudyConfigError naming the offending key. The plume model is recorded explicitly and must be SimplifiedGasKinetics; no model registry is introduced. n_firings has ONE meaning in both sweep modes -- entries contributed by each pose -- and reuses the rule defined with the JFH-generation code so the prescribed and dynamics-driven paths cannot drift apart. Co-Authored-By: Claude Opus 5 --- pyrpod/mdao/study_config.py | 647 ++++++++++++++++++++++++++++++++++++ 1 file changed, 647 insertions(+) create mode 100644 pyrpod/mdao/study_config.py diff --git a/pyrpod/mdao/study_config.py b/pyrpod/mdao/study_config.py new file mode 100644 index 0000000..465dda1 --- /dev/null +++ b/pyrpod/mdao/study_config.py @@ -0,0 +1,647 @@ +""" +YAML study configuration for prescribed plume/target validation sweeps. + +A study configuration is a thin, declarative layer ON TOP OF an existing +PyRPOD case directory: the case's ``config.ini`` keeps owning the vehicle, +thruster, plume-model and target assets (so every existing case and public +API keeps working unchanged), while the YAML file adds only what a trade +study needs and a ``config.ini`` cannot express -- swept plate angles and +source distances, prescribed firing counts, the moment reference point, +coefficient normalization, and output/plot/reference settings. + +The schema is deliberately explicit. Nothing is inferred that the caller +did not write down: + +* the plume model is recorded by name and must be ``SimplifiedGasKinetics`` + (this branch adds no plume-model registry); +* coefficients are computed only when every normalization input is present + (see :class:`Normalization`); otherwise they are reported as unavailable; +* ``n_firings`` means the exact number of entries written to the Jet Firing + History -- a mismatch against an explicitly prescribed firing list is an + error, never a silent truncation. + +Example +------- +>>> cfg = StudyConfig.from_yaml('case/.../flat_plate_baseline.yaml') +>>> cfg.study_name +'cai2016_flat_plate_baseline' +""" + +from __future__ import annotations + +import os +from dataclasses import dataclass, field, replace +from typing import Any, Mapping, Sequence + +import numpy as np +import yaml +from numpy.typing import NDArray + +from pyrpod.rpod.approach_maneuvers import ( + validate_n_firings as _validate_n_firings, +) + +#: The single plume model this study workflow supports (hard scope constraint). +SUPPORTED_PLUME_MODEL = "SimplifiedGasKinetics" + +#: Default unit labels carried into the result metadata. PyRPOD works in SI +#: throughout; recording them makes an exported result self-describing. +DEFAULT_UNITS: dict[str, str] = { + "length": "m", + "area": "m^2", + "force": "N", + "moment": "N*m", + "pressure": "Pa", + "shear_stress": "Pa", + "heat_flux": "W/m^2", + "time": "s", + "angle": "deg", +} + + +class StudyConfigError(ValueError): + """Raised when a study configuration is missing or internally inconsistent.""" + + +def _require(mapping: Mapping[str, Any], key: str, context: str) -> Any: + if key not in mapping: + raise StudyConfigError( + f"missing required key {key!r} in {context}") + return mapping[key] + + +def _as_vector(value: Any, context: str) -> NDArray[np.float64]: + """Parse a 3-element vector, failing loudly on anything else.""" + try: + vector = np.asarray(value, dtype=float).reshape(-1) + except (TypeError, ValueError) as exc: + raise StudyConfigError( + f"{context}: expected three numbers, got {value!r}") from exc + if vector.size != 3: + raise StudyConfigError( + f"{context}: expected three numbers, got {value!r}") + return vector + + +def _as_float_list(value: Any, context: str) -> list[float]: + if value is None or isinstance(value, (str, bytes)) or not isinstance( + value, Sequence): + raise StudyConfigError(f"{context}: expected a list of numbers, " + f"got {value!r}") + try: + return [float(v) for v in value] + except (TypeError, ValueError) as exc: + raise StudyConfigError(f"{context}: expected a list of numbers, " + f"got {value!r}") from exc + + +@dataclass(frozen=True) +class ComponentSpec: + """One target component: a named subset of the target mesh's faces. + + A component is selected either by ``face_indices`` (explicit) or by + ``bounds`` (an axis-aligned box on face centroids, in the case's global + frame). ``selector: all`` -- the default -- takes the whole mesh, which + is the right answer for the single-plate and single-cylinder targets this + branch ships. Nothing here assumes a flat plate. + + Attributes + ---------- + name : str + Component identifier, reported with every result row. + selector : str + ``'all'``, ``'face_indices'`` or ``'bounds'``. + face_indices : tuple of int + Explicit face indices when ``selector == 'face_indices'``. + bounds : dict + ``{'min': [x, y, z], 'max': [x, y, z]}`` when ``selector == 'bounds'``. + """ + + name: str + selector: str = "all" + face_indices: tuple[int, ...] = () + bounds: dict[str, tuple[float, float, float]] = field(default_factory=dict) + + @classmethod + def from_mapping(cls, data: Mapping[str, Any]) -> "ComponentSpec": + name = str(_require(data, "name", "target.components entry")) + if "face_indices" in data and data["face_indices"] is not None: + indices = tuple(int(i) for i in data["face_indices"]) + if not indices: + raise StudyConfigError( + f"target component {name!r}: face_indices is empty") + return cls(name=name, selector="face_indices", face_indices=indices) + if "bounds" in data and data["bounds"] is not None: + raw = data["bounds"] + lo = _as_vector(_require(raw, "min", f"component {name!r} bounds"), + f"component {name!r} bounds.min") + hi = _as_vector(_require(raw, "max", f"component {name!r} bounds"), + f"component {name!r} bounds.max") + if np.any(hi < lo): + raise StudyConfigError( + f"target component {name!r}: bounds.max must be >= " + "bounds.min componentwise") + bounds = {"min": (float(lo[0]), float(lo[1]), float(lo[2])), + "max": (float(hi[0]), float(hi[1]), float(hi[2]))} + return cls(name=name, selector="bounds", bounds=bounds) + selector = str(data.get("selector", "all")) + if selector != "all": + raise StudyConfigError( + f"target component {name!r}: selector {selector!r} needs " + "'face_indices' or 'bounds' to be supplied") + return cls(name=name, selector="all") + + +@dataclass(frozen=True) +class TargetSpec: + """Target geometry description and its component breakdown. + + ``reference_point`` is the geometric reference the sweep is built about + (the plate center, the cylinder axis midpoint, ...). ``normal`` and + ``tangent`` define the plane the approach angle is swept in: angle 0 + places the plume source on ``normal``, positive angles rotate it toward + ``tangent``. They are geometry properties, not physics, so a curved + target simply supplies the axes its sweep should use. + """ + + geometry_id: str + reference_point: NDArray[np.float64] + normal: NDArray[np.float64] + tangent: NDArray[np.float64] + components: tuple[ComponentSpec, ...] + + @classmethod + def from_mapping(cls, data: Mapping[str, Any], + default_geometry_id: str) -> "TargetSpec": + normal = _as_vector(data.get("normal", [0.0, 0.0, 1.0]), + "target.normal") + tangent = _as_vector(data.get("tangent", [1.0, 0.0, 0.0]), + "target.tangent") + for label, vec in (("normal", normal), ("tangent", tangent)): + if not np.isfinite(vec).all() or np.linalg.norm(vec) == 0.0: + raise StudyConfigError( + f"target.{label} must be a finite non-zero vector") + normal = normal / np.linalg.norm(normal) + tangent = tangent / np.linalg.norm(tangent) + if abs(float(normal @ tangent)) > 1e-8: + raise StudyConfigError( + "target.normal and target.tangent must be orthogonal " + f"(dot product {float(normal @ tangent):.3g})") + + raw_components = data.get("components") + if raw_components is None: + components: tuple[ComponentSpec, ...] = ( + ComponentSpec(name="target"),) + else: + if not isinstance(raw_components, Sequence) or not raw_components: + raise StudyConfigError( + "target.components must be a non-empty list") + components = tuple(ComponentSpec.from_mapping(entry) + for entry in raw_components) + names = [component.name for component in components] + if len(set(names)) != len(names): + raise StudyConfigError( + f"target.components names must be unique, got {names}") + + return cls( + geometry_id=str(data.get("geometry_id") or default_geometry_id), + reference_point=_as_vector(data.get("reference_point", + [0.0, 0.0, 0.0]), + "target.reference_point"), + normal=normal, + tangent=tangent, + components=components, + ) + + +@dataclass(frozen=True) +class PrescribedFiringSpec: + """One explicitly prescribed firing: pose, active thrusters, duration.""" + + position: NDArray[np.float64] + dcm: NDArray[np.float64] + thrusters: tuple[int, ...] + duration_s: float + + @classmethod + def from_mapping(cls, data: Mapping[str, Any], index: int, + default_thrusters: Sequence[int], + default_duration_s: float) -> "PrescribedFiringSpec": + context = f"firings[{index}]" + position = _as_vector(_require(data, "position", context), + f"{context}.position") + dcm = np.asarray(_require(data, "dcm", context), dtype=float) + if dcm.shape != (3, 3): + raise StudyConfigError( + f"{context}.dcm must be a 3x3 matrix, got shape {dcm.shape}") + thrusters = tuple(int(t) for t in + data.get("thrusters", default_thrusters)) + if not thrusters: + raise StudyConfigError(f"{context}.thrusters must not be empty") + duration = float(data.get("duration_s", default_duration_s)) + if duration <= 0.0: + raise StudyConfigError( + f"{context}.duration_s must be positive, got {duration}") + return cls(position=position, dcm=dcm, thrusters=thrusters, + duration_s=duration) + + +#: Sweep execution modes. ``per_case`` writes one Jet Firing History per +#: angle-distance combination (isolated poses); ``single_jfh`` writes ONE +#: history spanning the whole sweep, so every pose's strikes come from a +#: single pipeline run and the cumulative fields form a sweep envelope. +SWEEP_MODES = ("per_case", "single_jfh") + + +@dataclass(frozen=True) +class SweepSpec: + """Parameter sweep and firing-count definition. + + Attributes + ---------- + plate_angles_deg : tuple of float + Approach angles swept in the target's (normal, tangent) plane; + 0 deg is head-on along the target normal. + source_distances : tuple of float + Plume-source distances from ``TargetSpec.reference_point``, in the + case's length units. + n_firings : int + Number of Jet Firing History entries contributed by EACH pose. In + ``per_case`` mode that is the exact length of every case's history; + in ``single_jfh`` mode the one shared history holds exactly + ``len(poses) * n_firings`` entries. Either way the count is exact. + mode : str + ``'per_case'`` (default) or ``'single_jfh'``; see :data:`SWEEP_MODES`. + firing_duration_s : float + Firing time recorded for each JFH entry. + thrusters : tuple of int + JFH thruster indices active in every generated firing. + firings : tuple of PrescribedFiringSpec + Explicitly prescribed firings, replacing the generated poses. In + ``per_case`` mode their count must equal ``n_firings``; in + ``single_jfh`` mode it must equal ``len(poses) * n_firings``, and + they are assigned to the poses in sweep order. + """ + + plate_angles_deg: tuple[float, ...] + source_distances: tuple[float, ...] + n_firings: int + firing_duration_s: float + thrusters: tuple[int, ...] + firings: tuple[PrescribedFiringSpec, ...] = () + mode: str = "per_case" + + @property + def poses(self) -> tuple[tuple[float, float], ...]: + """(plate angle, source distance) pairs in execution order. + + Distance-major, matching the committed sweep-JFH generators: all + angles at the first distance, then all angles at the next. + """ + return tuple((angle, distance) + for distance in self.source_distances + for angle in self.plate_angles_deg) + + @property + def total_firings(self) -> int: + """Total JFH entries across the whole sweep.""" + return len(self.poses) * self.n_firings + + @classmethod + def from_mapping(cls, data: Mapping[str, Any]) -> "SweepSpec": + angles = tuple(_as_float_list( + data.get("plate_angles_deg", [0.0]), "sweep.plate_angles_deg")) + distances = tuple(_as_float_list( + data.get("source_distances", []), "sweep.source_distances")) + if not angles: + raise StudyConfigError("sweep.plate_angles_deg must not be empty") + if not distances: + raise StudyConfigError("sweep.source_distances must not be empty") + if any(d <= 0.0 for d in distances): + raise StudyConfigError( + "sweep.source_distances must all be positive") + + mode = str(data.get("mode", "per_case")) + if mode not in SWEEP_MODES: + raise StudyConfigError( + f"sweep.mode must be one of {list(SWEEP_MODES)}, got {mode!r}") + + n_firings = validate_n_firings(data.get("n_firings", 1)) + duration = float(data.get("firing_duration_s", 1.0)) + if duration <= 0.0: + raise StudyConfigError( + f"sweep.firing_duration_s must be positive, got {duration}") + thrusters = tuple(int(t) for t in data.get("thrusters", [1])) + if not thrusters: + raise StudyConfigError("sweep.thrusters must not be empty") + + raw_firings = data.get("firings") + firings: tuple[PrescribedFiringSpec, ...] = () + if raw_firings: + firings = tuple( + PrescribedFiringSpec.from_mapping(entry, i, thrusters, duration) + for i, entry in enumerate(raw_firings)) + n_poses = len(angles) * len(distances) + expected = n_firings if mode == "per_case" else n_poses * n_firings + if len(firings) != expected: + where = ("per case" if mode == "per_case" + else f"across {n_poses} poses") + raise StudyConfigError( + f"sweep.n_firings is {n_firings} ({expected} entries " + f"{where} in {mode!r} mode) but {len(firings)} explicit " + "firings were supplied; the counts must agree exactly") + + return cls(plate_angles_deg=angles, source_distances=distances, + n_firings=n_firings, firing_duration_s=duration, + thrusters=thrusters, firings=firings, mode=mode) + + +def validate_n_firings(value: Any) -> int: + """Validate a requested firing count: a positive integer, nothing else. + + ``n_firings`` is the exact number of Jet Firing History entries a case + writes, so zero, negative, fractional and non-numeric values are rejected + rather than coerced. The rule itself lives with the JFH-generation code + (:func:`pyrpod.rpod.approach_maneuvers.validate_n_firings`) so the + prescribed and dynamics-driven paths cannot drift apart; only the raised + error type is specialized here. + """ + try: + return _validate_n_firings(value) + except ValueError as exc: + raise StudyConfigError(str(exc)) from exc + + +@dataclass(frozen=True) +class Normalization: + """Coefficient normalization inputs. + + Coefficients are computed only when the configuration supplies every + value a given coefficient needs; nothing is invented or defaulted: + + * force coefficients need ``reference_area`` and ``dynamic_pressure``; + * moment coefficients additionally need ``reference_length``; + * surface-load coefficients (pressure/shear/heat flux) need + ``dynamic_pressure`` and, for heat flux, ``reference_heat_flux``. + """ + + reference_area: float | None = None + reference_length: float | None = None + dynamic_pressure: float | None = None + reference_heat_flux: float | None = None + + @classmethod + def from_mapping(cls, data: Mapping[str, Any] | None) -> "Normalization": + if not data: + return cls() + + def optional(key: str) -> float | None: + value = data.get(key) + if value is None: + return None + number = float(value) + if not np.isfinite(number) or number <= 0.0: + raise StudyConfigError( + f"loads.normalization.{key} must be a positive finite " + f"number, got {value!r}") + return number + + return cls(reference_area=optional("reference_area"), + reference_length=optional("reference_length"), + dynamic_pressure=optional("dynamic_pressure"), + reference_heat_flux=optional("reference_heat_flux")) + + @property + def has_force_inputs(self) -> bool: + return (self.reference_area is not None + and self.dynamic_pressure is not None) + + @property + def has_moment_inputs(self) -> bool: + return self.has_force_inputs and self.reference_length is not None + + def to_dict(self) -> dict[str, float | None]: + return {"reference_area": self.reference_area, + "reference_length": self.reference_length, + "dynamic_pressure": self.dynamic_pressure, + "reference_heat_flux": self.reference_heat_flux} + + +@dataclass(frozen=True) +class LoadsSpec: + """Surface-load integration settings.""" + + moment_reference_point: NDArray[np.float64] + normalization: Normalization + + @classmethod + def from_mapping(cls, data: Mapping[str, Any] | None) -> "LoadsSpec": + data = data or {} + return cls( + moment_reference_point=_as_vector( + data.get("moment_reference_point", [0.0, 0.0, 0.0]), + "loads.moment_reference_point"), + normalization=Normalization.from_mapping( + data.get("normalization")), + ) + + +@dataclass(frozen=True) +class OutputSpec: + """Output artifact settings (VTK, machine-readable summary, plots).""" + + write_vtk: bool = True + vtk_subdir: str = "vtk" + summary_csv: str = "case_results.csv" + summary_metadata: str = "study_metadata.json" + write_plots: bool = False + plots_subdir: str = "plots" + + @classmethod + def from_mapping(cls, data: Mapping[str, Any] | None) -> "OutputSpec": + data = data or {} + vtk = data.get("vtk") or {} + summary = data.get("summary") or {} + plots = data.get("plots") or {} + return cls( + write_vtk=bool(vtk.get("enabled", True)), + vtk_subdir=str(vtk.get("subdir", "vtk")), + summary_csv=str(summary.get("csv", "case_results.csv")), + summary_metadata=str(summary.get("metadata", + "study_metadata.json")), + write_plots=bool(plots.get("enabled", False)), + plots_subdir=str(plots.get("subdir", "plots")), + ) + + +@dataclass(frozen=True) +class ReferenceSpec: + """Optional external reference-data location (see mdao.reference_data).""" + + path: str | None = None + label: str | None = None + + @classmethod + def from_mapping(cls, data: Mapping[str, Any] | None) -> "ReferenceSpec": + data = data or {} + path = data.get("path") + return cls(path=str(path) if path else None, + label=str(data["label"]) if data.get("label") else None) + + +@dataclass(frozen=True) +class StudyConfig: + """A complete, validated study configuration. + + Attributes + ---------- + study_name : str + Identifier recorded with every result row. + case_dir : str + Existing PyRPOD case directory (owns ``config.ini`` and the STL / TCD + / plume assets). Always ends with a path separator, as the rest of + PyRPOD expects. + output_dir : str + Directory the study writes its artifacts to. + plume_model : str + Always ``'SimplifiedGasKinetics'`` in this branch. + source_path : str + Path the configuration was read from (configuration provenance). + """ + + study_name: str + description: str + case_dir: str + output_dir: str + target: TargetSpec + sweep: SweepSpec + loads: LoadsSpec + output: OutputSpec + reference: ReferenceSpec + thruster_id: str | None = None + plume_model: str = SUPPORTED_PLUME_MODEL + plume_model_parameters: dict[str, Any] = field(default_factory=dict) + coordinate_system: str = "case global frame" + units: dict[str, str] = field(default_factory=lambda: dict(DEFAULT_UNITS)) + source_path: str = "" + + # ------------------------------------------------------------------ load + @classmethod + def from_yaml(cls, path: str | os.PathLike[str]) -> "StudyConfig": + """Read and validate a YAML study configuration.""" + path = os.fspath(path) + if not os.path.isfile(path): + raise StudyConfigError(f"study configuration not found: {path!r}") + with open(path, "r", encoding="utf-8") as handle: + data = yaml.safe_load(handle) + if not isinstance(data, Mapping): + raise StudyConfigError( + f"study configuration {path!r} must be a YAML mapping") + return cls.from_mapping(data, source_path=path) + + @classmethod + def from_mapping(cls, data: Mapping[str, Any], + source_path: str = "") -> "StudyConfig": + """Validate an already-parsed configuration mapping.""" + study = data.get("study") or {} + if not isinstance(study, Mapping): + raise StudyConfigError("'study' section must be a mapping") + + name = str(_require(study, "name", "study")) + raw_case_dir = str(_require(study, "case_dir", "study")) + base_dir = os.path.dirname(os.path.abspath(source_path)) if \ + source_path else os.getcwd() + case_dir = _resolve_dir(raw_case_dir, base_dir) + if not os.path.isdir(case_dir): + raise StudyConfigError( + f"study.case_dir does not exist: {case_dir!r}") + if not os.path.isfile(os.path.join(case_dir, "config.ini")): + raise StudyConfigError( + f"study.case_dir {case_dir!r} has no config.ini; a study " + "always runs on top of an existing PyRPOD case") + case_dir = case_dir.rstrip("\\/") + os.sep + + raw_output_dir = str(study.get( + "output_dir", os.path.join(case_dir, "results", "studies", name))) + output_dir = _resolve_dir(raw_output_dir, base_dir) + + plume = data.get("plume_model") or {} + model_name = str(plume.get("name", SUPPORTED_PLUME_MODEL)) + if model_name != SUPPORTED_PLUME_MODEL: + raise StudyConfigError( + f"plume_model.name must be {SUPPORTED_PLUME_MODEL!r} " + f"(this study workflow supports exactly one model), got " + f"{model_name!r}") + model_parameters = dict(plume.get("parameters") or {}) + + thruster = data.get("thruster") or {} + thruster_id = thruster.get("id") + + target_data = data.get("target") or {} + default_geometry = _default_geometry_id(case_dir) + + metadata = data.get("metadata") or {} + units = dict(DEFAULT_UNITS) + units.update({str(k): str(v) + for k, v in (metadata.get("units") or {}).items()}) + + return cls( + study_name=name, + description=str(study.get("description", "")), + case_dir=case_dir, + output_dir=output_dir, + target=TargetSpec.from_mapping(target_data, default_geometry), + sweep=SweepSpec.from_mapping(data.get("sweep") or {}), + loads=LoadsSpec.from_mapping(data.get("loads")), + output=OutputSpec.from_mapping(data.get("output")), + reference=ReferenceSpec.from_mapping(data.get("reference")), + thruster_id=str(thruster_id) if thruster_id else None, + plume_model=model_name, + plume_model_parameters=model_parameters, + coordinate_system=str(metadata.get("coordinate_system", + "case global frame")), + units=units, + source_path=os.path.abspath(source_path) if source_path else "", + ) + + # ------------------------------------------------------------- accessors + @property + def n_cases(self) -> int: + """Number of angle x distance cases in this study.""" + return (len(self.sweep.plate_angles_deg) + * len(self.sweep.source_distances)) + + def with_output_dir(self, output_dir: str) -> "StudyConfig": + """Copy of this configuration writing to a different directory.""" + return replace(self, output_dir=os.path.abspath(output_dir)) + + def provenance(self) -> dict[str, Any]: + """Configuration provenance recorded in the study metadata.""" + return { + "study_name": self.study_name, + "config_path": self.source_path, + "case_dir": os.path.abspath(self.case_dir), + "output_dir": os.path.abspath(self.output_dir), + "plume_model": self.plume_model, + "plume_model_parameters": dict(self.plume_model_parameters), + "coordinate_system": self.coordinate_system, + "units": dict(self.units), + } + + +def _resolve_dir(path: str, base_dir: str) -> str: + """Resolve a configured directory relative to the config file location.""" + if os.path.isabs(path): + return os.path.abspath(path) + return os.path.abspath(os.path.join(base_dir, path)) + + +def _default_geometry_id(case_dir: str) -> str: + """Target geometry id taken from the case's own ``[tv] stl`` entry.""" + import configparser + + config = configparser.ConfigParser() + config.read(os.path.join(case_dir, "config.ini")) + if config.has_option("tv", "stl"): + return str(config["tv"]["stl"]) + return "unknown" From 0a2fe4f6e0a7b346b50148a1e47c0a747c570aa8 Mon Sep 17 00:00:00 2001 From: andytorrestb Date: Mon, 27 Jul 2026 22:53:38 -0500 Subject: [PATCH 03/14] mdao: add prescribed firing-plan generation Validation sweeps do not fly a trajectory: the firing poses are prescribed, so this module owns their generation and the Jet Firing History they are written to, leaving the dynamics-based approach-maneuver workflow untouched. The pose convention -- source on the arc of radius L about the target reference point, aimed at it, DCM first column carrying the thruster axis and second column the sweep plane's binormal -- reproduces the committed sweep JFH generators to file precision, for both the flat and inclined cases. Two builders: one pose's firings (exactly n_firings entries), and the whole sweep as one continuous sequence (exactly len(poses) * n_firings entries, firing times running through). Each firing records the pose it realizes, so a history spanning many poses stays keyed to the sweep grid. Co-Authored-By: Claude Opus 5 --- pyrpod/mdao/firing_plan.py | 290 +++++++++++++++++++++++++++++++++++++ 1 file changed, 290 insertions(+) create mode 100644 pyrpod/mdao/firing_plan.py diff --git a/pyrpod/mdao/firing_plan.py b/pyrpod/mdao/firing_plan.py new file mode 100644 index 0000000..28c91ab --- /dev/null +++ b/pyrpod/mdao/firing_plan.py @@ -0,0 +1,290 @@ +""" +Prescribed Jet Firing History generation for plume/target validation studies. + +Validation sweeps do not fly a trajectory: the firing poses are PRESCRIBED +(placed by the engineer, or generated from a swept approach angle and source +distance about a stationary target) and the Jet Firing History is written +from them directly. This module owns that generation, so the dynamics-based +approach-maneuver workflow in ``pyrpod.rpod.approach_maneuvers`` stays +untouched. + +``n_firings`` semantics +----------------------- +``n_firings`` is the EXACT number of entries written to the JFH. The count is +validated before anything is generated (:func:`validate_n_firings`), the +generated sequence is asserted to match it (:func:`build_case_firings`), and +an explicitly prescribed firing list whose length disagrees is an error -- +never a silent truncation or extension. + +Pose convention +--------------- +For a target reference point ``C`` with outward normal ``n_hat`` (pointing +toward the plume source side) and in-plane tangent ``t_hat``: + + d_hat(alpha) = cos(alpha) * n_hat + sin(alpha) * t_hat + source position = C + L * d_hat(alpha) + thruster axis = -d_hat(alpha) (aimed at C) + +``alpha = 0`` is head-on. The JFH DCM is built with the thruster axis as its +first COLUMN, which is what the strike pipeline reads as the plume normal +(``dcm.T`` rows, with an identity thruster DCM). This reproduces the pose +convention of the committed sweep-JFH generators +(``case/plume/plume_flat_plate_sweep/jfh/generate_sweep_jfh.py``) exactly, +including the binormal choice ``cross(n_hat, t_hat)`` for the DCM's second +column. +""" + +from __future__ import annotations + +import os +from dataclasses import dataclass +from typing import Iterable, Sequence + +import numpy as np +from numpy.typing import NDArray + +from pyrpod.mdao.study_config import ( + PrescribedFiringSpec, + StudyConfigError, + SweepSpec, + TargetSpec, + validate_n_firings, +) + +__all__ = [ + "Firing", + "build_case_firings", + "build_sweep_firings", + "pose_for", + "validate_n_firings", + "write_jfh_file", +] + + +@dataclass(frozen=True) +class Firing: + """One JFH entry: pose, active thrusters, firing duration and start time. + + Attributes + ---------- + position : np.ndarray + Visiting-vehicle (plume-source) position in the case's global frame. + dcm : np.ndarray + 3x3 direction cosine matrix; its first column is the thruster axis. + thrusters : tuple of int + JFH thruster indices active during this firing (1-based). + duration_s : float + Firing time; the strike pipeline multiplies heat flux by this to get + the per-firing heat-flux load. + start_time_s : float + Elapsed time at the start of this firing (informational; the strike + calculation does not read it). + plate_angle_deg, source_distance : float or None + The swept parameters this firing realizes. Set when the firing came + from a sweep pose, so a history spanning many poses stays keyed to + the sweep grid; None for firings whose pose was prescribed outright + without a sweep parameterization. + pose_index : int or None + Index of the firing's pose in the sweep's execution order. + """ + + position: NDArray[np.float64] + dcm: NDArray[np.float64] + thrusters: tuple[int, ...] + duration_s: float = 1.0 + start_time_s: float = 0.0 + plate_angle_deg: float | None = None + source_distance: float | None = None + pose_index: int | None = None + + +def pose_for(alpha_deg: float, distance: float, + reference_point: Sequence[float] | NDArray[np.float64], + normal: Sequence[float] | NDArray[np.float64], + tangent: Sequence[float] | NDArray[np.float64], + ) -> tuple[NDArray[np.float64], NDArray[np.float64]]: + """Plume-source pose for one (approach angle, distance) combination. + + Parameters + ---------- + alpha_deg : float + Approach angle in degrees; 0 is head-on along ``normal``, positive + angles rotate the source toward ``tangent``. + distance : float + Distance from ``reference_point`` to the plume source. + reference_point, normal, tangent : array-like + Target geometry axes (see :class:`pyrpod.mdao.study_config.TargetSpec`). + + Returns + ------- + (np.ndarray, np.ndarray) + The source position and the 3x3 DCM whose first column is the + thruster axis (aimed back at ``reference_point``). + """ + center = np.asarray(reference_point, dtype=float) + n_hat = np.asarray(normal, dtype=float) + t_hat = np.asarray(tangent, dtype=float) + n_hat = n_hat / np.linalg.norm(n_hat) + t_hat = t_hat / np.linalg.norm(t_hat) + + alpha = np.deg2rad(float(alpha_deg)) + d_hat = np.cos(alpha) * n_hat + np.sin(alpha) * t_hat + position = center + float(distance) * d_hat + axis = -d_hat # aimed at the target reference point + + # Right-handed triad completed with the sweep plane's binormal, matching + # the committed sweep-JFH generators (which hard-code the global Y axis + # for their X-Z sweep plane -- the same vector this expression yields). + binormal = np.cross(n_hat, t_hat) + binormal = binormal / np.linalg.norm(binormal) + dcm = np.column_stack([axis, binormal, np.cross(axis, binormal)]) + return position, dcm + + +def build_case_firings(sweep: SweepSpec, target: TargetSpec, + alpha_deg: float, distance: float, + pose_index: int = 0, + start_time_s: float = 0.0) -> list[Firing]: + """Build EXACTLY ``sweep.n_firings`` firings for one sweep pose. + + When the configuration prescribes firings explicitly they are used + verbatim (their count is validated against ``n_firings`` when the + configuration is parsed). Otherwise the swept pose is generated from + ``alpha_deg`` / ``distance`` and repeated for ``n_firings`` successive + firing intervals -- one JFH entry per requested firing. + + ``pose_index`` and ``start_time_s`` place this pose inside a longer + sequence; they matter only when many poses share one history (see + :func:`build_sweep_firings`). + + Raises + ------ + StudyConfigError + If the generated sequence length would differ from ``n_firings``. + """ + n_firings = validate_n_firings(sweep.n_firings) + + if sweep.firings: + specs = sweep.firings + if sweep.mode == "single_jfh": + # One shared history: each pose takes its own contiguous slice of + # the prescribed sequence, in sweep order. + offset = pose_index * n_firings + specs = sweep.firings[offset:offset + n_firings] + firings = [ + _from_spec(spec, index, alpha_deg, distance, pose_index, + start_time_s) + for index, spec in enumerate(specs) + ] + else: + position, dcm = pose_for(alpha_deg, distance, + target.reference_point, target.normal, + target.tangent) + firings = [ + Firing(position=position, dcm=dcm, thrusters=sweep.thrusters, + duration_s=sweep.firing_duration_s, + start_time_s=start_time_s + + index * sweep.firing_duration_s, + plate_angle_deg=float(alpha_deg), + source_distance=float(distance), + pose_index=pose_index) + for index in range(n_firings) + ] + + if len(firings) != n_firings: + raise StudyConfigError( + f"requested n_firings={n_firings} but built {len(firings)} " + "firings; the two must agree exactly") + return firings + + +def build_sweep_firings(sweep: SweepSpec, + target: TargetSpec) -> list[Firing]: + """Build the WHOLE sweep as one firing sequence. + + Every pose of ``sweep.poses`` contributes ``sweep.n_firings`` entries, in + execution order, with firing times running continuously across the + sequence. The result is the single Jet Firing History of a ``single_jfh`` + study, and its length is exactly ``sweep.total_firings``. + + Raises + ------ + StudyConfigError + If the assembled sequence length would differ from that product. + """ + firings: list[Firing] = [] + elapsed = 0.0 + for pose_index, (angle, distance) in enumerate(sweep.poses): + pose_firings = build_case_firings(sweep, target, angle, distance, + pose_index=pose_index, + start_time_s=elapsed) + firings.extend(pose_firings) + elapsed += sum(firing.duration_s for firing in pose_firings) + + if len(firings) != sweep.total_firings: + raise StudyConfigError( + f"sweep of {len(sweep.poses)} poses x n_firings=" + f"{sweep.n_firings} must produce {sweep.total_firings} JFH " + f"entries, built {len(firings)}") + return firings + + +def _from_spec(spec: PrescribedFiringSpec, index: int, + alpha_deg: float | None = None, + distance: float | None = None, + pose_index: int | None = None, + start_time_s: float = 0.0) -> Firing: + return Firing(position=np.asarray(spec.position, dtype=float), + dcm=np.asarray(spec.dcm, dtype=float), + thrusters=spec.thrusters, + duration_s=spec.duration_s, + start_time_s=start_time_s + index * spec.duration_s, + plate_angle_deg=(None if alpha_deg is None + else float(alpha_deg)), + source_distance=(None if distance is None + else float(distance)), + pose_index=pose_index) + + +def write_jfh_file(path: str | os.PathLike[str], + firings: Iterable[Firing]) -> int: + """Write firings to a JFH file readable by ``JetFiringHistory.read_jfh``. + + The emitted format matches the committed case generators byte-for-byte in + structure (header line with the firing count, an unused second line, then + one row per firing: index, dt, t, unused column, nine DCM values, three + position values, uncertainty factor, thruster count, thruster indices). + + Returns + ------- + int + The number of entries written, which is the number of firings the + JFH will report. + """ + firings = list(firings) + if not firings: + raise StudyConfigError( + "refusing to write a JFH with zero firings; n_firings must be a " + "positive integer") + + path = os.fspath(path) + parent = os.path.dirname(path) + if parent: + os.makedirs(parent, exist_ok=True) + + lines = [f"offseted {len(firings)} 0", + " 0.000 0.000 0.000"] + for index, firing in enumerate(firings, start=1): + dcm = np.asarray(firing.dcm, dtype=float).reshape(3, 3) + dcm_str = " ".join(f"{value:.6e}" for value in dcm.ravel()) + xyz_str = " ".join(f"{value:.9g}" + for value in np.asarray(firing.position, + dtype=float).ravel()) + thruster_str = " ".join(str(int(t)) for t in firing.thrusters) + lines.append( + f" {index} {firing.start_time_s:g} {firing.duration_s:g} 0 " + f"{dcm_str} {xyz_str} 1 {len(firing.thrusters)} {thruster_str}") + + with open(path, "w", newline="\n", encoding="utf-8") as handle: + handle.write("\n".join(lines) + "\n") + return len(firings) From 006130503c5746b3317547d0d7d8d29524ef2625 Mon Sep 17 00:00:00 2001 From: andytorrestb Date: Mon, 27 Jul 2026 22:53:58 -0500 Subject: [PATCH 04/14] mdao: add per-component surface-load integration Turns the strike pipeline's per-face pressure, shear and heat-flux arrays into the integrated quantities a validation study reports: pressure, shear and combined force vectors, the corresponding moments about a user-defined reference point, center of pressure, peak surface loads, affected area, component heat load and nondimensional coefficients. Per-face data are not reduced away here -- they keep flowing to the VTK writer unchanged. Conventions are stated explicitly and follow the pipeline's own face-selection rule: pressure acts along -n_hat, shear along the tangential projection of the radial flow direction. The center of pressure is defined as the point on the resultant's line of action closest to the moment reference, since a unique 3D center of pressure does not exist in general; the irreducible force-parallel couple is reported separately, and zero-load and near-cancellation cases return no center of pressure with a status naming the reason rather than a misleadingly large value. An always-defined pressure-weighted centroid is reported alongside. Coefficients are computed only for the normalizations whose inputs are supplied in full; nothing is defaulted or invented. Co-Authored-By: Claude Opus 5 --- pyrpod/mdao/surface_loads.py | 386 +++++++++++++++++++++++++++++++++++ 1 file changed, 386 insertions(+) create mode 100644 pyrpod/mdao/surface_loads.py diff --git a/pyrpod/mdao/surface_loads.py b/pyrpod/mdao/surface_loads.py new file mode 100644 index 0000000..0d7767c --- /dev/null +++ b/pyrpod/mdao/surface_loads.py @@ -0,0 +1,386 @@ +""" +Integration of per-face plume-impingement fields into component-level loads. + +The strike pipeline (:mod:`pyrpod.plume.PlumeStrikeCalculator`) produces +per-face pressure, shear-stress and heat-flux arrays. This module turns those +arrays into the integrated quantities an engineer reports for a target +component -- force, moment, center of pressure, peak surface loads, affected +area and (only when normalization inputs are supplied) nondimensional +coefficients. Per-face data are never discarded here; they keep flowing to +the VTK writer unchanged. + +Sign and direction conventions +------------------------------ +These follow the strike pipeline's own face-selection convention, in which a +face is struck only when its mesh unit normal ``n_hat`` opposes the plume +direction (``n_hat . u_plume < 0``), i.e. ``n_hat`` points back toward the +plume source: + +* **Pressure** acts INTO the surface, along ``-n_hat``: + ``F_p,i = -p_i * A_i * n_hat_i``. +* **Shear** acts along the tangential projection of the local flow direction + ``u_hat_i`` (the collisionless flow is radial, so ``u_hat_i`` is the unit + vector from the plume source to the face centroid): + ``t_hat_i = normalize(u_hat_i - (u_hat_i . n_in_i) * n_in_i)``, + ``n_in_i = -n_hat_i``, and ``F_s,i = tau_i * A_i * t_hat_i``. + Faces at normal incidence have no tangential direction and contribute no + shear force (their shear magnitude is zero there anyway). +* **Moment** about the user-supplied reference point ``r_ref``: + ``M_ref = sum (r_i - r_ref) x F_i``. + +Center of pressure +------------------ +A unique three-dimensional center of pressure does not exist in general: the +resultant of a distributed load is a force plus a couple, and only the +component of the moment perpendicular to the force can be represented by +shifting the force's line of action. The definition used here is stated +explicitly: + + ``r_cop = r_ref + (F x M_ref) / |F|^2`` + +which is the point ON THE LINE OF ACTION of the resultant force that lies +CLOSEST to the moment reference point. The moment component parallel to the +force is reported separately as ``residual_couple`` -- it cannot be removed +by any choice of ``r_cop``. When the resultant force is zero or numerically +insignificant against the summed face-force magnitudes (near-total +cancellation, e.g. a closed target loaded symmetrically), no center of +pressure is returned and the status field records why; no misleadingly large +value is ever produced. + +The auxiliary ``pressure_weighted_centroid``, +``sum(p_i A_i r_i) / sum(p_i A_i)``, is always defined for a non-zero +pressure load and coincides with the classical center of pressure for a +planar component under unidirectional pressure. +""" + +from __future__ import annotations + +from dataclasses import dataclass, field +from typing import Any, Mapping, Sequence + +import numpy as np +from numpy.typing import NDArray + +from pyrpod.mdao.study_config import ComponentSpec, Normalization + +__all__ = [ + "ComponentLoads", + "face_areas", + "flow_directions", + "integrate_component_loads", + "select_component_faces", +] + +#: Relative floor below which a resultant force is treated as cancellation +#: noise rather than a meaningful resultant (see the module docstring). +FORCE_SIGNIFICANCE_RATIO = 1e-6 + + +def face_areas(vectors: NDArray[np.float64]) -> NDArray[np.float64]: + """Per-face triangle areas for an (N, 3, 3) array of face vertices.""" + vectors = np.asarray(vectors, dtype=float) + v0, v1, v2 = vectors[:, 0], vectors[:, 1], vectors[:, 2] + return 0.5 * np.linalg.norm(np.cross(v1 - v0, v2 - v0), axis=1) + + +def flow_directions(centroids: NDArray[np.float64], + source_position: Sequence[float] | NDArray[np.float64], + ) -> NDArray[np.float64]: + """Unit vectors from the plume source to each face centroid. + + The collisionless plume flow is radial from the nozzle exit, so this is + the local flow direction at each face -- the same quantity the strike + calculator uses for the true incidence angle. + """ + centroids = np.asarray(centroids, dtype=float) + relative = centroids - np.asarray(source_position, dtype=float).reshape(3) + norms = np.linalg.norm(relative, axis=1) + directions = np.zeros_like(relative) + np.divide(relative, norms[:, None], out=directions, + where=(norms > 0.0)[:, None]) + return directions + + +def select_component_faces(component: ComponentSpec, + centroids: NDArray[np.float64], + n_faces: int) -> NDArray[np.int64]: + """Face indices belonging to one target component. + + Raises + ------ + ValueError + If an explicit face index is out of range, or a bounds selector + matches no face (a silently empty component would produce + meaningless zero loads). + """ + if component.selector == "all": + return np.arange(n_faces, dtype=np.int64) + + if component.selector == "face_indices": + indices = np.asarray(component.face_indices, dtype=np.int64) + if indices.min() < 0 or indices.max() >= n_faces: + raise ValueError( + f"component {component.name!r}: face index out of range for a " + f"{n_faces}-face mesh") + return indices + + lo = np.asarray(component.bounds["min"], dtype=float) + hi = np.asarray(component.bounds["max"], dtype=float) + inside = np.all((centroids >= lo) & (centroids <= hi), axis=1) + indices = np.nonzero(inside)[0].astype(np.int64) + if indices.size == 0: + raise ValueError( + f"component {component.name!r}: bounds selector matched no face " + "of the target mesh") + return indices + + +@dataclass(frozen=True) +class ComponentLoads: + """Integrated loads for one target component and one firing. + + All vectors are in the case's global frame and SI units (N, N*m, m). + ``coefficients`` is empty when the configuration supplies no (or + incomplete) normalization inputs -- coefficients are never invented. + """ + + component: str + n_faces: int + n_struck_faces: int + total_area: float + affected_area: float + + pressure_force: NDArray[np.float64] + shear_force: NDArray[np.float64] + force: NDArray[np.float64] + force_magnitude: float + + moment_reference_point: NDArray[np.float64] + pressure_moment: NDArray[np.float64] + shear_moment: NDArray[np.float64] + moment: NDArray[np.float64] + moment_magnitude: float + + center_of_pressure: NDArray[np.float64] | None + center_of_pressure_status: str + residual_couple: float + pressure_weighted_centroid: NDArray[np.float64] | None + + max_pressure: float + max_shear_stress: float + max_heat_flux: float + total_heat_load: float + + coefficients: dict[str, float] = field(default_factory=dict) + + @property + def has_coefficients(self) -> bool: + return bool(self.coefficients) + + +def integrate_component_loads( + *, + component_name: str, + face_indices: NDArray[np.int64] | Sequence[int], + centroids: NDArray[np.float64], + unit_normals: NDArray[np.float64], + areas: NDArray[np.float64], + pressures: NDArray[np.float64], + shear_stresses: NDArray[np.float64], + heat_fluxes: NDArray[np.float64], + strikes: NDArray[np.float64] | None, + moment_reference_point: Sequence[float] | NDArray[np.float64], + source_position: Sequence[float] | NDArray[np.float64] | None = None, + flow_unit_vectors: NDArray[np.float64] | None = None, + normalization: Normalization | None = None, +) -> ComponentLoads: + """Integrate per-face fields into one component's resultant loads. + + Parameters + ---------- + component_name : str + Reported component identifier. + face_indices : array-like of int + Indices of the faces belonging to this component (see + :func:`select_component_faces`). + centroids, unit_normals, areas : np.ndarray + Per-face geometry of the FULL target mesh, indexed by ``face_indices``. + pressures, shear_stresses, heat_fluxes : np.ndarray + Per-face load fields of the full mesh for one firing (Pa, Pa, W/m^2). + strikes : np.ndarray or None + Per-face strike counts, used for the affected-area tally. When None, + faces carrying any non-zero load are counted instead. + moment_reference_point : array-like + User-defined point the moment is taken about. + source_position : array-like, optional + Plume-source position, used to derive the local flow directions the + shear force acts along. Required unless ``flow_unit_vectors`` is given. + flow_unit_vectors : np.ndarray, optional + Precomputed per-face flow directions (see :func:`flow_directions`). + normalization : Normalization, optional + Coefficient normalization inputs. Coefficients are computed only for + the ones whose inputs are complete. + + Returns + ------- + ComponentLoads + """ + indices = np.asarray(face_indices, dtype=np.int64) + if indices.size == 0: + raise ValueError( + f"component {component_name!r} selects no faces; refusing to " + "report zero loads for an empty component") + + centroid = np.asarray(centroids, dtype=float)[indices] + normal = np.asarray(unit_normals, dtype=float)[indices] + area = np.asarray(areas, dtype=float)[indices] + pressure = np.asarray(pressures, dtype=float)[indices] + shear = np.asarray(shear_stresses, dtype=float)[indices] + heat_flux = np.asarray(heat_fluxes, dtype=float)[indices] + + if flow_unit_vectors is None: + if source_position is None: + raise ValueError( + "integrate_component_loads needs either source_position or " + "flow_unit_vectors to orient the shear force") + flow = flow_directions(centroid, source_position) + else: + flow = np.asarray(flow_unit_vectors, dtype=float)[indices] + + # --- forces ----------------------------------------------------------- + # Pressure pushes into the surface (along -n_hat, the plume direction for + # every struck face); shear drags along the tangential flow projection. + pressure_face_force = -(pressure * area)[:, None] * normal + + inward = -normal + cos_incidence = np.einsum("ij,ij->i", flow, inward) + tangential = flow - cos_incidence[:, None] * inward + tangential_norm = np.linalg.norm(tangential, axis=1) + tangent_hat = np.zeros_like(tangential) + np.divide(tangential, tangential_norm[:, None], out=tangent_hat, + where=(tangential_norm > 1e-12)[:, None]) + shear_face_force = (shear * area)[:, None] * tangent_hat + + pressure_force = pressure_face_force.sum(axis=0) + shear_force = shear_face_force.sum(axis=0) + total_force = pressure_force + shear_force + + # --- moments ---------------------------------------------------------- + reference = np.asarray(moment_reference_point, dtype=float).reshape(3) + arm = centroid - reference + pressure_moment = np.cross(arm, pressure_face_force).sum(axis=0) + shear_moment = np.cross(arm, shear_face_force).sum(axis=0) + total_moment = pressure_moment + shear_moment + + # --- center of pressure ---------------------------------------------- + face_force_scale = float(np.sum(np.linalg.norm( + pressure_face_force + shear_face_force, axis=1))) + force_magnitude = float(np.linalg.norm(total_force)) + cop, cop_status, residual_couple = _center_of_pressure( + total_force, total_moment, reference, force_magnitude, + face_force_scale) + + pressure_weight = pressure * area + weight_sum = float(np.sum(pressure_weight)) + pressure_centroid = ( + (pressure_weight[:, None] * centroid).sum(axis=0) / weight_sum + if weight_sum > 0.0 else None) + + # --- surface-field statistics ---------------------------------------- + if strikes is None: + affected = (pressure != 0.0) | (shear != 0.0) | (heat_flux != 0.0) + else: + affected = np.asarray(strikes, dtype=float)[indices] > 0.0 + + loads = ComponentLoads( + component=component_name, + n_faces=int(indices.size), + n_struck_faces=int(np.count_nonzero(affected)), + total_area=float(np.sum(area)), + affected_area=float(np.sum(area[affected])), + pressure_force=pressure_force, + shear_force=shear_force, + force=total_force, + force_magnitude=force_magnitude, + moment_reference_point=reference, + pressure_moment=pressure_moment, + shear_moment=shear_moment, + moment=total_moment, + moment_magnitude=float(np.linalg.norm(total_moment)), + center_of_pressure=cop, + center_of_pressure_status=cop_status, + residual_couple=residual_couple, + pressure_weighted_centroid=pressure_centroid, + max_pressure=float(np.max(pressure)) if pressure.size else 0.0, + max_shear_stress=float(np.max(shear)) if shear.size else 0.0, + max_heat_flux=float(np.max(heat_flux)) if heat_flux.size else 0.0, + total_heat_load=float(np.sum(heat_flux * area)), + coefficients={}, + ) + return _with_coefficients(loads, normalization) + + +def _center_of_pressure( + force: NDArray[np.float64], moment: NDArray[np.float64], + reference: NDArray[np.float64], force_magnitude: float, + face_force_scale: float, +) -> tuple[NDArray[np.float64] | None, str, float]: + """Line-of-action center of pressure, or an explicit reason there is none. + + See the module docstring for the definition. Returns + ``(point_or_None, status, residual_couple)`` where ``status`` is one of + ``'ok'``, ``'zero_load'`` (nothing is loaded) or ``'ill_conditioned'`` + (the resultant is cancellation noise against the summed face forces). + """ + if force_magnitude == 0.0: + return None, "zero_load" if face_force_scale == 0.0 else \ + "ill_conditioned", float(np.linalg.norm(moment)) + if face_force_scale > 0.0 and ( + force_magnitude < FORCE_SIGNIFICANCE_RATIO * face_force_scale): + return None, "ill_conditioned", float(np.linalg.norm(moment)) + + direction = force / force_magnitude + residual_couple = float(np.dot(moment, direction)) + point = reference + np.cross(force, moment) / force_magnitude ** 2 + return point, "ok", abs(residual_couple) + + +def _with_coefficients(loads: ComponentLoads, + normalization: Normalization | None) -> ComponentLoads: + """Attach nondimensional coefficients, but only the fully normalizable ones.""" + if normalization is None: + return loads + + coefficients: dict[str, float] = {} + q_dyn = normalization.dynamic_pressure + area_ref = normalization.reference_area + length_ref = normalization.reference_length + + if q_dyn is not None: + coefficients["Cp_max"] = loads.max_pressure / q_dyn + coefficients["Cf_max"] = loads.max_shear_stress / q_dyn + if normalization.reference_heat_flux is not None: + coefficients["Cq_max"] = (loads.max_heat_flux + / normalization.reference_heat_flux) + if normalization.has_force_inputs: + # mypy: has_force_inputs guarantees both are set + denominator = float(q_dyn) * float(area_ref) # type: ignore[arg-type] + coefficients["CF"] = loads.force_magnitude / denominator + coefficients["CFx"] = float(loads.force[0]) / denominator + coefficients["CFy"] = float(loads.force[1]) / denominator + coefficients["CFz"] = float(loads.force[2]) / denominator + if normalization.has_moment_inputs: + denominator = (float(q_dyn) * float(area_ref) # type: ignore[arg-type] + * float(length_ref)) # type: ignore[arg-type] + coefficients["CM"] = loads.moment_magnitude / denominator + coefficients["CMx"] = float(loads.moment[0]) / denominator + coefficients["CMy"] = float(loads.moment[1]) / denominator + coefficients["CMz"] = float(loads.moment[2]) / denominator + + if not coefficients: + return loads + return ComponentLoads(**{**_as_kwargs(loads), "coefficients": coefficients}) + + +def _as_kwargs(loads: ComponentLoads) -> Mapping[str, Any]: + return {f: getattr(loads, f) for f in loads.__dataclass_fields__} From 102ac21330bc412f4b1c4bf6bbb49794457abacc Mon Sep 17 00:00:00 2001 From: andytorrestb Date: Mon, 27 Jul 2026 22:53:58 -0500 Subject: [PATCH 05/14] mdao: add the structured result schema One record per case, component and firing, carrying everything needed to reproduce and later compare the calculation: identifiers, geometry and mesh information, coordinate system and units, the prescribed source pose, the swept parameters, the plume model and its parameters, the integrated loads, the surface-field peaks, the optional coefficients, the artifact paths and the code/configuration provenance. Two machine-readable artifacts, both in formats the repository already uses and with no new dependency: a flat CSV (vectors expanded to _x/_y/_z columns, one column per coefficient, so it is directly plottable) and a JSON metadata document carrying study-level provenance plus the nested per-case records -- the exchange format an externally generated dataset is transformed into for comparison. Co-Authored-By: Claude Opus 5 --- pyrpod/mdao/study_results.py | 335 +++++++++++++++++++++++++++++++++++ 1 file changed, 335 insertions(+) create mode 100644 pyrpod/mdao/study_results.py diff --git a/pyrpod/mdao/study_results.py b/pyrpod/mdao/study_results.py new file mode 100644 index 0000000..aefce4e --- /dev/null +++ b/pyrpod/mdao/study_results.py @@ -0,0 +1,335 @@ +""" +Structured results for prescribed plume/target validation studies. + +One :class:`CaseResult` records everything needed to REPRODUCE and later +COMPARE a single (case, component, firing) evaluation: the study/case +identity, the geometry and mesh it ran on, the coordinate system and units, +the prescribed plume-source pose, the swept parameters, the plume model and +its parameters, the integrated loads, the surface-field peaks, the optional +coefficients, the VTK artifact path, and the code/configuration provenance. + +Two machine-readable artifacts are written, both in formats already used by +the repository and with no new dependency: + +* a flat CSV of every case row (one row per case x component x firing), for + spreadsheets, plotting and quick diffing; +* a JSON metadata document carrying the study-level provenance plus the + nested per-case records, which is the exchange format an externally + generated dataset (DSMC, experiment, another code) is later transformed + into for :mod:`pyrpod.mdao.reference_data`. +""" + +from __future__ import annotations + +import csv +import json +import os +import subprocess +from dataclasses import asdict, dataclass, field +from datetime import datetime, timezone +from typing import Any, Iterable, Sequence + +import numpy as np +from numpy.typing import NDArray + +from pyrpod.mdao.surface_loads import ComponentLoads + +__all__ = ["CaseResult", "StudyResults", "code_version"] + + +def code_version(repo_dir: str | os.PathLike[str] | None = None) -> str: + """Short git commit of the working tree, or ``'unknown'``. + + Provenance is best-effort: a study must still run in an exported + source tree with no git metadata available. + """ + try: + completed = subprocess.run( + ["git", "rev-parse", "--short", "HEAD"], + cwd=os.fspath(repo_dir) if repo_dir else None, + capture_output=True, text=True, timeout=10, check=False) + except (OSError, subprocess.SubprocessError): + return "unknown" + if completed.returncode != 0: + return "unknown" + return completed.stdout.strip() or "unknown" + + +def _vector(value: Sequence[float] | NDArray[np.float64] | None + ) -> list[float] | None: + if value is None: + return None + return [float(v) for v in np.asarray(value, dtype=float).reshape(-1)] + + +@dataclass +class CaseResult: + """Result record for one case, component and firing. + + Attributes are grouped as: identity, geometry, pose/sweep parameters, + model, integrated loads, surface-field statistics, coefficients, and + artifacts/provenance. Every field is plain data (str, float, list) so the + record serializes to CSV and JSON without custom encoders. + """ + + # --- identity --------------------------------------------------------- + study_name: str + case_id: str + component: str + firing_id: int + + # --- geometry --------------------------------------------------------- + geometry_id: str + mesh_faces: int + component_faces: int + component_area: float + coordinate_system: str + units: dict[str, str] + + # --- prescribed pose and swept parameters ----------------------------- + plume_source_position: list[float] + plume_source_orientation: list[float] + target_normal: list[float] + target_tangent: list[float] + target_reference_point: list[float] + plate_angle_deg: float + source_distance: float + firing_duration_s: float + thrusters: list[int] + + # --- model ------------------------------------------------------------ + plume_model: str + plume_model_parameters: dict[str, Any] + + # --- integrated loads ------------------------------------------------- + pressure_force: list[float] + shear_force: list[float] + force: list[float] + force_magnitude: float + moment_reference_point: list[float] + pressure_moment: list[float] + shear_moment: list[float] + moment: list[float] + moment_magnitude: float + center_of_pressure: list[float] | None + center_of_pressure_status: str + residual_couple: float + pressure_weighted_centroid: list[float] | None + + # --- surface-field statistics ---------------------------------------- + max_pressure: float + max_shear_stress: float + max_heat_flux: float + total_heat_load: float + affected_area: float + struck_faces: int + + # --- optional coefficients ------------------------------------------- + coefficients: dict[str, float] = field(default_factory=dict) + coefficients_available: bool = False + + # --- artifacts and provenance ---------------------------------------- + vtk_path: str | None = None + jfh_path: str | None = None + config_path: str = "" + case_dir: str = "" + code_version: str = "unknown" + generated_at: str = "" + + # ------------------------------------------------------------ builders + @classmethod + def from_loads(cls, loads: ComponentLoads, **metadata: Any) -> "CaseResult": + """Build a result record from integrated loads plus study metadata.""" + return cls( + component=loads.component, + component_faces=loads.n_faces, + component_area=loads.total_area, + pressure_force=_vector(loads.pressure_force) or [], + shear_force=_vector(loads.shear_force) or [], + force=_vector(loads.force) or [], + force_magnitude=loads.force_magnitude, + moment_reference_point=_vector(loads.moment_reference_point) or [], + pressure_moment=_vector(loads.pressure_moment) or [], + shear_moment=_vector(loads.shear_moment) or [], + moment=_vector(loads.moment) or [], + moment_magnitude=loads.moment_magnitude, + center_of_pressure=_vector(loads.center_of_pressure), + center_of_pressure_status=loads.center_of_pressure_status, + residual_couple=loads.residual_couple, + pressure_weighted_centroid=_vector( + loads.pressure_weighted_centroid), + max_pressure=loads.max_pressure, + max_shear_stress=loads.max_shear_stress, + max_heat_flux=loads.max_heat_flux, + total_heat_load=loads.total_heat_load, + affected_area=loads.affected_area, + struck_faces=loads.n_struck_faces, + coefficients=dict(loads.coefficients), + coefficients_available=loads.has_coefficients, + **metadata, + ) + + # --------------------------------------------------------- serialization + def to_dict(self) -> dict[str, Any]: + """Nested dictionary form (JSON metadata document).""" + return asdict(self) + + def to_row(self) -> dict[str, Any]: + """Flat dictionary form (CSV row). + + Vectors are expanded to ``_x/_y/_z`` columns and coefficients to + one column each, so the CSV is directly plottable. + """ + row: dict[str, Any] = { + "study_name": self.study_name, + "case_id": self.case_id, + "component": self.component, + "firing_id": self.firing_id, + "geometry_id": self.geometry_id, + "mesh_faces": self.mesh_faces, + "component_faces": self.component_faces, + "component_area": self.component_area, + "coordinate_system": self.coordinate_system, + "plate_angle_deg": self.plate_angle_deg, + "source_distance": self.source_distance, + "firing_duration_s": self.firing_duration_s, + "thrusters": " ".join(str(t) for t in self.thrusters), + "plume_model": self.plume_model, + } + _expand(row, "plume_source_position", self.plume_source_position) + _expand(row, "target_normal", self.target_normal) + _expand(row, "target_tangent", self.target_tangent) + _expand(row, "moment_reference_point", self.moment_reference_point) + _expand(row, "pressure_force", self.pressure_force) + _expand(row, "shear_force", self.shear_force) + _expand(row, "force", self.force) + row["force_magnitude"] = self.force_magnitude + _expand(row, "pressure_moment", self.pressure_moment) + _expand(row, "shear_moment", self.shear_moment) + _expand(row, "moment", self.moment) + row["moment_magnitude"] = self.moment_magnitude + _expand(row, "center_of_pressure", self.center_of_pressure) + row["center_of_pressure_status"] = self.center_of_pressure_status + row["residual_couple"] = self.residual_couple + _expand(row, "pressure_weighted_centroid", + self.pressure_weighted_centroid) + row.update({ + "max_pressure": self.max_pressure, + "max_shear_stress": self.max_shear_stress, + "max_heat_flux": self.max_heat_flux, + "total_heat_load": self.total_heat_load, + "affected_area": self.affected_area, + "struck_faces": self.struck_faces, + "coefficients_available": self.coefficients_available, + }) + for name, value in sorted(self.coefficients.items()): + row[f"coeff_{name}"] = value + row.update({ + "vtk_path": self.vtk_path or "", + "jfh_path": self.jfh_path or "", + "config_path": self.config_path, + "case_dir": self.case_dir, + "code_version": self.code_version, + "generated_at": self.generated_at, + }) + return row + + # ------------------------------------------------------------ accessors + def quantity(self, name: str) -> float | list[float] | None: + """Look up a comparable quantity by name (see reference_data). + + Supports the integrated loads, peaks, areas and coefficients under + stable names, returning None when the quantity is unavailable for + this result (e.g. a coefficient with no normalization inputs). + """ + direct: dict[str, float | list[float] | None] = { + "force": self.force, + "force_magnitude": self.force_magnitude, + "pressure_force": self.pressure_force, + "shear_force": self.shear_force, + "moment": self.moment, + "moment_magnitude": self.moment_magnitude, + "center_of_pressure": self.center_of_pressure, + "max_pressure": self.max_pressure, + "max_shear_stress": self.max_shear_stress, + "max_heat_flux": self.max_heat_flux, + "total_heat_load": self.total_heat_load, + "affected_area": self.affected_area, + } + if name in direct: + return direct[name] + if name in self.coefficients: + return self.coefficients[name] + return None + + +def _expand(row: dict[str, Any], name: str, + vector: Sequence[float] | None) -> None: + for axis, index in (("x", 0), ("y", 1), ("z", 2)): + row[f"{name}_{axis}"] = (float(vector[index]) + if vector is not None else "") + + +@dataclass +class StudyResults: + """All case results of one study run, plus study-level provenance.""" + + study_name: str + output_dir: str + provenance: dict[str, Any] + cases: list[CaseResult] = field(default_factory=list) + summary_csv_path: str | None = None + metadata_path: str | None = None + plot_paths: list[str] = field(default_factory=list) + + def __len__(self) -> int: + return len(self.cases) + + def __iter__(self) -> Iterable[CaseResult]: # type: ignore[override] + return iter(self.cases) + + def rows(self) -> list[dict[str, Any]]: + return [case.to_row() for case in self.cases] + + def for_component(self, component: str) -> list[CaseResult]: + return [case for case in self.cases if case.component == component] + + # --------------------------------------------------------------- output + def write_csv(self, path: str | os.PathLike[str]) -> str: + """Write the flat per-case summary CSV; returns the path written.""" + rows = self.rows() + if not rows: + raise ValueError("no case results to write") + columns: list[str] = [] + for row in rows: + for key in row: + if key not in columns: + columns.append(key) + path = os.fspath(path) + os.makedirs(os.path.dirname(os.path.abspath(path)), exist_ok=True) + with open(path, "w", encoding="utf-8", newline="") as handle: + writer = csv.DictWriter(handle, fieldnames=columns) + writer.writeheader() + for row in rows: + writer.writerow(row) + self.summary_csv_path = path + return path + + def write_metadata(self, path: str | os.PathLike[str]) -> str: + """Write the structured JSON metadata document; returns the path.""" + path = os.fspath(path) + os.makedirs(os.path.dirname(os.path.abspath(path)), exist_ok=True) + document = { + "schema": "pyrpod.plume_validation_study/1", + "study_name": self.study_name, + "generated_at": datetime.now(timezone.utc).isoformat( + timespec="seconds"), + "provenance": self.provenance, + "n_cases": len(self.cases), + "cases": [case.to_dict() for case in self.cases], + } + with open(path, "w", encoding="utf-8", newline="\n") as handle: + json.dump(document, handle, indent=2, sort_keys=False) + handle.write("\n") + self.metadata_path = path + return path From 3d6583877ba88c893c46ab9f6900c80094a29912 Mon Sep 17 00:00:00 2001 From: andytorrestb Date: Mon, 27 Jul 2026 22:53:58 -0500 Subject: [PATCH 06/14] mdao: add generic external-reference comparison Compares study results against independently generated data WITHOUT knowing its origin: a reference record is named quantities attached to matching keys, loadable from CSV, JSON or YAML, and whether the numbers came from a DSMC solver, an analytical solution, an experiment or another code changes nothing. No producer-specific importer exists here. Metrics are applied only where they are mathematically meaningful: absolute error, relative error (None when the reference is zero, never infinity), normalized RMSE, peak-value error, integrated-load error (compared as vectors, so a load of the right magnitude pointing the wrong way is an error) and center-of-pressure displacement. A quantity the reference supplies but the result does not is reported as missing, and a case with no matching record is listed as unmatched -- nothing is fabricated or defaulted. The CSV layout is exactly what StudyResults.write_csv emits, so an external producer can be transformed into it column-for-column. Co-Authored-By: Claude Opus 5 --- pyrpod/mdao/reference_data.py | 467 ++++++++++++++++++++++++++++++++++ 1 file changed, 467 insertions(+) create mode 100644 pyrpod/mdao/reference_data.py diff --git a/pyrpod/mdao/reference_data.py b/pyrpod/mdao/reference_data.py new file mode 100644 index 0000000..455347c --- /dev/null +++ b/pyrpod/mdao/reference_data.py @@ -0,0 +1,467 @@ +""" +Generic external reference-data comparison. + +This module compares PyRPOD study results against INDEPENDENTLY GENERATED +data without knowing where that data came from. A reference record is just a +set of named quantities attached to matching keys (component, plate angle, +source distance, or an explicit case id); whether those numbers came from a +DSMC solver, an analytical solution, a wind-tunnel test or another numerical +method is irrelevant to everything here, and no importer for any particular +producer exists in this module. + +Supported reference formats -- all already used by the repository, no new +dependency: + +* **CSV** -- one row per record. Key columns (``case_id``, ``component``, + ``plate_angle_deg``, ``source_distance``) are recognized by name; every + other numeric column becomes a quantity. Columns named ``_x``, + ``_y`` and ``_z`` are assembled into the vector quantity + ````. This is exactly the layout + :meth:`pyrpod.mdao.study_results.StudyResults.write_csv` emits, so an + external producer can be transformed into it column-for-column. +* **JSON / YAML** -- ``{'label': ..., 'source': ..., 'units': {...}, + 'records': [{'key': {...}, 'quantities': {...}}, ...]}``. + +Metrics +------- +Only mathematically applicable metrics are computed for a given quantity: +absolute error and relative error for scalars and vectors, normalized RMSE +for sampled arrays, peak-value error for arrays, integrated-load error for +the force/moment resultants, and center-of-pressure displacement for the +center of pressure. A quantity absent from the reference is reported as +``missing_reference`` -- never defaulted, never fabricated. +""" + +from __future__ import annotations + +import csv +import json +import math +import os +from dataclasses import dataclass, field +from typing import Any, Iterable, Mapping, Sequence + +import numpy as np +from numpy.typing import NDArray + +from pyrpod.mdao.study_results import CaseResult + +__all__ = [ + "ComparisonReport", + "QuantityComparison", + "ReferenceDataset", + "ReferenceRecord", + "absolute_error", + "center_of_pressure_displacement", + "compare_case", + "compare_results", + "integrated_load_error", + "load_reference_dataset", + "normalized_rmse", + "peak_error", + "relative_error", +] + +#: Key columns recognized when a reference record is matched to a case. +KEY_FIELDS = ("case_id", "component", "plate_angle_deg", "source_distance") + +#: Tolerance used when matching numeric keys (angles in degrees, distances in +#: the case's length unit). +KEY_TOLERANCE = 1e-6 + + +# --------------------------------------------------------------------- metrics +def _as_array(value: Any) -> NDArray[np.float64]: + return np.asarray(value, dtype=float).reshape(-1) + + +def absolute_error(candidate: Any, reference: Any) -> float: + """|candidate - reference|; the Euclidean norm for vector quantities.""" + diff = _as_array(candidate) - _as_array(reference) + return float(np.linalg.norm(diff)) + + +def relative_error(candidate: Any, reference: Any) -> float | None: + """Absolute error normalized by the reference magnitude. + + Returns None when the reference magnitude is zero, where a relative + error is undefined rather than infinite. + """ + scale = float(np.linalg.norm(_as_array(reference))) + if scale == 0.0: + return None + return absolute_error(candidate, reference) / scale + + +def normalized_rmse(candidate: Any, reference: Any, + norm: str = "range") -> float | None: + """RMSE of a sampled array normalized by the reference range or mean. + + ``norm='range'`` (default) divides by ``max(ref) - min(ref)``; ``'mean'`` + divides by ``|mean(ref)|``. Returns None when the chosen normalizer is + zero, or when the arrays have different lengths (nothing is resampled or + interpolated here). + """ + values = _as_array(candidate) + target = _as_array(reference) + if values.size != target.size or values.size == 0: + return None + rmse = float(np.sqrt(np.mean((values - target) ** 2))) + if norm == "mean": + scale = abs(float(np.mean(target))) + else: + scale = float(np.max(target) - np.min(target)) + if scale == 0.0: + return None + return rmse / scale + + +def peak_error(candidate: Any, reference: Any) -> tuple[float, float | None]: + """Absolute and relative error of the peak (maximum) sampled value.""" + peak_candidate = float(np.max(_as_array(candidate))) + peak_reference = float(np.max(_as_array(reference))) + absolute = abs(peak_candidate - peak_reference) + if peak_reference == 0.0: + return absolute, None + return absolute, absolute / abs(peak_reference) + + +def integrated_load_error(candidate: Any, reference: Any + ) -> tuple[float, float | None]: + """Absolute and relative error of an integrated load (force or moment). + + Vector loads are compared as vectors (the error is the norm of the + difference, not the difference of the norms), so a load of the right + magnitude pointing the wrong way is reported as an error. + """ + absolute = absolute_error(candidate, reference) + return absolute, relative_error(candidate, reference) + + +def center_of_pressure_displacement(candidate: Any, reference: Any, + reference_length: float | None = None + ) -> tuple[float, float | None]: + """Distance between two centers of pressure, optionally normalized.""" + distance = absolute_error(candidate, reference) + if reference_length in (None, 0.0): + return distance, None + return distance, distance / float(reference_length) # type: ignore[arg-type] + + +# ------------------------------------------------------------------ dataset +@dataclass(frozen=True) +class ReferenceRecord: + """One reference data point: matching keys plus named quantities. + + ``quantities`` values may be scalars, vectors (lists of three numbers) or + sampled arrays; the comparison layer selects metrics accordingly. + """ + + key: dict[str, Any] + quantities: dict[str, Any] + + def matches(self, result: CaseResult) -> bool: + """Whether this record describes the given case result.""" + for name, value in self.key.items(): + if value is None: + continue + actual = getattr(result, name, None) + if actual is None: + return False + if isinstance(value, (int, float)) and not isinstance(value, bool): + if not math.isclose(float(actual), float(value), + abs_tol=KEY_TOLERANCE, rel_tol=1e-9): + return False + elif str(actual) != str(value): + return False + return True + + +@dataclass +class ReferenceDataset: + """A set of reference records with provenance, from any producer. + + Attributes + ---------- + label : str + Human-readable name of the dataset, reported with every comparison + (e.g. ``'Cai 2016 exact'``, ``'DSMC run 12'``). The comparison code + never behaves differently based on this string. + source : str + Where the data came from (path, DOI, run id). + units : dict + Unit labels declared by the producer, recorded for traceability. + records : list of ReferenceRecord + """ + + label: str = "reference" + source: str = "" + units: dict[str, str] = field(default_factory=dict) + records: list[ReferenceRecord] = field(default_factory=list) + + def __len__(self) -> int: + return len(self.records) + + def match(self, result: CaseResult) -> ReferenceRecord | None: + """First record matching the given case result, or None.""" + for record in self.records: + if record.matches(result): + return record + return None + + def quantity_names(self) -> list[str]: + names: list[str] = [] + for record in self.records: + for name in record.quantities: + if name not in names: + names.append(name) + return names + + +def load_reference_dataset(path: str | os.PathLike[str], + label: str | None = None) -> ReferenceDataset: + """Load a reference dataset from CSV, JSON or YAML (by file extension).""" + path = os.fspath(path) + if not os.path.isfile(path): + raise FileNotFoundError(f"reference data not found: {path!r}") + suffix = os.path.splitext(path)[1].lower() + if suffix == ".csv": + dataset = _load_csv(path) + elif suffix in (".json",): + dataset = _load_structured(json.loads( + _read_text(path)), path) + elif suffix in (".yaml", ".yml"): + import yaml + + dataset = _load_structured(yaml.safe_load(_read_text(path)), path) + else: + raise ValueError( + f"unsupported reference-data format {suffix!r}; use .csv, .json " + "or .yaml") + if label: + dataset.label = label + return dataset + + +def _read_text(path: str) -> str: + with open(path, "r", encoding="utf-8") as handle: + return handle.read() + + +def _load_csv(path: str) -> ReferenceDataset: + records: list[ReferenceRecord] = [] + with open(path, "r", encoding="utf-8", newline="") as handle: + for raw in csv.DictReader(handle): + key: dict[str, Any] = {} + scalars: dict[str, float] = {} + for column, text in raw.items(): + if column is None or text is None or text == "": + continue + if column in KEY_FIELDS: + key[column] = _maybe_number(text) + continue + number = _maybe_number(text) + if isinstance(number, float): + scalars[column] = number + records.append(ReferenceRecord(key=key, + quantities=_assemble_vectors(scalars))) + return ReferenceDataset(label=os.path.basename(path), source=path, + records=records) + + +def _load_structured(data: Any, path: str) -> ReferenceDataset: + if not isinstance(data, Mapping): + raise ValueError( + f"reference data {path!r} must be a mapping with a 'records' list") + raw_records = data.get("records") + if not isinstance(raw_records, Sequence): + raise ValueError(f"reference data {path!r} has no 'records' list") + records = [] + for index, entry in enumerate(raw_records): + if not isinstance(entry, Mapping): + raise ValueError( + f"reference data {path!r}: record {index} is not a mapping") + records.append(ReferenceRecord( + key=dict(entry.get("key") or {}), + quantities=dict(entry.get("quantities") or {}))) + return ReferenceDataset( + label=str(data.get("label", os.path.basename(path))), + source=str(data.get("source", path)), + units={str(k): str(v) for k, v in (data.get("units") or {}).items()}, + records=records) + + +def _maybe_number(text: str) -> Any: + try: + return float(text) + except (TypeError, ValueError): + return text + + +def _assemble_vectors(scalars: Mapping[str, float]) -> dict[str, Any]: + """Fold ``_x/_y/_z`` column triplets into vector quantities.""" + quantities: dict[str, Any] = {} + consumed: set[str] = set() + for name in scalars: + if not name.endswith("_x"): + continue + stem = name[:-2] + triplet = [f"{stem}_x", f"{stem}_y", f"{stem}_z"] + if all(component in scalars for component in triplet): + quantities[stem] = [scalars[component] for component in triplet] + consumed.update(triplet) + for name, value in scalars.items(): + if name not in consumed: + quantities[name] = value + return quantities + + +# --------------------------------------------------------------- comparison +@dataclass(frozen=True) +class QuantityComparison: + """Comparison of one quantity for one case against one reference record.""" + + case_id: str + component: str + quantity: str + status: str + candidate: Any = None + reference: Any = None + absolute_error: float | None = None + relative_error: float | None = None + normalized_rmse: float | None = None + peak_absolute_error: float | None = None + peak_relative_error: float | None = None + displacement: float | None = None + + def to_row(self) -> dict[str, Any]: + return { + "case_id": self.case_id, + "component": self.component, + "quantity": self.quantity, + "status": self.status, + "candidate": _scalar_or_text(self.candidate), + "reference": _scalar_or_text(self.reference), + "absolute_error": self.absolute_error, + "relative_error": self.relative_error, + "normalized_rmse": self.normalized_rmse, + "peak_absolute_error": self.peak_absolute_error, + "peak_relative_error": self.peak_relative_error, + "displacement": self.displacement, + } + + +def _scalar_or_text(value: Any) -> Any: + if value is None: + return "" + array = np.asarray(value, dtype=float).reshape(-1) + if array.size == 1: + return float(array[0]) + return " ".join(f"{v:.9g}" for v in array) + + +def compare_case(result: CaseResult, record: ReferenceRecord, + reference_length: float | None = None, + ) -> list[QuantityComparison]: + """Compare every quantity the reference record supplies for one case.""" + comparisons: list[QuantityComparison] = [] + for name, reference_value in record.quantities.items(): + candidate = result.quantity(name) + if candidate is None: + comparisons.append(QuantityComparison( + case_id=result.case_id, component=result.component, + quantity=name, status="missing_candidate", + reference=reference_value)) + continue + comparisons.append(_compare_quantity( + result, name, candidate, reference_value, reference_length)) + return comparisons + + +def _compare_quantity(result: CaseResult, name: str, candidate: Any, + reference: Any, reference_length: float | None, + ) -> QuantityComparison: + candidate_array = _as_array(candidate) + reference_array = _as_array(reference) + if candidate_array.size != reference_array.size: + return QuantityComparison( + case_id=result.case_id, component=result.component, quantity=name, + status="shape_mismatch", candidate=candidate, reference=reference) + + absolute = absolute_error(candidate_array, reference_array) + relative = relative_error(candidate_array, reference_array) + nrmse = (normalized_rmse(candidate_array, reference_array) + if candidate_array.size > 3 else None) + peak_absolute: float | None = None + peak_relative: float | None = None + if candidate_array.size > 1: + peak_absolute, peak_relative = peak_error(candidate_array, + reference_array) + displacement: float | None = None + if name == "center_of_pressure" and candidate_array.size == 3: + displacement, _normalized = center_of_pressure_displacement( + candidate_array, reference_array, reference_length) + + return QuantityComparison( + case_id=result.case_id, component=result.component, quantity=name, + status="compared", candidate=candidate, reference=reference, + absolute_error=absolute, relative_error=relative, + normalized_rmse=nrmse, peak_absolute_error=peak_absolute, + peak_relative_error=peak_relative, displacement=displacement) + + +@dataclass +class ComparisonReport: + """All quantity comparisons for a study run against one dataset.""" + + label: str + source: str + comparisons: list[QuantityComparison] = field(default_factory=list) + unmatched_cases: list[str] = field(default_factory=list) + + def __len__(self) -> int: + return len(self.comparisons) + + def rows(self) -> list[dict[str, Any]]: + return [comparison.to_row() for comparison in self.comparisons] + + def for_quantity(self, name: str) -> list[QuantityComparison]: + return [c for c in self.comparisons if c.quantity == name] + + def max_relative_error(self, name: str) -> float | None: + """Largest relative error recorded for a quantity, or None.""" + values = [c.relative_error for c in self.for_quantity(name) + if c.relative_error is not None] + return max(values) if values else None + + def write_csv(self, path: str | os.PathLike[str]) -> str: + rows = self.rows() + if not rows: + raise ValueError("no comparisons to write") + path = os.fspath(path) + os.makedirs(os.path.dirname(os.path.abspath(path)), exist_ok=True) + with open(path, "w", encoding="utf-8", newline="") as handle: + writer = csv.DictWriter(handle, fieldnames=list(rows[0].keys())) + writer.writeheader() + writer.writerows(rows) + return path + + +def compare_results(results: Iterable[CaseResult], dataset: ReferenceDataset, + reference_length: float | None = None) -> ComparisonReport: + """Compare a study's case results against a reference dataset. + + Cases with no matching reference record are listed in + ``ComparisonReport.unmatched_cases`` rather than being silently dropped + or compared against invented values. + """ + report = ComparisonReport(label=dataset.label, source=dataset.source) + for result in results: + record = dataset.match(result) + if record is None: + report.unmatched_cases.append(f"{result.case_id}/{result.component}") + continue + report.comparisons.extend( + compare_case(result, record, + reference_length=reference_length)) + return report From 00ac6a100cb9456061c92aad6ca4f3093abd0219 Mon Sep 17 00:00:00 2001 From: andytorrestb Date: Mon, 27 Jul 2026 22:54:23 -0500 Subject: [PATCH 07/14] mdao: add optional sweep trend plots Force, moment and heat flux versus plate angle and source distance, center-of-pressure travel, and a study-versus-reference error summary. Every figure is optional output: the module is imported lazily by the study engines and pins matplotlib's non-interactive Agg backend, so no run and no automated test needs graphical output. Series are keyed on the swept pose rather than on the case, so both study engines plot identically. Co-Authored-By: Claude Opus 5 --- pyrpod/mdao/study_plots.py | 231 +++++++++++++++++++++++++++++++++++++ 1 file changed, 231 insertions(+) create mode 100644 pyrpod/mdao/study_plots.py diff --git a/pyrpod/mdao/study_plots.py b/pyrpod/mdao/study_plots.py new file mode 100644 index 0000000..05a2671 --- /dev/null +++ b/pyrpod/mdao/study_plots.py @@ -0,0 +1,231 @@ +""" +Optional parameter-sweep trend plots for plume validation studies. + +Every function here is OPTIONAL output: nothing in the study engine, and no +automated test, requires a figure to be produced. The module is imported +lazily by :meth:`pyrpod.mdao.plume_validation.PlumeValidationStudy.plot` and +selects matplotlib's non-interactive Agg backend on import, so a headless or +CI run never opens a window (the same artifact-control behavior the RPOD +integration and verification tests use). + +Produced figures (written to the study's ``plots`` subdirectory): + +* ``force_vs_angle.png`` -- resultant force magnitude vs plate angle +* ``moment_vs_angle.png`` -- resultant moment magnitude vs plate angle +* ``heat_flux_vs_angle.png`` -- peak surface heat flux vs plate angle +* ``force_vs_distance.png`` -- resultant force magnitude vs source distance +* ``moment_vs_distance.png`` -- resultant moment magnitude vs source distance +* ``center_of_pressure.png`` -- center-of-pressure travel across the sweep +* ``reference_comparison.png`` -- study vs external reference, when a + comparison report is supplied +""" + +from __future__ import annotations + +import os +from typing import Iterable, Sequence + +import matplotlib +import numpy as np + +matplotlib.use("Agg") + +import matplotlib.pyplot as plt # noqa: E402 (backend must be set first) + +from pyrpod.mdao.reference_data import ComparisonReport # noqa: E402 +from pyrpod.mdao.study_results import CaseResult, StudyResults # noqa: E402 + +__all__ = ["plot_sweep_trends"] + +#: Categorical palette used across the sweep plots (fixed slot order, the +#: same CVD-safe palette the existing sweep study script uses). +PALETTE = ["#2a78d6", "#008300", "#e87ba4", "#eda100", "#1baf7a"] + + +def _distances(cases: Sequence[CaseResult]) -> list[float]: + return sorted({case.source_distance for case in cases}) + + +def _angles(cases: Sequence[CaseResult]) -> list[float]: + return sorted({case.plate_angle_deg for case in cases}) + + +def _series(cases: Sequence[CaseResult], fixed: str, value: float, + variable: str, quantity: str) -> tuple[list[float], list[float]]: + selected = [case for case in cases + if getattr(case, fixed) == value] + selected.sort(key=lambda case: getattr(case, variable)) + xs = [float(getattr(case, variable)) for case in selected] + ys = [float(getattr(case, quantity)) for case in selected] + return xs, ys + + +def _line_figure(path: str, title: str, xlabel: str, ylabel: str, + series: Iterable[tuple[str, list[float], list[float]]], + log_y: bool = False) -> str: + figure, axes = plt.subplots(figsize=(7, 4.5)) + for color, (label, xs, ys) in zip(PALETTE * 8, series): + axes.plot(xs, ys, "-o", color=color, lw=2, ms=4, label=label) + axes.set_xlabel(xlabel) + axes.set_ylabel(ylabel) + axes.set_title(title, fontsize=10) + if log_y: + axes.set_yscale("log") + axes.grid(True, color="0.9", lw=0.8) + axes.legend(fontsize=8) + figure.savefig(path, dpi=200, bbox_inches="tight") + plt.close(figure) + return path + + +def plot_sweep_trends(results: StudyResults, out_dir: str, + comparison: ComparisonReport | None = None, + component: str | None = None) -> list[str]: + """Write the sweep trend figures; returns the paths written. + + Parameters + ---------- + results : StudyResults + Executed study results. + out_dir : str + Directory to write the figures to (created when missing). + comparison : ComparisonReport, optional + When supplied, an extra study-vs-reference figure is written. + component : str, optional + Restrict the plots to one target component. Defaults to the first + component present in the results. + """ + cases = list(results.cases) + if not cases: + return [] + component = component or cases[0].component + cases = [case for case in cases if case.component == component] + # One point per swept pose: keep the first record of each (angle, + # distance). This is the first firing of each case when every pose has + # its own history, and the first firing of each pose's slice when one + # history spans the sweep -- both engines plot identically. + seen: set[tuple[float, float]] = set() + unique: list[CaseResult] = [] + for case in cases: + pose = (case.plate_angle_deg, case.source_distance) + if pose in seen: + continue + seen.add(pose) + unique.append(case) + cases = unique + if not cases: + return [] + + os.makedirs(out_dir, exist_ok=True) + written: list[str] = [] + + distances = _distances(cases) + angles = _angles(cases) + + for quantity, filename, ylabel, title in ( + ("force_magnitude", "force_vs_angle.png", "|F| (N)", + "Resultant force vs plate angle"), + ("moment_magnitude", "moment_vs_angle.png", "|M| (N*m)", + "Resultant moment vs plate angle"), + ("max_heat_flux", "heat_flux_vs_angle.png", + "peak heat flux (W/m^2)", "Peak heat flux vs plate angle")): + series = [(f"d = {distance:g}", + *_series(cases, "source_distance", distance, + "plate_angle_deg", quantity)) + for distance in distances] + written.append(_line_figure( + os.path.join(out_dir, filename), f"{title} ({component})", + "plate angle (deg), 0 = head-on", ylabel, series)) + + if len(distances) > 1: + for quantity, filename, ylabel, title in ( + ("force_magnitude", "force_vs_distance.png", "|F| (N)", + "Resultant force vs source distance"), + ("moment_magnitude", "moment_vs_distance.png", "|M| (N*m)", + "Resultant moment vs source distance")): + series = [(f"alpha = {angle:g} deg", + *_series(cases, "plate_angle_deg", angle, + "source_distance", quantity)) + for angle in angles] + written.append(_line_figure( + os.path.join(out_dir, filename), f"{title} ({component})", + "source distance", ylabel, series)) + + cop_path = _plot_center_of_pressure(cases, out_dir, component) + if cop_path: + written.append(cop_path) + + if comparison is not None and len(comparison): + written.append(_plot_reference_comparison(comparison, out_dir)) + + return written + + +def _plot_center_of_pressure(cases: Sequence[CaseResult], out_dir: str, + component: str) -> str | None: + """Center-of-pressure travel, projected on the two axes that vary most. + + Cases whose center of pressure is unavailable (zero load, or a resultant + that is cancellation noise) simply do not appear. + """ + located = [(case, np.asarray(case.center_of_pressure, dtype=float)) + for case in cases if case.center_of_pressure is not None] + if not located: + return None + + points = np.vstack([point for _, point in located]) + spread = points.max(axis=0) - points.min(axis=0) + order = sorted(range(3), key=lambda axis: float(spread[axis]), + reverse=True)[:2] + first, second = sorted(order) + labels = "xyz" + + distances = sorted({case.source_distance for case, _ in located}) + figure, axes = plt.subplots(figsize=(6, 5)) + for color, distance in zip(PALETTE * 8, distances): + subset = sorted((entry for entry in located + if entry[0].source_distance == distance), + key=lambda entry: entry[0].plate_angle_deg) + axes.plot([point[first] for _, point in subset], + [point[second] for _, point in subset], + "-o", color=color, lw=1.6, ms=4, label=f"d = {distance:g}") + axes.set_xlabel(f"CoP {labels[first]} (m)") + axes.set_ylabel(f"CoP {labels[second]} (m)") + axes.set_title(f"Center-of-pressure travel across the sweep ({component})", + fontsize=10) + axes.grid(True, color="0.9", lw=0.8) + axes.legend(fontsize=8) + path = os.path.join(out_dir, "center_of_pressure.png") + figure.savefig(path, dpi=200, bbox_inches="tight") + plt.close(figure) + return path + + +def _plot_reference_comparison(comparison: ComparisonReport, + out_dir: str) -> str: + """Bar chart of the relative error per compared quantity.""" + quantities: dict[str, list[float]] = {} + for entry in comparison.comparisons: + if entry.relative_error is None: + continue + quantities.setdefault(entry.quantity, []).append(entry.relative_error) + + figure, axes = plt.subplots(figsize=(7, 4.5)) + names = sorted(quantities) + maxima = [max(quantities[name]) for name in names] + means = [sum(quantities[name]) / len(quantities[name]) for name in names] + positions = range(len(names)) + axes.bar([p - 0.2 for p in positions], maxima, width=0.4, + color=PALETTE[0], label="max") + axes.bar([p + 0.2 for p in positions], means, width=0.4, + color=PALETTE[1], label="mean") + axes.set_xticks(list(positions)) + axes.set_xticklabels(names, rotation=30, ha="right", fontsize=8) + axes.set_ylabel("relative error") + axes.set_title(f"Study vs {comparison.label}", fontsize=10) + axes.grid(True, axis="y", color="0.9", lw=0.8) + axes.legend(fontsize=8) + path = os.path.join(out_dir, "reference_comparison.png") + figure.savefig(path, dpi=200, bbox_inches="tight") + plt.close(figure) + return path From 041eb577ebd918faf2ac0cb2cdaf937d567f7fd2 Mon Sep 17 00:00:00 2001 From: andytorrestb Date: Mon, 27 Jul 2026 22:54:22 -0500 Subject: [PATCH 08/14] mdao: add the shared study runtime The plumbing every study engine needs, in one place: build the existing PyRPOD case objects (failing early on an unsupported plume model or an unknown thruster id), precompute the target geometry and component face selection once, read a generated JFH back through the normal JetFiringHistory parser, run one history through the plume-strike calculation, and turn a firing's per-face arrays into per-component result records. Strike execution has two paths: with VTK output enabled the full PlumeStrikeEstimationStudy.jfh_plume_strikes() pipeline runs with the target vehicle's output root redirected, so studies that execute several histories cannot overwrite one another's artifacts; with it disabled the pipeline's own per-firing core runs instead -- identical numbers, no files. Living here rather than on a base class lets the two study engines be siblings with no duplicated logic and no inheritance between them. Co-Authored-By: Claude Opus 5 --- pyrpod/mdao/study_runtime.py | 327 +++++++++++++++++++++++++++++++++++ 1 file changed, 327 insertions(+) create mode 100644 pyrpod/mdao/study_runtime.py diff --git a/pyrpod/mdao/study_runtime.py b/pyrpod/mdao/study_runtime.py new file mode 100644 index 0000000..0841154 --- /dev/null +++ b/pyrpod/mdao/study_runtime.py @@ -0,0 +1,327 @@ +""" +Shared runtime for prescribed plume/target validation studies. + +Both study styles -- :class:`pyrpod.mdao.plume_validation.PlumeValidationStudy` +(one Jet Firing History per angle-distance case) and +:class:`pyrpod.mdao.parameter_sweep.ParameterSweepStudy` (one history spanning +the whole sweep) -- need exactly the same plumbing: build the case objects, +precompute the target geometry, read back a generated JFH, run the strike +calculation, and turn one firing's per-face arrays into result records. That +plumbing lives here so the two classes stay siblings with no duplicated logic +and no inheritance between them. + +Nothing in this module decides HOW a sweep is decomposed into histories; that +is the studies' own business. +""" + +from __future__ import annotations + +import logging +import os +from dataclasses import dataclass +from datetime import datetime, timezone +from typing import Any + +import numpy as np +from numpy.typing import NDArray + +from pyrpod.mdao.firing_plan import Firing +from pyrpod.mdao.study_config import StudyConfig +from pyrpod.mdao.study_results import CaseResult, code_version +from pyrpod.mdao.surface_loads import ( + face_areas, + flow_directions, + integrate_component_loads, + select_component_faces, +) +from pyrpod.mission import MissionEnvironment +from pyrpod.plume.PlumeStrikeCalculator import ( + compute_face_centroids, + compute_plume_strikes, +) +from pyrpod.rpod import JetFiringHistory, PlumeStrikeEstimationStudy +from pyrpod.util.io.fs import ensure_dir +from pyrpod.vehicle import TargetVehicle, VisitingVehicle + +logger = logging.getLogger(__name__) + +__all__ = [ + "CaseAssets", + "TargetGeometry", + "build_case_results", + "component_envelope", + "compute_strikes", + "load_case_assets", + "read_generated_jfh", + "study_provenance", + "utc_timestamp", +] + +#: Cumulative per-face fields the strike pipeline accumulates across a single +#: Jet Firing History. They are a sweep ENVELOPE only when one history spans +#: the sweep; with one history per case they simply restate that case. +CUMULATIVE_FIELDS = ("cum_strikes", "max_pressures", "max_shears", + "cum_heat_flux_load") + + +@dataclass +class CaseAssets: + """The existing PyRPOD objects a study runs on.""" + + target_vehicle: Any + visiting_vehicle: Any + environment: Any + + +@dataclass +class TargetGeometry: + """Per-face target geometry, computed once and reused by every firing.""" + + mesh: Any + centroids: NDArray[np.float64] + normals: NDArray[np.float64] + areas: NDArray[np.float64] + n_faces: int + components: list[tuple[str, NDArray[np.int64]]] + + @classmethod + def from_config(cls, config: StudyConfig, mesh: Any) -> "TargetGeometry": + centroids = compute_face_centroids(mesh.vectors) + n_faces = int(len(mesh.vectors)) + components = [ + (component.name, + select_component_faces(component, centroids, n_faces)) + for component in config.target.components + ] + return cls(mesh=mesh, centroids=centroids, + normals=mesh.get_unit_normals(), + areas=face_areas(mesh.vectors), n_faces=n_faces, + components=components) + + +def utc_timestamp() -> str: + return datetime.now(timezone.utc).isoformat(timespec="seconds") + + +def load_case_assets(config: StudyConfig) -> CaseAssets: + """Construct the existing PyRPOD case objects for a study. + + Fails early and specifically when the case does not support the study + workflow: an unsupported plume model, or a configured thruster id that + the case's thruster configuration file does not define. + """ + case_dir = config.case_dir + + target_vehicle = TargetVehicle.TargetVehicle(case_dir) + target_vehicle.set_stl() + + visiting_vehicle = VisitingVehicle.VisitingVehicle(case_dir) + visiting_vehicle.set_thruster_config() + visiting_vehicle.set_thruster_metrics() + + thruster_id = config.thruster_id + if thruster_id is not None: + available = list(getattr(visiting_vehicle, "thruster_data", {}) or {}) + if thruster_id not in available: + raise ValueError( + f"thruster.id {thruster_id!r} is not configured in the case's " + f"thruster configuration file; available: {available}") + + environment = MissionEnvironment.MissionEnvironment(case_dir) + configured_model = environment.config["pm"]["kinetics"] + if configured_model != "Simplified": + raise ValueError( + f"case {case_dir!r} configures plume kinetics " + f"{configured_model!r}; this study workflow requires the " + "SimplifiedGasKinetics model ([pm] kinetics = Simplified)") + return CaseAssets(target_vehicle=target_vehicle, + visiting_vehicle=visiting_vehicle, + environment=environment) + + +def study_provenance(config: StudyConfig, geometry: TargetGeometry, + **extra: Any) -> dict[str, Any]: + """Study-level provenance recorded in the metadata document.""" + provenance = config.provenance() + provenance.update({ + "code_version": code_version( + os.path.dirname(os.path.abspath(__file__))), + "mesh_faces": geometry.n_faces, + "geometry_id": config.target.geometry_id, + "sweep_mode": config.sweep.mode, + "n_cases": config.n_cases, + "n_firings_per_pose": config.sweep.n_firings, + "total_firings": config.sweep.total_firings, + "components": [name for name, _ in geometry.components], + "known_limitations": [ + "plume shadowing, occlusion and back-facing geometry are not " + "modeled (existing pipeline face-selection behavior)", + ], + }) + provenance.update(extra) + return provenance + + +def read_generated_jfh(config: StudyConfig, jfh_name: str) -> Any: + """Load a generated JFH through the normal JetFiringHistory reader. + + The JFH object keeps the case's ``config.ini`` (so every other setting is + the case's own) but resolves its JFH asset from the study's own output + tree, which is where generated firing histories live. ``case_dir`` is the + documented asset-resolution root of ``JetFiringHistory``, so pointing it + at the study output directory is exactly its intended use. + """ + jfh = JetFiringHistory.JetFiringHistory(config.case_dir) + jfh.case_dir = config.output_dir.rstrip("\\/") + os.sep + if not jfh.config.has_section("jfh"): + jfh.config.add_section("jfh") + jfh.config.set("jfh", "jfh", jfh_name) + jfh.read_jfh() + if getattr(jfh, "JFH", None) is None: + raise ValueError( + f"failed to read generated JFH {jfh_name!r} from " + f"{os.path.join(config.output_dir, 'jfh')}") + return jfh + + +def compute_strikes(config: StudyConfig, assets: CaseAssets, + geometry: TargetGeometry, jfh: Any, + output_root: str) -> tuple[dict[str, dict[str, Any]], + list[str]]: + """Run one Jet Firing History through the plume-strike calculation. + + With VTK output enabled the full + ``PlumeStrikeEstimationStudy.jfh_plume_strikes()`` pipeline runs, writing + the standard per-face strike files to + ``/results/strikes/firing-.vtu``; the target vehicle's + output root is redirected there for the duration of the run so studies + that execute several histories cannot overwrite one another's artifacts. + With VTK disabled the pipeline's own per-firing core runs instead -- + identical numbers, no files. + + Returns + ------- + (dict, list of str) + The per-firing field dictionary keyed ``'1'``..``'N'``, and the + per-firing VTK paths (empty when VTK output is disabled). + """ + if not config.output.write_vtk: + firing_data: dict[str, dict[str, Any]] = {} + for index in range(len(jfh.JFH)): + step = {"thrusters": jfh.JFH[index]["thrusters"], + "xyz": np.array(jfh.JFH[index]["xyz"]), + "dcm": np.array(jfh.JFH[index]["dcm"]), + "t": float(jfh.JFH[index]["t"])} + result = compute_plume_strikes( + geometry.mesh, geometry.normals, assets.visiting_vehicle, + step, assets.environment, face_centroids=geometry.centroids) + firing_data[str(index + 1)] = dict(result) + return firing_data, [] + + study = PlumeStrikeEstimationStudy.PlumeStrikeEstimationStudy( + assets.environment) + study.study_init(jfh, assets.target_vehicle, assets.visiting_vehicle) + + target_vehicle = assets.target_vehicle + original_case_dir = target_vehicle.case_dir + redirected = output_root.rstrip("\\/") + os.sep + ensure_dir(os.path.join(redirected, "results", "strikes")) + try: + target_vehicle.case_dir = redirected + firing_data = study.jfh_plume_strikes() + finally: + target_vehicle.case_dir = original_case_dir + + vtk_paths = [os.path.join(redirected, "results", "strikes", + f"firing-{index}.vtu") + for index in range(len(jfh.JFH))] + return firing_data, vtk_paths + + +def build_case_results(config: StudyConfig, geometry: TargetGeometry, + firing: Firing, per_face: dict[str, Any], *, + case_id: str, firing_id: int, + plate_angle_deg: float, source_distance: float, + jfh_path: str, vtk_path: str | None, + code_version_id: str, timestamp: str, + ) -> list[CaseResult]: + """Integrate one firing's per-face fields into one record per component.""" + flow = flow_directions(geometry.centroids, firing.position) + records: list[CaseResult] = [] + for component_name, face_indices in geometry.components: + loads = integrate_component_loads( + component_name=component_name, + face_indices=face_indices, + centroids=geometry.centroids, unit_normals=geometry.normals, + areas=geometry.areas, + pressures=per_face["pressures"], + shear_stresses=per_face["shear_stress"], + heat_fluxes=per_face["heat_flux_rate"], + strikes=per_face["strikes"], + moment_reference_point=config.loads.moment_reference_point, + flow_unit_vectors=flow, + normalization=config.loads.normalization) + + records.append(CaseResult.from_loads( + loads, + study_name=config.study_name, + case_id=case_id, + firing_id=firing_id, + geometry_id=config.target.geometry_id, + mesh_faces=geometry.n_faces, + coordinate_system=config.coordinate_system, + units=dict(config.units), + plume_source_position=[float(v) for v in firing.position], + plume_source_orientation=[ + float(v) for v in np.asarray(firing.dcm).ravel()], + target_normal=[float(v) for v in config.target.normal], + target_tangent=[float(v) for v in config.target.tangent], + target_reference_point=[ + float(v) for v in config.target.reference_point], + plate_angle_deg=float(plate_angle_deg), + source_distance=float(source_distance), + firing_duration_s=float(firing.duration_s), + thrusters=[int(t) for t in firing.thrusters], + plume_model=config.plume_model, + plume_model_parameters=dict(config.plume_model_parameters), + vtk_path=vtk_path, + jfh_path=jfh_path, + config_path=config.source_path, + case_dir=os.path.abspath(config.case_dir), + code_version=code_version_id, + generated_at=timestamp)) + return records + + +def component_envelope(geometry: TargetGeometry, + cumulative: dict[str, Any], + ) -> dict[str, dict[str, float]]: + """Worst-case per-face statistics over a whole firing history. + + Only meaningful when ONE history spans the sweep: the strike pipeline's + cumulative arrays then hold the maximum pressure and shear any pose + produced on each face, the summed heat-flux load, and the number of times + each face was struck. Reported per component so a multi-component target + keeps them separate. + """ + envelope: dict[str, dict[str, float]] = {} + for component_name, face_indices in geometry.components: + indices = np.asarray(face_indices, dtype=np.int64) + areas = geometry.areas[indices] + strikes = np.asarray(cumulative["cum_strikes"], dtype=float)[indices] + pressures = np.asarray(cumulative["max_pressures"], + dtype=float)[indices] + shears = np.asarray(cumulative["max_shears"], dtype=float)[indices] + heat_load = np.asarray(cumulative["cum_heat_flux_load"], + dtype=float)[indices] + struck = strikes > 0.0 + envelope[component_name] = { + "max_pressure": float(np.max(pressures)), + "max_shear_stress": float(np.max(shears)), + "max_heat_flux_load": float(np.max(heat_load)), + "total_strike_events": float(np.sum(strikes)), + "unique_struck_faces": int(np.count_nonzero(struck)), + "swept_affected_area": float(np.sum(areas[struck])), + "component_area": float(np.sum(areas)), + } + return envelope From fa43ec09e40ff77146d9b7fa0b65b511fc37ad15 Mon Sep 17 00:00:00 2001 From: andytorrestb Date: Mon, 27 Jul 2026 22:54:22 -0500 Subject: [PATCH 09/14] mdao: add the per-case study engine PlumeValidationStudy runs a configured sweep with ONE Jet Firing History per angle-distance combination (sweep.mode: per_case, the default): every pose is an independent case with its own history, strike run and artifacts, so no pose's results can be polluted by another's. That is what a validation matrix compared pose-by-pose against reference data wants. Per case it builds the prescribed firings, writes and re-reads the case's JFH, runs the plume-strike calculation, integrates per-component loads, exports the per-face VTK fields, records structured results and optionally compares against external reference data. Known limitation, documented here and in the emitted provenance: face selection is the pipeline's existing behavior, so plume shadowing, occlusion and self-shadowing of concave targets are not modeled. Co-Authored-By: Claude Opus 5 --- pyrpod/mdao/plume_validation.py | 227 ++++++++++++++++++++++++++++++++ 1 file changed, 227 insertions(+) create mode 100644 pyrpod/mdao/plume_validation.py diff --git a/pyrpod/mdao/plume_validation.py b/pyrpod/mdao/plume_validation.py new file mode 100644 index 0000000..81c6e76 --- /dev/null +++ b/pyrpod/mdao/plume_validation.py @@ -0,0 +1,227 @@ +""" +Per-case engine for prescribed plume/target validation studies. + +:class:`PlumeValidationStudy` runs a configured sweep with ONE Jet Firing +History per angle-distance combination: every pose is an independent case +with its own history, its own strike run and its own artifacts, so no pose's +results can be polluted by another's. It is the ``sweep.mode: per_case`` +(default) engine behind the package-level API + +>>> from pyrpod.mdao.TradeStudy import TradeStudy +>>> study = TradeStudy.from_config('case.yaml') +>>> results = study.run() + +The sibling :class:`pyrpod.mdao.parameter_sweep.ParameterSweepStudy` runs the +same configuration as a SINGLE history spanning every pose +(``sweep.mode: single_jfh``); both share the plumbing in +:mod:`pyrpod.mdao.study_runtime`, so their per-firing numbers are identical +and only the decomposition into histories differs. + +Per case this engine does exactly seven things: build the prescribed firings, +write the case's Jet Firing History, run the plume-strike calculation, +integrate per-component loads, export the per-face VTK fields, record a +structured result, and (optionally) compare against external reference data. + +Existing PyRPOD objects do all the domain work -- ``JetFiringHistory``, +``TargetVehicle``, ``VisitingVehicle``, ``MissionEnvironment`` and +``PlumeStrikeEstimationStudy`` -- so a study inherits the pipeline's +validation, logging, plume physics and VTK conventions rather than +reimplementing them. + +Known limitation +---------------- +Face selection is the pipeline's existing one: a face is struck when it lies +inside the plume wedge/radius and its normal faces the source. Plume +shadowing, occlusion by intervening geometry and self-shadowing of concave +targets are NOT modeled, so a face hidden behind other geometry still +receives load. This branch does not change that behavior; see +``docs/plume_validation_study.md``. +""" + +from __future__ import annotations + +import logging +import os +import time + +from pyrpod.mdao import firing_plan, study_runtime +from pyrpod.mdao.reference_data import ( + ComparisonReport, + ReferenceDataset, + compare_results, + load_reference_dataset, +) +from pyrpod.mdao.study_config import StudyConfig +from pyrpod.mdao.study_results import StudyResults +from pyrpod.mdao.study_runtime import CaseAssets, TargetGeometry +from pyrpod.util.io.fs import ensure_dir + +logger = logging.getLogger(__name__) + +__all__ = ["PlumeValidationStudy", "case_id_for"] + + +def case_id_for(index: int, plate_angle_deg: float, + source_distance: float) -> str: + """Stable, filesystem-safe case identifier.""" + angle = f"{plate_angle_deg:.1f}".replace("-", "m").replace(".", "p") + distance = f"{source_distance:.4g}".replace("-", "m").replace(".", "p") + return f"case{index:03d}_alpha{angle}_d{distance}" + + +class PlumeValidationStudy: + """Prescribed validation sweep with one Jet Firing History per case.""" + + def __init__(self, config: StudyConfig) -> None: + self.config = config + self.results: StudyResults | None = None + self.comparison: ComparisonReport | None = None + + # ------------------------------------------------------------ lifecycle + @classmethod + def from_config(cls, path: str | os.PathLike[str] + ) -> "PlumeValidationStudy": + """Build a study from a YAML configuration file.""" + return cls(StudyConfig.from_yaml(path)) + + # ------------------------------------------------------------------ run + def run(self, write_outputs: bool = True) -> StudyResults: + """Execute every angle x distance case and return structured results. + + Parameters + ---------- + write_outputs : bool, optional + When True (default) the summary CSV, the JSON metadata document + and the optional plots are written to the study output directory. + The per-case VTK export is controlled by the configuration's own + ``output.vtk.enabled`` flag. + + Returns + ------- + StudyResults + One :class:`~pyrpod.mdao.study_results.CaseResult` per case, + component and firing. + """ + config = self.config + started = time.perf_counter() + ensure_dir(config.output_dir) + + assets = study_runtime.load_case_assets(config) + geometry = TargetGeometry.from_config(config, assets.target_vehicle.mesh) + + results = StudyResults( + study_name=config.study_name, output_dir=config.output_dir, + provenance=study_runtime.study_provenance(config, geometry)) + + logger.info("Plume validation study started: study=%s mode=per_case " + "cases=%d firings_per_case=%d components=%d " + "target_faces=%d", config.study_name, config.n_cases, + config.sweep.n_firings, len(geometry.components), + geometry.n_faces) + + for index, (angle, distance) in enumerate(config.sweep.poses): + self._run_case(case_id=case_id_for(index, angle, distance), + pose_index=index, plate_angle_deg=angle, + source_distance=distance, results=results, + assets=assets, geometry=geometry) + + if write_outputs: + self.write_outputs(results) + + self.results = results + if config.reference.path: + self.comparison = self.compare( + load_reference_dataset(config.reference.path, + label=config.reference.label)) + if write_outputs and len(self.comparison): + self.comparison.write_csv(os.path.join( + config.output_dir, "reference_comparison.csv")) + + logger.info("Plume validation study completed: study=%s cases=%d " + "records=%d wall_time_s=%.2f output=%s", + config.study_name, config.n_cases, len(results), + time.perf_counter() - started, config.output_dir) + return results + + # ------------------------------------------------------------ one case + def _run_case(self, *, case_id: str, pose_index: int, + plate_angle_deg: float, source_distance: float, + results: StudyResults, assets: CaseAssets, + geometry: TargetGeometry) -> None: + config = self.config + firings = firing_plan.build_case_firings( + config.sweep, config.target, plate_angle_deg, source_distance, + pose_index=pose_index) + + jfh_path = os.path.join(config.output_dir, "jfh", f"{case_id}.A") + n_written = firing_plan.write_jfh_file(jfh_path, firings) + if n_written != config.sweep.n_firings: + raise ValueError( + f"case {case_id}: wrote {n_written} JFH entries for a " + f"requested n_firings of {config.sweep.n_firings}") + + jfh = study_runtime.read_generated_jfh(config, f"{case_id}.A") + if len(jfh.JFH) != config.sweep.n_firings: + raise ValueError( + f"case {case_id}: JFH file {jfh_path!r} reports " + f"{len(jfh.JFH)} firings, expected {config.sweep.n_firings}") + + firing_data, vtk_paths = study_runtime.compute_strikes( + config, assets, geometry, jfh, + output_root=os.path.join(config.output_dir, "cases", case_id)) + + timestamp = study_runtime.utc_timestamp() + version = str(results.provenance.get("code_version", "unknown")) + + for firing_index, firing in enumerate(firings): + results.cases.extend(study_runtime.build_case_results( + config, geometry, firing, firing_data[str(firing_index + 1)], + case_id=case_id, firing_id=firing_index + 1, + plate_angle_deg=plate_angle_deg, + source_distance=source_distance, jfh_path=jfh_path, + vtk_path=vtk_paths[firing_index] if vtk_paths else None, + code_version_id=version, timestamp=timestamp)) + + # -------------------------------------------------------------- outputs + def write_outputs(self, results: StudyResults) -> None: + """Write the summary CSV, the JSON metadata and any configured plots.""" + config = self.config + results.write_csv(os.path.join(config.output_dir, + config.output.summary_csv)) + results.write_metadata(os.path.join(config.output_dir, + config.output.summary_metadata)) + if config.output.write_plots: + results.plot_paths = self.plot(results) + + # ----------------------------------------------------------- comparison + def compare(self, dataset: ReferenceDataset) -> ComparisonReport: + """Compare this study's results against a reference dataset. + + The dataset is opaque: DSMC, analytical, experimental or another + code, all handled identically (see :mod:`pyrpod.mdao.reference_data`). + """ + if self.results is None: + raise RuntimeError("run() the study before comparing results") + report = compare_results( + self.results.cases, dataset, + reference_length=self.config.loads.normalization.reference_length) + self.comparison = report + return report + + # ---------------------------------------------------------------- plots + def plot(self, results: StudyResults | None = None) -> list[str]: + """Generate the optional parameter-sweep trend plots. + + Plot generation is entirely optional -- automated tests never require + it -- and is imported lazily so a headless run pays no matplotlib + cost unless plots were asked for. + """ + from pyrpod.mdao import study_plots + + results = results or self.results + if results is None: + raise RuntimeError("run() the study before plotting results") + plot_dir = os.path.join(self.config.output_dir, + self.config.output.plots_subdir) + return study_plots.plot_sweep_trends(results, plot_dir, + comparison=self.comparison) From 97141dde0fdb6cc467d2b414490a3133742ef134 Mon Sep 17 00:00:00 2001 From: andytorrestb Date: Mon, 27 Jul 2026 22:54:22 -0500 Subject: [PATCH 10/14] mdao: add ParameterSweepStudy, a whole sweep in one JFH ParameterSweepStudy runs the sweep as ONE case driven by ONE Jet Firing History (sweep.mode: single_jfh): every pose contributes its firings to a single history, the strike pipeline runs once over all of them, and the strike files form a single results/strikes/firing-.vtu series -- the repository's existing sweep convention, scrubbable in ParaView as a time sequence. Because the history is shared, the pipeline's cumulative fields finally mean something: max_pressures / max_shears become the worst load any pose put on each face, cum_strikes a coverage map and cum_heat_flux_load an accumulated dose. Their per-component summary is exposed as `study.envelope` and recorded in the metadata; with the full pipeline disabled those arrays do not exist and the envelope is reported empty rather than reconstructed. A separate class rather than a mode flag inside the per-case engine: the two decompositions answer different questions and their outputs differ in kind, not just in layout. Everything they share lives in study_runtime, so per-firing numbers are identical between engines and only the decomposition differs. n_firings keeps its per-pose meaning, and the single history is required to hold exactly len(poses) * n_firings entries -- checked after writing and again after reading the file back. Co-Authored-By: Claude Opus 5 --- pyrpod/mdao/parameter_sweep.py | 237 +++++++++++++++++++++++++++++++++ 1 file changed, 237 insertions(+) create mode 100644 pyrpod/mdao/parameter_sweep.py diff --git a/pyrpod/mdao/parameter_sweep.py b/pyrpod/mdao/parameter_sweep.py new file mode 100644 index 0000000..40fe43b --- /dev/null +++ b/pyrpod/mdao/parameter_sweep.py @@ -0,0 +1,237 @@ +""" +Single-history engine for prescribed plume/target parameter sweeps. + +:class:`ParameterSweepStudy` runs a whole angle x distance sweep as ONE case +driven by ONE Jet Firing History: every pose contributes ``n_firings`` +consecutive entries to a single history, the strike pipeline runs once over +all of them, and every strike file therefore belongs to the same run. It is +the ``sweep.mode: single_jfh`` engine behind the package-level API + +>>> from pyrpod.mdao.TradeStudy import TradeStudy +>>> study = TradeStudy.from_config('sweep.yaml') # mode: single_jfh +>>> results = study.run() + +Why a second engine rather than a flag +-------------------------------------- +The two decompositions answer different questions, and their outputs differ +in kind, not just in layout: + +* :class:`pyrpod.mdao.plume_validation.PlumeValidationStudy` (per case) keeps + poses independent -- the right thing when each pose is a separate + experiment, as in a validation matrix compared pose-by-pose against + reference data; +* this class treats the sweep as a single firing history -- the right thing + when the poses form one sequence. The pipeline's cumulative fields then + mean something: ``max_pressures`` / ``max_shears`` become the worst load + each face saw ANYWHERE in the sweep, ``cum_strikes`` a coverage map, and + ``cum_heat_flux_load`` the accumulated dose. Those envelopes are recorded + in the study metadata and live in the VTK files, which land in one + ``results/strikes/firing-.vtu`` series that ParaView can scrub as a + time sequence. + +Everything else is shared with the per-case engine through +:mod:`pyrpod.mdao.study_runtime`: the same case objects, the same pose +generation, the same load integration, the same result schema, the same +reference comparison and plots. Per-firing numbers are identical between the +two engines -- only the decomposition into histories differs. + +``n_firings`` semantics +----------------------- +``n_firings`` stays the number of entries EACH POSE contributes, so the same +configuration can be run either way. The single history is required to hold +exactly ``len(poses) * n_firings`` entries, checked after writing and again +after reading the file back. +""" + +from __future__ import annotations + +import logging +import os +import time + +from pyrpod.mdao import firing_plan, study_runtime +from pyrpod.mdao.reference_data import ( + ComparisonReport, + ReferenceDataset, + compare_results, + load_reference_dataset, +) +from pyrpod.mdao.study_config import StudyConfig +from pyrpod.mdao.study_results import StudyResults +from pyrpod.mdao.study_runtime import TargetGeometry +from pyrpod.util.io.fs import ensure_dir + +logger = logging.getLogger(__name__) + +__all__ = ["ParameterSweepStudy"] + + +class ParameterSweepStudy: + """Prescribed parameter sweep executed as one Jet Firing History.""" + + def __init__(self, config: StudyConfig) -> None: + self.config = config + self.results: StudyResults | None = None + self.comparison: ComparisonReport | None = None + #: Worst-case per-face statistics over the whole sweep, per component + #: (populated by run() when the full pipeline path is used). + self.envelope: dict[str, dict[str, float]] = {} + + # ------------------------------------------------------------ lifecycle + @classmethod + def from_config(cls, path: str | os.PathLike[str] + ) -> "ParameterSweepStudy": + """Build a sweep study from a YAML configuration file.""" + return cls(StudyConfig.from_yaml(path)) + + @property + def case_id(self) -> str: + """The single case identifier every result row carries.""" + return self.config.study_name + + @property + def jfh_name(self) -> str: + """File name of the sweep's one Jet Firing History.""" + return f"{self.config.study_name}.A" + + # ------------------------------------------------------------------ run + def run(self, write_outputs: bool = True) -> StudyResults: + """Execute the whole sweep as one firing history. + + Returns + ------- + StudyResults + One :class:`~pyrpod.mdao.study_results.CaseResult` per firing and + component, each keyed to the pose it realizes (``plate_angle_deg`` + and ``source_distance``) and sharing one ``case_id``. + """ + config = self.config + started = time.perf_counter() + ensure_dir(config.output_dir) + + assets = study_runtime.load_case_assets(config) + geometry = TargetGeometry.from_config(config, assets.target_vehicle.mesh) + + # 1. One history for the whole sweep, exactly poses x n_firings long. + firings = firing_plan.build_sweep_firings(config.sweep, config.target) + jfh_path = os.path.join(config.output_dir, "jfh", self.jfh_name) + n_written = firing_plan.write_jfh_file(jfh_path, firings) + expected = config.sweep.total_firings + if n_written != expected: + raise ValueError( + f"sweep JFH {jfh_path!r} was written with {n_written} entries " + f"but {len(config.sweep.poses)} poses x n_firings=" + f"{config.sweep.n_firings} requires exactly {expected}") + + jfh = study_runtime.read_generated_jfh(config, self.jfh_name) + if len(jfh.JFH) != expected: + raise ValueError( + f"sweep JFH {jfh_path!r} reports {len(jfh.JFH)} firings, " + f"expected {expected}") + + logger.info("Parameter sweep study started: study=%s mode=single_jfh " + "poses=%d firings_per_pose=%d jfh_entries=%d " + "components=%d target_faces=%d", config.study_name, + len(config.sweep.poses), config.sweep.n_firings, expected, + len(geometry.components), geometry.n_faces) + + # 2. One strike run over every firing in that history. + firing_data, vtk_paths = study_runtime.compute_strikes( + config, assets, geometry, jfh, output_root=config.output_dir) + + # 3. One result record per firing and component, keyed to its pose. + results = StudyResults( + study_name=config.study_name, output_dir=config.output_dir, + provenance=study_runtime.study_provenance( + config, geometry, jfh_path=os.path.abspath(jfh_path), + case_id=self.case_id)) + timestamp = study_runtime.utc_timestamp() + version = str(results.provenance.get("code_version", "unknown")) + + for firing_index, firing in enumerate(firings): + per_face = firing_data[str(firing_index + 1)] + results.cases.extend(study_runtime.build_case_results( + config, geometry, firing, per_face, + case_id=self.case_id, firing_id=firing_index + 1, + plate_angle_deg=_pose_value(firing.plate_angle_deg), + source_distance=_pose_value(firing.source_distance), + jfh_path=jfh_path, + vtk_path=vtk_paths[firing_index] if vtk_paths else None, + code_version_id=version, timestamp=timestamp)) + + # 4. The sweep envelope, which only a shared history can produce. + last = firing_data[str(len(firings))] + if all(name in last for name in study_runtime.CUMULATIVE_FIELDS): + self.envelope = study_runtime.component_envelope(geometry, last) + results.provenance["sweep_envelope"] = self.envelope + else: + self.envelope = {} + logger.info("Sweep envelope unavailable: the cumulative fields " + "come from the full strike pipeline, which is " + "disabled when output.vtk.enabled is false.") + + if write_outputs: + self.write_outputs(results) + + self.results = results + if config.reference.path: + self.comparison = self.compare( + load_reference_dataset(config.reference.path, + label=config.reference.label)) + if write_outputs and len(self.comparison): + self.comparison.write_csv(os.path.join( + config.output_dir, "reference_comparison.csv")) + + logger.info("Parameter sweep study completed: study=%s firings=%d " + "records=%d wall_time_s=%.2f output=%s", + config.study_name, expected, len(results), + time.perf_counter() - started, config.output_dir) + return results + + # -------------------------------------------------------------- outputs + def write_outputs(self, results: StudyResults) -> None: + """Write the summary CSV, the JSON metadata and any configured plots.""" + config = self.config + results.write_csv(os.path.join(config.output_dir, + config.output.summary_csv)) + results.write_metadata(os.path.join(config.output_dir, + config.output.summary_metadata)) + if config.output.write_plots: + results.plot_paths = self.plot(results) + + # ----------------------------------------------------------- comparison + def compare(self, dataset: ReferenceDataset) -> ComparisonReport: + """Compare this sweep's results against a reference dataset. + + Records are matched on the swept parameters exactly as in the + per-case engine, so the same reference data compares against either. + """ + if self.results is None: + raise RuntimeError("run() the study before comparing results") + report = compare_results( + self.results.cases, dataset, + reference_length=self.config.loads.normalization.reference_length) + self.comparison = report + return report + + # ---------------------------------------------------------------- plots + def plot(self, results: StudyResults | None = None) -> list[str]: + """Generate the optional parameter-sweep trend plots.""" + from pyrpod.mdao import study_plots + + results = results or self.results + if results is None: + raise RuntimeError("run() the study before plotting results") + plot_dir = os.path.join(self.config.output_dir, + self.config.output.plots_subdir) + return study_plots.plot_sweep_trends(results, plot_dir, + comparison=self.comparison) + + +def _pose_value(value: float | None) -> float: + """A firing's swept parameter; NaN when its pose was prescribed outright. + + Explicitly prescribed firings that carry no sweep parameterization are + recorded as NaN rather than as a made-up angle or distance. + """ + return float("nan") if value is None else float(value) From e7b1a4c008f4e79a6483b2c8db5dd2f93d37de59 Mon Sep 17 00:00:00 2001 From: andytorrestb Date: Mon, 27 Jul 2026 22:54:42 -0500 Subject: [PATCH 11/14] mdao: repair TradeStudy and expose the package-level study API Repairs, in the legacy dynamics-driven sweep path -- each of which raised before any physics ran: - init_trade_study constructed PlumeStrikeEstimationStudy.RPOD(case_dir), a class that does not exist, and passed a case directory where the study takes a MissionEnvironment. It now builds the environment and the real PlumeStrikeEstimationStudy; - the sweeps called jfh_plume_strikes(trade_study=True), which that method's (parallel, workers) signature rejects with TypeError. They now use the per-firing arrays it returns; - print_mission_report read impingement maxima off study attributes that are never set. It now summarizes them from those returned arrays, keeping the old attribute path for existing callers; - interpret_mission_report read self.rpod.config, which does not exist; the configuration lives on the study's MissionEnvironment. Adds the package-level prescribed-validation API: study = TradeStudy.from_config('case.yaml') results = study.run() with compare() and plot() alongside it, and an optional output_dir override so a committed configuration can be run into a scratch location. The engine is selected from sweep.mode -- per_case binds PlumeValidationStudy, single_jfh binds ParameterSweepStudy -- so the API and the configuration file are the same either way. The class stays a thin facade over the pyrpod.mdao modules, and its legacy constructor keeps its original behavior. Co-Authored-By: Claude Opus 5 --- pyrpod/mdao/TradeStudy.py | 232 +++++++++++++++++++++++++++++++++----- 1 file changed, 205 insertions(+), 27 deletions(-) diff --git a/pyrpod/mdao/TradeStudy.py b/pyrpod/mdao/TradeStudy.py index 89890b3..7be2e94 100644 --- a/pyrpod/mdao/TradeStudy.py +++ b/pyrpod/mdao/TradeStudy.py @@ -1,3 +1,25 @@ +""" +Trade studies over RPOD / plume-impingement configurations. + +Two entry points live here: + +* the legacy dynamics-driven sweeps (``run_axial_overshoot_sweep``, + ``run_surface_cant_sweep``, ``run_multi_var_sweep``), which vary an + approach velocity and/or a thruster cant angle and report a mission + summary per configuration; +* the package-level PRESCRIBED validation API, + + >>> study = TradeStudy.from_config('flat_plate_baseline.yaml') + >>> results = study.run() + + which runs a YAML-configured plume/target validation sweep through + :class:`pyrpod.mdao.plume_validation.PlumeValidationStudy` and returns + structured results. The heavy lifting lives in the small modules of + ``pyrpod.mdao`` (study configuration, firing plan, surface-load + integration, results, reference comparison, plots), so this class stays a + thin façade over them. +""" + from __future__ import annotations import os @@ -8,11 +30,22 @@ from typing import Any import numpy as np +from numpy.typing import NDArray import pandas as pd import matplotlib.pyplot as plt from pyrpod.rpod import JetFiringHistory, PlumeStrikeEstimationStudy from pyrpod.mdao import SweepConfig +from pyrpod.mdao.parameter_sweep import ParameterSweepStudy +from pyrpod.mdao.plume_validation import PlumeValidationStudy +from pyrpod.mdao.reference_data import ( + ComparisonReport, + ReferenceDataset, + load_reference_dataset, +) +from pyrpod.mdao.study_config import StudyConfig +from pyrpod.mdao.study_results import StudyResults +from pyrpod.mission import MissionEnvironment from pyrpod.util.io.fs import ensure_dir import configparser @@ -22,6 +55,21 @@ logger = logging.getLogger(__name__) +def _engine_for( + study_config: StudyConfig, +) -> PlumeValidationStudy | ParameterSweepStudy: + """Select the study engine named by the configuration's sweep mode. + + ``per_case`` (the default) gives every angle-distance combination its own + Jet Firing History and strike run; ``single_jfh`` runs the whole sweep as + one history, so every firing's strikes come from the same pipeline run + and the cumulative fields form a sweep envelope. + """ + if study_config.sweep.mode == "single_jfh": + return ParameterSweepStudy(study_config) + return PlumeValidationStudy(study_config) + + class TradeStudy(): # Built by init_trade_study and read by every sweep; max_v0 is recorded # by the sweeps and read back by print_mission_report. Declared bare so @@ -29,31 +77,107 @@ class TradeStudy(): rpod: Any max_v0: Any - def __init__(self, case_dir: str) -> None: + def __init__(self, case_dir: str, + study_config: StudyConfig | None = None) -> None: self.case_dir = case_dir config = configparser.ConfigParser() config.read(self.case_dir + "config.ini") self.config = config + # Present only for a study built by from_config(); the legacy + # constructor keeps its original single-argument behavior. The engine + # is selected by the configuration's sweep mode: one Jet Firing + # History per case, or one spanning the whole sweep. + self.study_config = study_config + self._validation_study = (_engine_for(study_config) + if study_config is not None else None) + + # ------------------------------------------------------- package API + @classmethod + def from_config(cls, path: str | os.PathLike[str], + output_dir: str | os.PathLike[str] | None = None, + ) -> "TradeStudy": + """Build a prescribed validation study from a YAML configuration. + + Parameters + ---------- + path : str or path-like + Path to a study configuration file (see + :class:`pyrpod.mdao.study_config.StudyConfig` and + ``docs/plume_validation_study.md``). + output_dir : str or path-like, optional + Overrides the configured output directory, so the same committed + configuration can be run into a scratch location (this is what + the automated tests do). + + Returns + ------- + TradeStudy + A study bound to the parsed configuration; call :meth:`run`. + Which engine backs it follows ``sweep.mode``: + ``per_case`` (default) binds + :class:`~pyrpod.mdao.plume_validation.PlumeValidationStudy`, + ``single_jfh`` binds + :class:`~pyrpod.mdao.parameter_sweep.ParameterSweepStudy`. + """ + config = StudyConfig.from_yaml(path) + if output_dir is not None: + config = config.with_output_dir(os.fspath(output_dir)) + return cls(config.case_dir, study_config=config) + + @property + def validation_study(self) -> PlumeValidationStudy | ParameterSweepStudy: + """The bound study engine (requires from_config).""" + if self._validation_study is None: + raise RuntimeError( + "this TradeStudy was not built from a study configuration; " + "use TradeStudy.from_config(path) for prescribed plume " + "validation sweeps") + return self._validation_study + + def run(self, write_outputs: bool = True) -> StudyResults: + """Run the configured sweep through its selected engine.""" + return self.validation_study.run(write_outputs=write_outputs) + + def compare(self, reference: ReferenceDataset | str | os.PathLike[str] + ) -> ComparisonReport: + """Compare the last run against external reference data. + + The reference may be a loaded :class:`ReferenceDataset` or a path to + a CSV / JSON / YAML file in the generic reference format; where the + data came from (DSMC, analytical, experiment) is irrelevant here. + """ + if not isinstance(reference, ReferenceDataset): + reference = load_reference_dataset(reference) + return self.validation_study.compare(reference) + + def plot(self) -> list[str]: + """Generate the optional sweep trend plots for the last run.""" + return self.validation_study.plot() + # ------------------------------------------------------ legacy sweeps def init_trade_study(self, lm: LogisticsModule, tv: TargetVehicle) -> None: """ Organizes data needed to kick off an RPOD trade study. - Mainly done by properly configuring an RPOD object + Mainly done by properly configuring a plume-strike study object. """ - # Save variable name for readability. + # Save variable name for readability. case_dir = self.case_dir # Instantiate JetFiringHistory object. jfh = JetFiringHistory.JetFiringHistory(case_dir) - # Instantiate RPOD object. - # The class in that module is PlumeStrikeEstimationStudy; there is no - # RPOD, so this raises AttributeError. Renaming it is a behavior - # change (see deferred observations), so it is flagged here. - rpod = PlumeStrikeEstimationStudy.RPOD(case_dir) # type: ignore[attr-defined] + # Instantiate the plume-strike study object. It takes a + # MissionEnvironment (MissionPlanner's constructor argument), not a + # case directory; the previous call named a class that does not + # exist in the module and raised AttributeError before any sweep + # could start. + environment = MissionEnvironment.MissionEnvironment(case_dir) + rpod = PlumeStrikeEstimationStudy.PlumeStrikeEstimationStudy( + environment) rpod.study_init(jfh, tv, lm) + self.environment = environment self.rpod = rpod def init_trade_study_case(self) -> None: @@ -70,19 +194,66 @@ def init_trade_study_case(self) -> None: self.rpod.jfh = jfh - def print_mission_report(self) -> None: + @staticmethod + def summarize_firing_data( + firing_data: Mapping[str, Mapping[str, NDArray[np.float64]]] + ) -> dict[str, float]: + """Reduce a plume-strike run's per-firing arrays to case maxima. + + ``jfh_plume_strikes`` returns the per-face arrays for every firing + (including the running cumulative ones); the mission report wants one + number per quantity per configuration, so they are reduced here + instead of being read off study attributes that the plume-strike + study never sets. + """ + summary = {"max_pressure": 0.0, "max_shear": 0.0, + "max_heat_rate": 0.0, "max_heat_load": 0.0, + "max_cum_heat_load": 0.0} + for fields in firing_data.values(): + for key, field_name in (("max_pressure", "pressures"), + ("max_shear", "shear_stress"), + ("max_heat_rate", "heat_flux_rate"), + ("max_heat_load", "heat_flux_load"), + ("max_cum_heat_load", + "cum_heat_flux_load")): + values = fields.get(field_name) + if values is None or len(values) == 0: + continue + summary[key] = max(summary[key], float(np.max(values))) + return summary + + def print_mission_report( + self, + firing_data: Mapping[str, Mapping[str, NDArray[np.float64]]] | None + = None, + ) -> None: + """Append one configuration's summary row to results/MissionReport.csv. + + ``firing_data`` is the mapping returned by + ``PlumeStrikeEstimationStudy.jfh_plume_strikes()``; when supplied the + impingement maxima are reduced from it. It is optional so an existing + caller that sets the maxima on the study object itself keeps working. + """ case_key = self.rpod.get_case_key() max_v0 = self.max_v0 - fuel_mass = self.rpod.fuel_mass - - max_pressure = self.rpod.max_pressure - max_shear = self.rpod.max_shear - max_heat_rate = self.rpod.max_heat_rate - max_heat_load = self.rpod.max_heat_load - max_cum_heat_load = self.rpod.max_cum_heat_load + fuel_mass = getattr(self.rpod, 'fuel_mass', '') + + if firing_data is not None: + summary = self.summarize_firing_data(firing_data) + max_pressure = summary['max_pressure'] + max_shear = summary['max_shear'] + max_heat_rate = summary['max_heat_rate'] + max_heat_load = summary['max_heat_load'] + max_cum_heat_load = summary['max_cum_heat_load'] + else: + max_pressure = self.rpod.max_pressure + max_shear = self.rpod.max_shear + max_heat_rate = self.rpod.max_heat_rate + max_heat_load = self.rpod.max_heat_load + max_cum_heat_load = self.rpod.max_cum_heat_load # Check if the file exists report_path = self.case_dir + 'results/MissionReport.csv' @@ -139,11 +310,15 @@ def interpret_mission_report(self) -> None: report_path = self.case_dir + 'results/MissionReport.csv' report_results = pd.read_csv(report_path) - max_pressure = float(self.rpod.config['tv']['normal_pressure']) - max_shear = float(self.rpod.config['tv']['shear_pressure']) - max_heat_rate = float(self.rpod.config['tv']['heat_flux']) - # max_heat_load = float(self.rpod.config['tv'][]) no such constraint - max_cum_heat_load = float(self.rpod.config['tv']['heat_flux_load']) + # The plume-strike study holds its configuration on its + # MissionEnvironment; `self.rpod.config` never existed and raised + # AttributeError here. + config = self.rpod.environment.config + max_pressure = float(config['tv']['normal_pressure']) + max_shear = float(config['tv']['shear_pressure']) + max_heat_rate = float(config['tv']['heat_flux']) + # max_heat_load = float(config['tv'][]) no such constraint + max_cum_heat_load = float(config['tv']['heat_flux_load']) plume_status = [] plume_failure_mode = [] @@ -215,9 +390,12 @@ def run_axial_overshoot_sweep(self, sweep_vars: Mapping[str, Any], # Reset JFH according to specific case. self.init_trade_study_case() self.rpod.graph_jfh(trade_study= True) - self.rpod.jfh_plume_strikes(trade_study = True) + # jfh_plume_strikes takes (parallel, workers); the previous + # trade_study keyword raised TypeError. The per-firing arrays it + # returns are what the mission report summarizes. + firing_data = self.rpod.jfh_plume_strikes() - self.print_mission_report() + self.print_mission_report(firing_data) logger.info("Trade-study config %d/%d completed: case_key=%s " "config_time_s=%.3f", i + 1, len(axial_overshoot), self.rpod.get_case_key(), @@ -275,9 +453,9 @@ def run_surface_cant_sweep(self, sweep_vars: Mapping[str, Any], self.init_trade_study_case() self.rpod.graph_jfh(trade_study= True) - self.rpod.jfh_plume_strikes(trade_study = True) + firing_data = self.rpod.jfh_plume_strikes() - self.print_mission_report() + self.print_mission_report(firing_data) self.interpret_mission_report() def run_multi_var_sweep(self, sweep_vars: Mapping[str, Any], @@ -325,7 +503,7 @@ def run_multi_var_sweep(self, sweep_vars: Mapping[str, Any], self.init_trade_study_case() self.rpod.graph_jfh(trade_study= True) - self.rpod.jfh_plume_strikes(trade_study = True) + firing_data = self.rpod.jfh_plume_strikes() - self.print_mission_report() + self.print_mission_report(firing_data) self.interpret_mission_report() \ No newline at end of file From 8ec09a100d4b5c369d68421a5bf7c2d0318b55d9 Mon Sep 17 00:00:00 2001 From: andytorrestb Date: Mon, 27 Jul 2026 22:54:42 -0500 Subject: [PATCH 12/14] case: add worked study configurations Four examples on existing case assets, all recording the plume model explicitly and taking moments about the target's own reference point: - flat_plate_baseline.yaml: one head-on firing at L/D = 4; - flat_plate_sweep.yaml: 19 angles x 5 distances as 95 independent cases; - flat_plate_sweep_single_jfh.yaml: the same sweep as one 95-firing history, differing from the previous file only in sweep.mode and the output names; - cylinder_baseline.yaml: the same machinery on a curved, closed target, to keep the architecture honest and to be where quantitative cylinder reference data plugs in later. The flat-plate files carry the Cai normalization derived from the case's own thruster definition file, so coefficients are available. The cylinder file supplies none, so its coefficients are correctly reported as unavailable rather than invented -- no cylinder reference data exists yet. Co-Authored-By: Claude Opus 5 --- .../study/cylinder_baseline.yaml | 69 ++++++++++++++ .../study/flat_plate_baseline.yaml | 86 ++++++++++++++++++ .../study/flat_plate_sweep.yaml | 76 ++++++++++++++++ .../study/flat_plate_sweep_single_jfh.yaml | 89 +++++++++++++++++++ 4 files changed, 320 insertions(+) create mode 100644 case/plume/plume_cylinder_sweep/study/cylinder_baseline.yaml create mode 100644 case/plume/plume_flat_plate_sweep/study/flat_plate_baseline.yaml create mode 100644 case/plume/plume_flat_plate_sweep/study/flat_plate_sweep.yaml create mode 100644 case/plume/plume_flat_plate_sweep/study/flat_plate_sweep_single_jfh.yaml diff --git a/case/plume/plume_cylinder_sweep/study/cylinder_baseline.yaml b/case/plume/plume_cylinder_sweep/study/cylinder_baseline.yaml new file mode 100644 index 0000000..250164e --- /dev/null +++ b/case/plume/plume_cylinder_sweep/study/cylinder_baseline.yaml @@ -0,0 +1,69 @@ +# PyRPOD prescribed plume-validation study: CYLINDER target, one firing. +# +# The same study machinery as the flat-plate examples, pointed at the curved, +# closed target of case/plume/plume_cylinder_sweep (radius 2 m, axis along X +# spanning x in [-7, 0], centroid (-3.5, 0, 0)). It exists to keep the study +# architecture honest -- nothing in it assumes a flat plate -- and to be the +# place quantitative cylinder reference data plugs into later, through the +# same generic reference-comparison interface the plate cases use. +# +# NOTE: as in the case itself, the cylinder results are a PIPELINE SMOKE +# CASE, not a validated physical answer: no cylinder reference data exists +# yet. No normalization inputs are supplied, so coefficients are correctly +# reported as unavailable rather than invented. + +study: + name: cylinder_baseline + description: >- + Single head-on firing against the cylinder target; architecture check and + future home of quantitative cylinder validation. + case_dir: .. + output_dir: ../results/studies/cylinder_baseline + +thruster: + id: T1 + +plume_model: + name: SimplifiedGasKinetics + parameters: + gas: argon + speed_ratio_S0: 2.0 + stagnation_temperature_K: 200.0 + nozzle_diameter_m: 1.0 + +target: + geometry_id: cylinder.stl + # Cylinder centroid; the sweep arc is built about it, exactly as the + # case's own jfh/generate_cylinder_sweep_jfh.py does. + reference_point: [-3.5, 0.0, 0.0] + normal: [0.0, 0.0, 1.0] + tangent: [1.0, 0.0, 0.0] + components: + - name: cylinder + selector: all + +sweep: + plate_angles_deg: [0.0] + # Orbit radius; larger than the ~4.03 m bounding sphere of the target. + source_distances: [6.0] + n_firings: 1 + firing_duration_s: 1.0 + thrusters: [1] + +loads: + moment_reference_point: [-3.5, 0.0, 0.0] + # No normalization block: coefficients are reported as unavailable. + +output: + vtk: + enabled: true + summary: + csv: case_results.csv + metadata: study_metadata.json + plots: + enabled: false + +metadata: + coordinate_system: >- + case global frame; cylinder axis along X, plume source swept in the X-Z + plane about the cylinder centroid diff --git a/case/plume/plume_flat_plate_sweep/study/flat_plate_baseline.yaml b/case/plume/plume_flat_plate_sweep/study/flat_plate_baseline.yaml new file mode 100644 index 0000000..3625a3c --- /dev/null +++ b/case/plume/plume_flat_plate_sweep/study/flat_plate_baseline.yaml @@ -0,0 +1,86 @@ +# PyRPOD prescribed plume-validation study: flat-plate BASELINE case. +# +# One head-on firing of the Cai 2016 round argon jet (D = 1 m, S0 = 2.0, +# T0 = 200 K, Tw = 300 K, fully diffuse) against the flat 8 m x 8 m plate of +# case/plume/plume_flat_plate_sweep, with the plume source at L = 4D on the +# plate normal. This is the single-pose reference point of the sweep in +# flat_plate_sweep.yaml. +# +# Run it: +# from pyrpod.mdao.TradeStudy import TradeStudy +# results = TradeStudy.from_config( +# 'case/plume/plume_flat_plate_sweep/study/flat_plate_baseline.yaml' +# ).run() +# +# The case's own config.ini still owns every asset (STL, TCF, TDF, plume +# model and gating geometry); this file adds only the study layer. + +study: + name: cai2016_flat_plate_baseline + description: >- + Head-on flat-plate plume impingement at L/D = 4, Cai 2016 Section 4 + conditions, run through the prescribed trade-study API. + case_dir: .. + output_dir: ../results/studies/flat_plate_baseline + +# One thruster, recorded explicitly. It must exist in the case's TCF. +thruster: + id: T1 + +# The only plume model this workflow supports. +plume_model: + name: SimplifiedGasKinetics + parameters: + # Recorded for provenance; the values themselves are read from the + # case's thruster definition file (tcd/tdf.csv, thruster type ARG). + gas: argon + speed_ratio_S0: 2.0 + stagnation_temperature_K: 200.0 + nozzle_diameter_m: 1.0 + +target: + geometry_id: flat_plate_transformed.stl + # Plate center; the sweep arc and the angle convention are built about it. + reference_point: [0.0, 0.0, 0.0] + # Outward normal (toward the plume-source side) and in-plane tilt axis. + normal: [0.0, 0.0, 1.0] + tangent: [1.0, 0.0, 0.0] + components: + - name: plate + selector: all + +sweep: + plate_angles_deg: [0.0] + source_distances: [4.0] + # Exactly this many entries are written to the case's Jet Firing History. + n_firings: 1 + firing_duration_s: 1.0 + thrusters: [1] + +loads: + # Moments are reported about this user-defined point (the plate center). + moment_reference_point: [0.0, 0.0, 0.0] + normalization: + # Plate area 8 m x 8 m and semi-length H0 = 4 m. + reference_area: 64.0 + reference_length: 4.0 + # Cai's normalization, derived from the case's own TDF entry: + # q_dyn = n0*m*U0^2/2 with n0 = 1e20 m^-3, U0 = 577.0684534784414 m/s, + # m = (R_universal / R_specific) / N_A, R_specific = 208.13 + # q_heat = n0*m*U0^3/2 + dynamic_pressure: 1.1044652197738332 + reference_heat_flux: 637.3520362956127 + +output: + vtk: + enabled: true + summary: + csv: case_results.csv + metadata: study_metadata.json + plots: + enabled: false + +metadata: + coordinate_system: >- + case global frame; plate in the X-Y plane centered at the origin, + plume source on the +Z side diff --git a/case/plume/plume_flat_plate_sweep/study/flat_plate_sweep.yaml b/case/plume/plume_flat_plate_sweep/study/flat_plate_sweep.yaml new file mode 100644 index 0000000..bfc7241 --- /dev/null +++ b/case/plume/plume_flat_plate_sweep/study/flat_plate_sweep.yaml @@ -0,0 +1,76 @@ +# PyRPOD prescribed plume-validation study: flat-plate ANGLE x DISTANCE sweep. +# +# The multi-angle / multi-distance analog of flat_plate_baseline.yaml: the +# stationary 8 m x 8 m plate of case/plume/plume_flat_plate_sweep is struck +# from 19 approach angles (-90 .. +90 deg in 10 deg steps, 0 = head-on) at 5 +# stand-off distances (L/D in {2, 4, 6, 8, 10}, D = 1 m) -- 95 cases, one +# firing each, matching the committed sweep JFH of that case. +# +# Run it: +# from pyrpod.mdao.TradeStudy import TradeStudy +# study = TradeStudy.from_config( +# 'case/plume/plume_flat_plate_sweep/study/flat_plate_sweep.yaml') +# results = study.run() +# study.plot() # optional trend figures +# +# Expect ~10368 faces x 95 poses; with VTK output enabled this writes one +# .vtu per case under results/studies/flat_plate_sweep/cases/. + +study: + name: cai2016_flat_plate_sweep + description: >- + Flat-plate plume impingement swept over 19 approach angles and 5 + stand-off distances, Cai 2016 Section 4 conditions. + case_dir: .. + output_dir: ../results/studies/flat_plate_sweep + +thruster: + id: T1 + +plume_model: + name: SimplifiedGasKinetics + parameters: + gas: argon + speed_ratio_S0: 2.0 + stagnation_temperature_K: 200.0 + nozzle_diameter_m: 1.0 + +target: + geometry_id: flat_plate_transformed.stl + reference_point: [0.0, 0.0, 0.0] + normal: [0.0, 0.0, 1.0] + tangent: [1.0, 0.0, 0.0] + components: + - name: plate + selector: all + +sweep: + plate_angles_deg: [-90.0, -80.0, -70.0, -60.0, -50.0, -40.0, -30.0, -20.0, + -10.0, 0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, + 80.0, 90.0] + source_distances: [2.0, 4.0, 6.0, 8.0, 10.0] + n_firings: 1 + firing_duration_s: 1.0 + thrusters: [1] + +loads: + moment_reference_point: [0.0, 0.0, 0.0] + normalization: + reference_area: 64.0 + reference_length: 4.0 + dynamic_pressure: 1.1044652197738332 + reference_heat_flux: 637.3520362956127 + +output: + vtk: + enabled: true + summary: + csv: sweep_results.csv + metadata: sweep_metadata.json + plots: + enabled: true + +metadata: + coordinate_system: >- + case global frame; plate in the X-Y plane centered at the origin, + plume source swept in the X-Z plane diff --git a/case/plume/plume_flat_plate_sweep/study/flat_plate_sweep_single_jfh.yaml b/case/plume/plume_flat_plate_sweep/study/flat_plate_sweep_single_jfh.yaml new file mode 100644 index 0000000..76b8640 --- /dev/null +++ b/case/plume/plume_flat_plate_sweep/study/flat_plate_sweep_single_jfh.yaml @@ -0,0 +1,89 @@ +# PyRPOD prescribed plume-validation study: flat-plate sweep as ONE JFH. +# +# Same physics, same poses and same normalization as flat_plate_sweep.yaml -- +# 19 approach angles x 5 stand-off distances against the 8 m x 8 m plate of +# case/plume/plume_flat_plate_sweep -- but executed as a SINGLE case driven by +# a SINGLE Jet Firing History holding all 95 firings (sweep.mode: single_jfh). +# +# What that buys, over the per-case form: +# * one strike run, so every firing's results come from the same pipeline +# invocation and the same cumulative accumulation; +# * one results/strikes/firing-.vtu series (i = 0 .. 94) that ParaView +# can scrub as a time sequence, matching the case's own sweep convention; +# * a sweep ENVELOPE in the metadata: per component, the worst pressure and +# shear any pose produced on any face, the accumulated heat-flux load, the +# coverage (unique struck faces / swept affected area). +# +# n_firings stays the count per POSE, so this file and flat_plate_sweep.yaml +# differ only in `mode` and the output names: 19 x 5 poses x 1 firing = a JFH +# of exactly 95 entries. +# +# Run it: +# from pyrpod.mdao.TradeStudy import TradeStudy +# study = TradeStudy.from_config( +# 'case/plume/plume_flat_plate_sweep/study/' +# 'flat_plate_sweep_single_jfh.yaml') +# results = study.run() # -> ParameterSweepStudy +# study.validation_study.envelope['plate'] + +study: + name: cai2016_flat_plate_sweep_single_jfh + description: >- + Flat-plate plume impingement over 19 approach angles and 5 stand-off + distances, executed as one 95-firing Jet Firing History. + case_dir: .. + output_dir: ../results/studies/flat_plate_sweep_single_jfh + +thruster: + id: T1 + +plume_model: + name: SimplifiedGasKinetics + parameters: + gas: argon + speed_ratio_S0: 2.0 + stagnation_temperature_K: 200.0 + nozzle_diameter_m: 1.0 + +target: + geometry_id: flat_plate_transformed.stl + reference_point: [0.0, 0.0, 0.0] + normal: [0.0, 0.0, 1.0] + tangent: [1.0, 0.0, 0.0] + components: + - name: plate + selector: all + +sweep: + # One Jet Firing History for the whole sweep. + mode: single_jfh + plate_angles_deg: [-90.0, -80.0, -70.0, -60.0, -50.0, -40.0, -30.0, -20.0, + -10.0, 0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, + 80.0, 90.0] + source_distances: [2.0, 4.0, 6.0, 8.0, 10.0] + # Entries contributed by each pose; the history holds 19 * 5 * 1 = 95. + n_firings: 1 + firing_duration_s: 1.0 + thrusters: [1] + +loads: + moment_reference_point: [0.0, 0.0, 0.0] + normalization: + reference_area: 64.0 + reference_length: 4.0 + dynamic_pressure: 1.1044652197738332 + reference_heat_flux: 637.3520362956127 + +output: + vtk: + enabled: true + summary: + csv: sweep_results.csv + metadata: sweep_metadata.json + plots: + enabled: true + +metadata: + coordinate_system: >- + case global frame; plate in the X-Y plane centered at the origin, + plume source swept in the X-Z plane From 0b2e61823d7b3df0566990f14180ca9efe6f72ba Mon Sep 17 00:00:00 2001 From: andytorrestb Date: Mon, 27 Jul 2026 22:58:29 -0500 Subject: [PATCH 13/14] tests: cover the prescribed plume-validation study Unit level: - mdao_unit_test_03, the exact meaning of n_firings: invalid counts are rejected; one requested firing produces one JFH entry and N produce exactly N, read back through JetFiringHistory; an explicit firing list that disagrees is an error; a whole-sweep sequence holds exactly poses x n_firings pose-tagged entries with continuous firing times; the generated poses reproduce all 95 firings of the committed flat-plate sweep JFH to file precision; and the dynamics-driven approach honors an exact count, reporting an unreachable one rather than shortening the history. - mdao_unit_test_04, surface-load integration on meshes whose loading has a known closed-form resultant: force, moment about a user-defined point, center of pressure (recovered from a linear pressure field, consistent with the reported moment, and unavailable with a stated reason for zero-load and cancelling loads), coefficients present and absent, component selection. - mdao_unit_test_05, configuration parsing and validation, including that the case's own config.ini keeps parsing unchanged and that n_firings keeps its per-pose meaning in both sweep modes. - mdao_unit_test_06, the comparison metrics and CSV / JSON / YAML datasets (compared identically regardless of origin), plus the result schema's serialization. Integration level: - mdao_integration_test_02, the baseline flat-plate case end to end through TradeStudy.from_config: artifacts, provenance, head-on physics, available coefficients, and agreement of the integrated normal load with the INDEPENDENT Cai 2016 exact reference, not with PyRPOD's own output. - mdao_integration_test_03, a reduced multi-angle / multi-distance sweep: mirror symmetry in +/- angle, load decay with distance, center-of-pressure travel, per-case VTK isolation, optional plots, and the same machinery on the cylinder target with coefficients correctly unavailable. - mdao_integration_test_04, the single-history engine: engine dispatch from the configuration, one JFH of exactly poses x n_firings entries, one record per firing keyed to its pose, a single results/strikes series, and the sweep envelope (absent, not fabricated, without the full pipeline). Its strongest assertion is equivalence -- per-firing loads match the per-case engine pose by pose. Measured against the exact reference: integrated normal load within 1.5% mean / 2.2% max, component heat load ~15% -- the documented accuracy of the Maxwellian engineering chain, asserted only as a loose envelope so a convention regression fails loudly while the known gap passes. Studies write into temporary directories, so a test run leaves no artifacts behind. The manifest entries and the regenerated tests/README.md land here rather than in a follow-up commit: the tooling suite checks that the manifest and pytest's collection agree, so the two must move together. Co-Authored-By: Claude Opus 5 --- tests/README.md | 13 +- tests/mdao/mdao_integration_test_02.py | 230 +++++++++++++++ tests/mdao/mdao_integration_test_03.py | 316 +++++++++++++++++++++ tests/mdao/mdao_integration_test_04.py | 255 +++++++++++++++++ tests/mdao/mdao_unit_test_03.py | 350 +++++++++++++++++++++++ tests/mdao/mdao_unit_test_04.py | 370 +++++++++++++++++++++++++ tests/mdao/mdao_unit_test_05.py | 351 +++++++++++++++++++++++ tests/mdao/mdao_unit_test_06.py | 324 ++++++++++++++++++++++ tests/test_manifest.yaml | 125 +++++++++ 9 files changed, 2331 insertions(+), 3 deletions(-) create mode 100644 tests/mdao/mdao_integration_test_02.py create mode 100644 tests/mdao/mdao_integration_test_03.py create mode 100644 tests/mdao/mdao_integration_test_04.py create mode 100644 tests/mdao/mdao_unit_test_03.py create mode 100644 tests/mdao/mdao_unit_test_04.py create mode 100644 tests/mdao/mdao_unit_test_05.py create mode 100644 tests/mdao/mdao_unit_test_06.py diff --git a/tests/README.md b/tests/README.md index 0ff3075..994cd47 100644 --- a/tests/README.md +++ b/tests/README.md @@ -32,9 +32,9 @@ Source of truth for the metadata below: [`test_manifest.yaml`](test_manifest.yam | Metric | Count | | --- | --- | -| Manifest entries | 93 | -| Collected by pytest (files) | 44 | -| Collected by pytest (test cases) | 204 | +| Manifest entries | 100 | +| Collected by pytest (files) | 51 | +| Collected by pytest (test cases) | 334 | | Manual verification scripts | 42 | | Placeholder tests | 12 | | Blocked tests | 5 | @@ -71,12 +71,19 @@ separately). | --- | --- | --- | --- | --- | | [`mdao_unit_test_01.py`](mdao/mdao_unit_test_01.py) | 1 | Placeholder MDAO unit test; the test body returns immediately and asserts nothing. | `placeholder` | — | | [`mdao_unit_test_02.py`](mdao/mdao_unit_test_02.py) | 1 | Intended to build an array of cant-angle-swept thruster configurations (symmetric pitch/yaw canting) and visualize each sweep step. The entire body is currently commented out, so the test asserts nothing. | `placeholder` | — | +| [`mdao_unit_test_03.py`](mdao/mdao_unit_test_03.py) | 23 | Unit tests for prescribed firing generation and the exact meaning of n_firings: invalid counts are rejected, N requested firings produce exactly N JFH entries (read back through JetFiringHistory), an explicit firing list that disagrees with n_firings is an error, a whole-sweep sequence holds exactly poses x n_firings pose-tagged entries, the generated pose convention reproduces the committed flat-plate sweep JFH, and the dynamics-driven approach honors an exact count. | `implemented` | — | +| [`mdao_unit_test_04.py`](mdao/mdao_unit_test_04.py) | 20 | Unit tests for per-component surface-load integration on meshes with analytically known loading: pressure and shear force integration, moments about a user-defined reference point, center of pressure (recovered, moment-consistent, and unavailable for zero-load and near-cancellation cases), coefficient calculation and omission, and component face selection. | `implemented` | — | +| [`mdao_unit_test_05.py`](mdao/mdao_unit_test_05.py) | 32 | Unit tests for the YAML study configuration: the committed flat-plate examples parse, paths resolve against the configuration file, validation rejects an unsupported plume model, an unknown sweep mode, invalid firing counts, mismatched firing lists and bad geometry or normalization inputs, the per-pose n_firings semantics hold in both sweep modes, and the case's own config.ini keeps parsing unchanged. | `implemented` | — | +| [`mdao_unit_test_06.py`](mdao/mdao_unit_test_06.py) | 22 | Unit tests for the generic external-reference comparison (absolute and relative error, normalized RMSE, peak error, integrated-load error, center-of-pressure displacement; CSV / JSON / YAML datasets compared identically regardless of origin; unmatched cases and missing quantities reported rather than fabricated) and for the structured result schema's CSV and JSON serialization. | `implemented` | — | #### Integration | Test file | Cases | Description | Development status | Reference | | --- | --- | --- | --- | --- | | [`mdao_integration_test_01.py`](mdao/mdao_integration_test_01.py) | 1 | Placeholder MDAO integration test; the test body returns immediately and asserts nothing. | `placeholder` | — | +| [`mdao_integration_test_02.py`](mdao/mdao_integration_test_02.py) | 8 | Baseline flat-plate case run end to end through the package-level TradeStudy.from_config API: one case with exactly one JFH entry, the standard per-face strike VTK at the advertised path, CSV and JSON summaries with full provenance, physically consistent head-on loads, available coefficients, and agreement of the integrated normal load with the independent Cai 2016 exact reference. | `implemented` | Cai 2016, Aerospace 3(4), 43, doi:10.3390/aerospace3040043, Eq. 15 | +| [`mdao_integration_test_03.py`](mdao/mdao_integration_test_03.py) | 13 | Multi-angle / multi-distance flat-plate sweep through the trade-study API (a reduced subset of the committed sweep configuration): one result per angle-distance case, mirror symmetry in +/- angle, load decay with distance, center-of-pressure travel, agreement with the independent Cai 2016 reference through the generic comparison interface, per-case VTK isolation, optional trend plots, and the same machinery running on the cylinder target with coefficients correctly unavailable. | `implemented` | Cai 2016, Aerospace 3(4), 43, doi:10.3390/aerospace3040043, Eq. 15 | +| [`mdao_integration_test_04.py`](mdao/mdao_integration_test_04.py) | 12 | The single-history sweep engine (sweep.mode: single_jfh): engine dispatch from the configuration, exactly one Jet Firing History holding poses x n_firings entries, one result record per firing keyed to its pose, per-firing equivalence with the per-case engine, a single results/strikes VTK series, the per-component sweep envelope (absent rather than fabricated without the full pipeline), and the same independent Cai 2016 reference comparison. | `implemented` | Cai 2016, Aerospace 3(4), 43, doi:10.3390/aerospace3040043, Eq. 15 | #### Verification diff --git a/tests/mdao/mdao_integration_test_02.py b/tests/mdao/mdao_integration_test_02.py new file mode 100644 index 0000000..2dd6570 --- /dev/null +++ b/tests/mdao/mdao_integration_test_02.py @@ -0,0 +1,230 @@ +# ======================== +# PyRPOD: tests/mdao/mdao_integration_test_02.py +# ======================== +# Baseline flat-plate case run end to end through the package-level trade +# study API: +# +# study = TradeStudy.from_config('.../flat_plate_baseline.yaml') +# results = study.run() +# +# The case is the Cai 2016 Section-4 geometry reframed flat +# (case/plume/plume_flat_plate_sweep): argon round jet, D = 1 m, S0 = 2.0, +# T0 = 200 K, Tw = 300 K, fully diffuse, 8 m x 8 m plate, source head-on at +# L = 4D. Checked here: +# +# * exactly one case x component x firing, with a Jet Firing History +# holding exactly the requested single entry; +# * the standard per-face strike VTK is written at the path the result +# record advertises, carrying the pipeline's own cell fields; +# * the machine-readable summary (CSV) and metadata (JSON) are written and +# carry enough provenance to reproduce the run; +# * the integrated loads are physically right for a head-on pose: force +# along the plume direction, no in-plane resultant, center of pressure at +# the plate center, every face struck; +# * the integrated normal load agrees with the INDEPENDENT Cai 2016 exact +# reference (pyrpod/plume/CaiImpingement2016.py, Eq. 15 quadrature) to +# within the documented Maxwellian-chain gap -- compared through the +# generic reference-data interface, not against PyRPOD's own output. +# +# The study writes into a temporary directory, so a test run leaves no +# artifacts in the repository. +# +# Run: python -m pytest mdao/mdao_integration_test_02.py -s (from tests/) + +import json +import os +import shutil +import tempfile +import unittest +from pathlib import Path + +import numpy as np + +from pyrpod.mdao.TradeStudy import TradeStudy +from pyrpod.mdao.reference_data import ReferenceDataset, ReferenceRecord +from pyrpod.plume import CaiImpingement2016 as cai + +_TESTS_DIR = Path(__file__).resolve().parents[1] +CASE_DIR = _TESTS_DIR.parent / 'case' / 'plume' / 'plume_flat_plate_sweep' +CONFIG_PATH = CASE_DIR / 'study' / 'flat_plate_baseline.yaml' + +# Paper conditions of the case (and of the committed configuration). +S_0 = 2.0 +EPS = 1.5 # Tw / T0 +R_0 = 0.5 # nozzle radius (m) +PLATE_SEMI = 4.0 # 8 m x 8 m plate +PLATE_AREA = 64.0 +Q_DYN = 1.1044652197738332 # n0*m*U0^2/2 (Pa) +Q_DYN_HEAT = 637.3520362956127 # n0*m*U0^3/2 (W/m^2) + + +def cai_reference(alpha_deg, distance): + """Exact Eq.-15 plate averages for one pose, as dimensional loads. + + The sweep angle maps onto the paper's inclination as + alpha_paper = 90 deg - |alpha|. The plate's faces all share the +Z + normal, so the exact pressure average CP becomes a pure normal force + -CP * q_dyn * S, and the average heat-flux coefficient CQ becomes a + component heat load CQ * q_heat * S. + """ + coefficients = cai.averaged_coefficients( + S_0, np.deg2rad(90.0 - abs(alpha_deg)), EPS, R_0, distance, + PLATE_SEMI, PLATE_SEMI) + return { + 'pressure_force': [0.0, 0.0, + -float(coefficients['CP']) * Q_DYN * PLATE_AREA], + 'total_heat_load': float(coefficients['CQ']) * Q_DYN_HEAT * PLATE_AREA, + } + + +class FlatPlateBaselineStudy(unittest.TestCase): + + @classmethod + def setUpClass(cls): + cls.output_dir = tempfile.mkdtemp(prefix='pyrpod_baseline_study_') + cls.study = TradeStudy.from_config(CONFIG_PATH, + output_dir=cls.output_dir) + cls.results = cls.study.run() + + @classmethod + def tearDownClass(cls): + shutil.rmtree(cls.output_dir, ignore_errors=True) + + # ------------------------------------------------------------ structure + def test_one_case_one_component_one_firing(self): + self.assertEqual(len(self.results), 1) + case = self.results.cases[0] + self.assertEqual(case.study_name, 'cai2016_flat_plate_baseline') + self.assertEqual(case.component, 'plate') + self.assertEqual(case.firing_id, 1) + self.assertEqual(case.plate_angle_deg, 0.0) + self.assertEqual(case.source_distance, 4.0) + self.assertEqual(case.plume_model, 'SimplifiedGasKinetics') + self.assertEqual(case.mesh_faces, case.component_faces) + + def test_jfh_holds_exactly_the_requested_firing_count(self): + case = self.results.cases[0] + self.assertTrue(os.path.isfile(case.jfh_path)) + lines = [line for line in + Path(case.jfh_path).read_text(encoding='utf-8').splitlines() + if line.strip()] + # header + unused second line + one firing row + self.assertEqual(len(lines), 3) + self.assertIn('1', lines[0].split()) + + def test_source_pose_is_head_on_at_the_configured_distance(self): + case = self.results.cases[0] + np.testing.assert_allclose(case.plume_source_position, + [0.0, 0.0, 4.0], atol=1e-9) + # DCM first column is the thruster axis: aimed back at the plate. + dcm = np.asarray(case.plume_source_orientation).reshape(3, 3) + np.testing.assert_allclose(dcm[:, 0], [0.0, 0.0, -1.0], atol=1e-9) + + # ------------------------------------------------------------ artifacts + def test_vtk_is_written_with_the_standard_per_face_fields(self): + case = self.results.cases[0] + self.assertIsNotNone(case.vtk_path) + self.assertTrue(os.path.isfile(case.vtk_path), + msg=f'missing VTK artifact {case.vtk_path}') + + head = Path(case.vtk_path).read_bytes()[:8000].decode('utf-8', + 'replace') + for field in ('strikes', 'cum_strikes', 'pressures', 'max_pressures', + 'shear_stress', 'max_shears', 'heat_flux_rate', + 'heat_flux_load', 'cum_heat_flux_load'): + self.assertIn(f'Name="{field}"', head, + msg=f'per-face field {field} missing from the VTK') + + def test_summary_csv_and_metadata_are_written(self): + csv_path = self.results.summary_csv_path + metadata_path = self.results.metadata_path + self.assertTrue(os.path.isfile(csv_path)) + self.assertTrue(os.path.isfile(metadata_path)) + + rows = Path(csv_path).read_text(encoding='utf-8').splitlines() + self.assertEqual(len(rows), 2) # header + one case + for column in ('force_z', 'moment_y', 'center_of_pressure_x', + 'max_pressure', 'affected_area', 'coeff_CF', + 'vtk_path'): + self.assertIn(column, rows[0]) + + document = json.loads(Path(metadata_path).read_text(encoding='utf-8')) + self.assertEqual(document['n_cases'], 1) + provenance = document['provenance'] + self.assertEqual(provenance['plume_model'], 'SimplifiedGasKinetics') + self.assertEqual(provenance['geometry_id'], + 'flat_plate_transformed.stl') + self.assertEqual(provenance['n_firings_per_pose'], 1) + self.assertEqual(provenance['total_firings'], 1) + self.assertEqual(provenance['sweep_mode'], 'per_case') + self.assertEqual(provenance['config_path'], str(CONFIG_PATH)) + self.assertIn('code_version', provenance) + self.assertIn('known_limitations', provenance) + self.assertIn('units', document['cases'][0]) + + # -------------------------------------------------------------- physics + def test_head_on_loads_are_physically_consistent(self): + case = self.results.cases[0] + force = np.asarray(case.force) + + # The plume travels along -Z here, so the plate is pushed along -Z. + self.assertLess(force[2], 0.0) + self.assertLess(abs(force[0]), 1e-6 * abs(force[2])) + self.assertLess(abs(force[1]), 1e-6 * abs(force[2])) + # Symmetric pose: no resultant moment about the plate center. + self.assertLess(case.moment_magnitude, 1e-6 * abs(force[2])) + + self.assertEqual(case.center_of_pressure_status, 'ok') + np.testing.assert_allclose(case.center_of_pressure, + [0.0, 0.0, 0.0], atol=1e-6) + + # The whole plate sits inside the gating wedge at this pose. + self.assertEqual(case.struck_faces, case.component_faces) + self.assertAlmostEqual(case.affected_area, PLATE_AREA, places=6) + self.assertGreater(case.max_pressure, 0.0) + self.assertGreater(case.max_heat_flux, 0.0) + + def test_coefficients_are_available_with_the_configured_normalization(self): + case = self.results.cases[0] + self.assertTrue(case.coefficients_available) + for name in ('CF', 'CFz', 'CM', 'Cp_max', 'Cf_max', 'Cq_max'): + self.assertIn(name, case.coefficients) + self.assertAlmostEqual( + case.coefficients['CF'], + case.force_magnitude / (Q_DYN * PLATE_AREA), places=12) + + # ------------------------------------------- independent Cai reference + def test_matches_the_independent_cai_reference(self): + dataset = ReferenceDataset( + label='Cai 2016 exact (Eq. 15)', + source='pyrpod/plume/CaiImpingement2016.py', + records=[ReferenceRecord( + key={'plate_angle_deg': 0.0, 'source_distance': 4.0}, + quantities=cai_reference(0.0, 4.0))]) + + report = self.study.compare(dataset) + self.assertEqual(report.unmatched_cases, []) + by_name = {entry.quantity: entry for entry in report.comparisons} + self.assertEqual(set(by_name), {'pressure_force', 'total_heat_load'}) + + force_error = by_name['pressure_force'].relative_error + heat_error = by_name['total_heat_load'].relative_error + print(f'[baseline] integrated normal load vs Cai 2016 exact: ' + f'{force_error:.3%} relative error') + print(f'[baseline] component heat load vs Cai 2016 exact: ' + f'{heat_error:.3%} relative error') + + # Documented accuracy of the Maxwellian engineering chain against the + # exact collisionless solution (plate-averaged CP tracks the + # reference to a few percent, the heat flux to ~15%). A convention + # or sign regression would be O(1) here. + self.assertLess(force_error, 0.08) + self.assertLess(heat_error, 0.25) + + report_path = Path(self.output_dir) / 'cai_reference_comparison.csv' + report.write_csv(report_path) + self.assertTrue(report_path.is_file()) + + +if __name__ == '__main__': + unittest.main() diff --git a/tests/mdao/mdao_integration_test_03.py b/tests/mdao/mdao_integration_test_03.py new file mode 100644 index 0000000..36f1588 --- /dev/null +++ b/tests/mdao/mdao_integration_test_03.py @@ -0,0 +1,316 @@ +# ======================== +# PyRPOD: tests/mdao/mdao_integration_test_03.py +# ======================== +# Multi-angle / multi-distance flat-plate sweep through the package-level +# trade study API, plus the architecture checks that keep the study honest +# for non-plate geometry. +# +# The committed sweep configuration +# (case/plume/plume_flat_plate_sweep/study/flat_plate_sweep.yaml) covers 19 +# angles x 5 distances; this test runs a REDUCED subset of the same +# configuration (5 angles x 2 distances = 10 cases) so the automated suite +# stays quick -- the exact-reference quadrature, not the pipeline, dominates +# the runtime. What is checked: +# +# * one structured result per angle x distance combination, each with +# exactly the requested number of JFH entries; +# * mirror symmetry in +/- angle, which a paper-frame convention error +# would break as O(1); +# * monotonic decay of the integrated load with source distance, and +# center-of-pressure travel with approach angle; +# * agreement with the INDEPENDENT Cai 2016 exact reference across the +# sweep, through the generic reference-data interface; +# * per-case VTK artifacts land in per-case directories (sweep cases never +# overwrite one another); +# * the optional trend plots can be generated (into a temporary directory, +# so no artifact reaches the repository); +# * the same machinery runs against the CYLINDER case -- nothing in the +# study assumes a flat target -- and correctly reports coefficients as +# unavailable when the configuration supplies no normalization inputs. +# +# Run: python -m pytest mdao/mdao_integration_test_03.py -s (from tests/) + +import os +import shutil +import tempfile +import unittest +from pathlib import Path + +import numpy as np +import yaml + +from pyrpod.mdao.TradeStudy import TradeStudy +from pyrpod.mdao.plume_validation import PlumeValidationStudy +from pyrpod.mdao.reference_data import ( + ReferenceDataset, + ReferenceRecord, + compare_results, +) +from pyrpod.mdao.study_config import StudyConfig +from pyrpod.plume import CaiImpingement2016 as cai + +_TESTS_DIR = Path(__file__).resolve().parents[1] +CASE_DIR = _TESTS_DIR.parent / 'case' / 'plume' / 'plume_flat_plate_sweep' +SWEEP_CONFIG = CASE_DIR / 'study' / 'flat_plate_sweep.yaml' +CYLINDER_CONFIG = (_TESTS_DIR.parent / 'case' / 'plume' + / 'plume_cylinder_sweep' / 'study' + / 'cylinder_baseline.yaml') + +ANGLES = [-40.0, -20.0, 0.0, 20.0, 40.0] +DISTANCES = [2.0, 6.0] + +S_0 = 2.0 +EPS = 1.5 +R_0 = 0.5 +PLATE_SEMI = 4.0 +PLATE_AREA = 64.0 +Q_DYN = 1.1044652197738332 +Q_DYN_HEAT = 637.3520362956127 + + +def reduced_sweep_config(output_dir, angles=ANGLES, distances=DISTANCES, + write_vtk=False, write_plots=False): + """The committed sweep configuration, restricted to a small subset.""" + data = yaml.safe_load(SWEEP_CONFIG.read_text(encoding='utf-8')) + data['sweep']['plate_angles_deg'] = list(angles) + data['sweep']['source_distances'] = list(distances) + data['output']['vtk']['enabled'] = write_vtk + data['output']['plots']['enabled'] = write_plots + data['study']['output_dir'] = str(output_dir) + return StudyConfig.from_mapping(data, source_path=str(SWEEP_CONFIG)) + + +def cai_reference_dataset(angles, distances): + """Exact Cai 2016 plate averages as dimensional reference records. + + Independent of PyRPOD's strike pipeline: the Eq.-15 quadrature of + pyrpod/plume/CaiImpingement2016.py evaluated at alpha_paper = 90 - |alpha|, + converted to a normal pressure force and a component heat load with the + case's own normalization. Mirror-invariant in +/- alpha, so it is cached + per |alpha|. + """ + cache = {} + records = [] + for distance in distances: + for angle in angles: + key = (abs(angle), distance) + if key not in cache: + cache[key] = cai.averaged_coefficients( + S_0, np.deg2rad(90.0 - abs(angle)), EPS, R_0, distance, + PLATE_SEMI, PLATE_SEMI) + coefficients = cache[key] + records.append(ReferenceRecord( + key={'plate_angle_deg': float(angle), + 'source_distance': float(distance), + 'component': 'plate'}, + quantities={ + 'pressure_force': [ + 0.0, 0.0, + -float(coefficients['CP']) * Q_DYN * PLATE_AREA], + 'total_heat_load': (float(coefficients['CQ']) + * Q_DYN_HEAT * PLATE_AREA)})) + return ReferenceDataset(label='Cai 2016 exact (Eq. 15)', + source='pyrpod/plume/CaiImpingement2016.py', + records=records) + + +class FlatPlateSweep(unittest.TestCase): + + @classmethod + def setUpClass(cls): + cls.output_dir = tempfile.mkdtemp(prefix='pyrpod_sweep_study_') + cls.config = reduced_sweep_config(cls.output_dir) + cls.study = PlumeValidationStudy(cls.config) + cls.results = cls.study.run() + cls.by_pose = {(case.plate_angle_deg, case.source_distance): case + for case in cls.results.cases} + + @classmethod + def tearDownClass(cls): + shutil.rmtree(cls.output_dir, ignore_errors=True) + + def test_every_angle_distance_combination_produced_a_result(self): + self.assertEqual(len(self.results), len(ANGLES) * len(DISTANCES)) + self.assertEqual(set(self.by_pose), + {(angle, distance) for angle in ANGLES + for distance in DISTANCES}) + for case in self.results.cases: + self.assertEqual(case.component, 'plate') + self.assertEqual(case.firing_id, 1) + self.assertEqual(case.study_name, 'cai2016_flat_plate_sweep') + + def test_each_case_wrote_its_own_jfh_with_the_exact_entry_count(self): + seen = set() + for case in self.results.cases: + self.assertTrue(os.path.isfile(case.jfh_path)) + self.assertNotIn(case.jfh_path, seen) + seen.add(case.jfh_path) + header = Path(case.jfh_path).read_text( + encoding='utf-8').splitlines()[0] + self.assertEqual(int(header.split()[1]), + self.config.sweep.n_firings) + + def test_results_are_mirror_symmetric_in_plus_minus_angle(self): + for distance in DISTANCES: + for angle in (20.0, 40.0): + positive = self.by_pose[(angle, distance)] + negative = self.by_pose[(-angle, distance)] + scale = abs(positive.force[2]) + + # Normal load and heat load are mirror invariant. + self.assertAlmostEqual(positive.force[2], negative.force[2], + delta=1e-6 * scale) + self.assertAlmostEqual(positive.total_heat_load, + negative.total_heat_load, + delta=1e-6 * positive.total_heat_load) + # The in-plane resultant mirrors. + self.assertAlmostEqual(positive.force[0], -negative.force[0], + delta=1e-6 * scale) + # ... and so does the center of pressure. + self.assertAlmostEqual(positive.center_of_pressure[0], + -negative.center_of_pressure[0], + delta=1e-6 * PLATE_SEMI) + + def test_out_of_plane_resultant_vanishes_by_symmetry(self): + for case in self.results.cases: + self.assertLess(abs(case.force[1]), 1e-6 * abs(case.force[2])) + + def test_load_decays_with_source_distance(self): + for angle in ANGLES: + near = self.by_pose[(angle, min(DISTANCES))] + far = self.by_pose[(angle, max(DISTANCES))] + self.assertLess(abs(far.force[2]), abs(near.force[2])) + self.assertLess(far.max_pressure, near.max_pressure) + self.assertLess(far.max_heat_flux, near.max_heat_flux) + + def test_center_of_pressure_moves_with_approach_angle(self): + for distance in DISTANCES: + head_on = self.by_pose[(0.0, distance)] + inclined = self.by_pose[(-40.0, distance)] + self.assertEqual(head_on.center_of_pressure_status, 'ok') + self.assertEqual(inclined.center_of_pressure_status, 'ok') + np.testing.assert_allclose(head_on.center_of_pressure[:2], + [0.0, 0.0], atol=1e-6) + self.assertGreater( + abs(inclined.center_of_pressure[0] + - head_on.center_of_pressure[0]), 1e-3) + + def test_sweep_agrees_with_the_independent_cai_reference(self): + dataset = cai_reference_dataset(ANGLES, DISTANCES) + report = compare_results(self.results.cases, dataset, + reference_length=PLATE_SEMI) + self.assertEqual(report.unmatched_cases, []) + self.assertEqual(len(report.comparisons), + 2 * len(ANGLES) * len(DISTANCES)) + + for quantity, envelope in (('pressure_force', 0.08), + ('total_heat_load', 0.25)): + errors = [entry.relative_error + for entry in report.for_quantity(quantity)] + worst, mean = max(errors), sum(errors) / len(errors) + print(f'[sweep] {quantity} vs Cai 2016 exact: max ' + f'{worst:.3%}, mean {mean:.3%} over {len(errors)} poses') + # Documented Maxwellian-chain gap, not a physics tolerance: a + # convention regression would be O(1) rather than a few percent. + self.assertLess(mean, envelope) + + report_path = Path(self.output_dir) / 'cai_reference_comparison.csv' + report.write_csv(report_path) + self.assertTrue(report_path.is_file()) + + def test_summary_artifacts_cover_every_case(self): + rows = Path(self.results.summary_csv_path).read_text( + encoding='utf-8').splitlines() + self.assertEqual(len(rows), len(self.results) + 1) + self.assertTrue(os.path.isfile(self.results.metadata_path)) + + +class SweepArtifacts(unittest.TestCase): + """VTK export and optional plots, on a two-case slice of the sweep.""" + + @classmethod + def setUpClass(cls): + cls.output_dir = tempfile.mkdtemp(prefix='pyrpod_sweep_artifacts_') + cls.config = reduced_sweep_config(cls.output_dir, + angles=[-30.0, 30.0], + distances=[4.0], write_vtk=True) + cls.study = PlumeValidationStudy(cls.config) + cls.results = cls.study.run() + + @classmethod + def tearDownClass(cls): + shutil.rmtree(cls.output_dir, ignore_errors=True) + + def test_each_case_writes_its_own_vtk(self): + paths = [case.vtk_path for case in self.results.cases] + self.assertEqual(len(paths), 2) + self.assertEqual(len(set(paths)), 2, + msg='sweep cases must not share a VTK path') + for path in paths: + self.assertTrue(os.path.isfile(path)) + self.assertTrue(Path(path).name.startswith('firing-')) + + def test_optional_plots_are_generated_on_request(self): + # Plot generation is optional: nothing above needed it. Here it is + # exercised explicitly, writing into the temporary output directory. + paths = self.study.plot() + self.assertTrue(paths) + for path in paths: + self.assertTrue(os.path.isfile(path)) + names = {Path(path).name for path in paths} + self.assertIn('force_vs_angle.png', names) + self.assertIn('moment_vs_angle.png', names) + self.assertIn('heat_flux_vs_angle.png', names) + + +class CylinderTargetSupport(unittest.TestCase): + """The study architecture does not assume a flat plate.""" + + @classmethod + def setUpClass(cls): + cls.output_dir = tempfile.mkdtemp(prefix='pyrpod_cylinder_study_') + cls.study = TradeStudy.from_config(CYLINDER_CONFIG, + output_dir=cls.output_dir) + cls.results = cls.study.run() + + @classmethod + def tearDownClass(cls): + shutil.rmtree(cls.output_dir, ignore_errors=True) + + def test_curved_closed_target_runs_end_to_end(self): + self.assertEqual(len(self.results), 1) + case = self.results.cases[0] + self.assertEqual(case.component, 'cylinder') + self.assertGreater(case.mesh_faces, 0) + self.assertGreater(case.struck_faces, 0) + self.assertLess(case.struck_faces, case.mesh_faces, + msg='a closed target must not be struck all over ' + 'from a single pose') + self.assertGreater(case.affected_area, 0.0) + self.assertLess(case.affected_area, case.component_area) + self.assertTrue(os.path.isfile(case.vtk_path)) + + def test_coefficients_are_unavailable_without_normalization_inputs(self): + case = self.results.cases[0] + self.assertFalse(case.coefficients_available) + self.assertEqual(case.coefficients, {}) + row = case.to_row() + self.assertNotIn('coeff_CF', row) + + def test_center_of_pressure_is_reported_with_a_status(self): + case = self.results.cases[0] + self.assertIn(case.center_of_pressure_status, + ('ok', 'ill_conditioned', 'zero_load')) + if case.center_of_pressure_status == 'ok': + # Whatever it is, it must reproduce the reported moment. + arm = (np.asarray(case.center_of_pressure) + - np.asarray(case.moment_reference_point)) + np.testing.assert_allclose( + np.cross(arm, case.force), case.moment, + atol=1e-6 * max(1.0, case.moment_magnitude)) + else: + self.assertIsNone(case.center_of_pressure) + + +if __name__ == '__main__': + unittest.main() diff --git a/tests/mdao/mdao_integration_test_04.py b/tests/mdao/mdao_integration_test_04.py new file mode 100644 index 0000000..6c8fa4e --- /dev/null +++ b/tests/mdao/mdao_integration_test_04.py @@ -0,0 +1,255 @@ +# ======================== +# PyRPOD: tests/mdao/mdao_integration_test_04.py +# ======================== +# The single-history sweep engine (pyrpod.mdao.parameter_sweep. +# ParameterSweepStudy, sweep.mode: single_jfh): the whole angle x distance +# sweep runs as ONE case driven by ONE Jet Firing History, so every firing's +# strikes come from the same pipeline run. +# +# The committed configuration +# (case/plume/plume_flat_plate_sweep/study/flat_plate_sweep_single_jfh.yaml) +# covers 19 angles x 5 distances; this test runs a REDUCED subset of it so +# the automated suite stays quick. What is checked: +# +# * TradeStudy.from_config dispatches on sweep.mode -- single_jfh binds the +# sweep engine, per_case the original one; +# * exactly ONE Jet Firing History is written, holding exactly +# len(poses) x n_firings entries, and one result record per firing; +# * every record shares one case_id but carries the pose it realizes, with +# firings numbered continuously through the history; +# * EQUIVALENCE: per-firing loads are identical to the per-case engine's, +# so the two decompositions are numerically interchangeable; +# * the strike VTK files form a single results/strikes/firing-.vtu +# series (the repository's sweep convention), not one folder per case; +# * the sweep envelope -- worst pressure/shear over every pose, accumulated +# heat-flux load, coverage -- is recorded per component, and is absent +# (not fabricated) when the full pipeline path is disabled; +# * the same independent Cai 2016 reference data compares against this +# engine's results unchanged. +# +# Run: python -m pytest mdao/mdao_integration_test_04.py -s (from tests/) + +import os +import shutil +import tempfile +import unittest +from pathlib import Path + +import numpy as np +import yaml + +from pyrpod.mdao.TradeStudy import TradeStudy +from pyrpod.mdao.parameter_sweep import ParameterSweepStudy +from pyrpod.mdao.plume_validation import PlumeValidationStudy +from pyrpod.mdao.reference_data import compare_results +from pyrpod.mdao.study_config import StudyConfig + +from mdao_integration_test_03 import ( # noqa: E402 (sibling test module) + ANGLES, + DISTANCES, + PLATE_SEMI, + cai_reference_dataset, +) + +_TESTS_DIR = Path(__file__).resolve().parents[1] +CASE_DIR = _TESTS_DIR.parent / 'case' / 'plume' / 'plume_flat_plate_sweep' +SINGLE_JFH_CONFIG = CASE_DIR / 'study' / 'flat_plate_sweep_single_jfh.yaml' +PER_CASE_CONFIG = CASE_DIR / 'study' / 'flat_plate_sweep.yaml' + + +def reduced_config(source, output_dir, angles=ANGLES, distances=DISTANCES, + write_vtk=False, n_firings=1): + """A committed configuration restricted to a small subset of poses.""" + data = yaml.safe_load(Path(source).read_text(encoding='utf-8')) + data['sweep']['plate_angles_deg'] = list(angles) + data['sweep']['source_distances'] = list(distances) + data['sweep']['n_firings'] = n_firings + data['output']['vtk']['enabled'] = write_vtk + data['output']['plots']['enabled'] = False + data['study']['output_dir'] = str(output_dir) + return StudyConfig.from_mapping(data, source_path=str(source)) + + +class EngineSelection(unittest.TestCase): + + def test_sweep_mode_selects_the_engine(self): + with tempfile.TemporaryDirectory() as tmp: + single = TradeStudy.from_config(SINGLE_JFH_CONFIG, output_dir=tmp) + per_case = TradeStudy.from_config(PER_CASE_CONFIG, output_dir=tmp) + self.assertIsInstance(single.validation_study, ParameterSweepStudy) + self.assertIsInstance(per_case.validation_study, PlumeValidationStudy) + self.assertEqual(single.study_config.sweep.mode, 'single_jfh') + self.assertEqual(per_case.study_config.sweep.mode, 'per_case') + + +class SingleHistorySweep(unittest.TestCase): + + @classmethod + def setUpClass(cls): + cls.output_dir = tempfile.mkdtemp(prefix='pyrpod_single_jfh_') + cls.config = reduced_config(SINGLE_JFH_CONFIG, cls.output_dir) + cls.study = ParameterSweepStudy(cls.config) + cls.results = cls.study.run() + + @classmethod + def tearDownClass(cls): + shutil.rmtree(cls.output_dir, ignore_errors=True) + + def test_one_history_holds_every_pose(self): + expected = len(ANGLES) * len(DISTANCES) + self.assertEqual(self.config.sweep.total_firings, expected) + self.assertEqual(len(self.results), expected) + + jfh_paths = {case.jfh_path for case in self.results.cases} + self.assertEqual(len(jfh_paths), 1, msg='the sweep must write ONE JFH') + jfh_path = jfh_paths.pop() + self.assertTrue(os.path.isfile(jfh_path)) + + lines = [line for line in + Path(jfh_path).read_text(encoding='utf-8').splitlines() + if line.strip()] + # header + unused second line + one row per firing + self.assertEqual(len(lines), expected + 2) + self.assertEqual(int(lines[0].split()[1]), expected) + + def test_records_share_a_case_but_carry_their_own_pose(self): + case_ids = {case.case_id for case in self.results.cases} + self.assertEqual(case_ids, {self.study.case_id}) + self.assertEqual([case.firing_id for case in self.results.cases], + list(range(1, len(self.results) + 1))) + self.assertEqual( + [(case.plate_angle_deg, case.source_distance) + for case in self.results.cases], + [pose for pose in self.config.sweep.poses]) + + def test_multiple_firings_per_pose_extend_the_same_history(self): + with tempfile.TemporaryDirectory() as tmp: + config = reduced_config(SINGLE_JFH_CONFIG, tmp, + angles=[-20.0, 20.0], distances=[4.0], + n_firings=3) + results = ParameterSweepStudy(config).run() + + self.assertEqual(len(results), 6) # 2 poses x 3 firings + self.assertEqual([case.plate_angle_deg for case in results.cases], + [-20.0] * 3 + [20.0] * 3) + self.assertEqual({case.jfh_path for case in results.cases}.__len__(), 1) + + def test_per_firing_loads_match_the_per_case_engine(self): + with tempfile.TemporaryDirectory() as tmp: + per_case = PlumeValidationStudy( + reduced_config(PER_CASE_CONFIG, tmp)).run() + + by_pose = {(case.plate_angle_deg, case.source_distance): case + for case in per_case.cases} + self.assertEqual(len(by_pose), len(self.results)) + + for case in self.results.cases: + reference = by_pose[(case.plate_angle_deg, case.source_distance)] + np.testing.assert_array_equal(case.force, reference.force) + np.testing.assert_array_equal(case.moment, reference.moment) + self.assertEqual(case.max_pressure, reference.max_pressure) + self.assertEqual(case.max_heat_flux, reference.max_heat_flux) + self.assertEqual(case.struck_faces, reference.struck_faces) + + def test_envelope_is_absent_without_the_full_pipeline(self): + # This run has VTK disabled, so the pipeline's cumulative arrays were + # never produced; the envelope must be empty rather than invented. + self.assertEqual(self.study.envelope, {}) + self.assertNotIn('sweep_envelope', self.results.provenance) + + def test_metadata_records_the_sweep_mode_and_history(self): + provenance = self.results.provenance + self.assertEqual(provenance['sweep_mode'], 'single_jfh') + self.assertEqual(provenance['n_firings_per_pose'], 1) + self.assertEqual(provenance['total_firings'], + len(ANGLES) * len(DISTANCES)) + self.assertTrue(os.path.isfile(provenance['jfh_path'])) + + def test_independent_cai_reference_compares_unchanged(self): + dataset = cai_reference_dataset(ANGLES, DISTANCES) + report = compare_results(self.results.cases, dataset, + reference_length=PLATE_SEMI) + self.assertEqual(report.unmatched_cases, []) + + errors = [entry.relative_error + for entry in report.for_quantity('pressure_force')] + worst, mean = max(errors), sum(errors) / len(errors) + print(f'[single-jfh] pressure_force vs Cai 2016 exact: max ' + f'{worst:.3%}, mean {mean:.3%} over {len(errors)} poses') + self.assertLess(mean, 0.08) + + +class SingleHistoryArtifacts(unittest.TestCase): + """VTK layout and the sweep envelope, on a three-pose slice.""" + + ANGLES = [-30.0, 0.0, 30.0] + + @classmethod + def setUpClass(cls): + cls.output_dir = tempfile.mkdtemp(prefix='pyrpod_single_jfh_vtk_') + cls.config = reduced_config(SINGLE_JFH_CONFIG, cls.output_dir, + angles=cls.ANGLES, distances=[4.0], + write_vtk=True) + cls.study = ParameterSweepStudy(cls.config) + cls.results = cls.study.run() + + @classmethod + def tearDownClass(cls): + shutil.rmtree(cls.output_dir, ignore_errors=True) + + def test_strikes_form_one_numbered_series(self): + strikes_dir = Path(self.output_dir) / 'results' / 'strikes' + self.assertTrue(strikes_dir.is_dir()) + self.assertEqual(sorted(p.name for p in strikes_dir.glob('*.vtu')), + [f'firing-{i}.vtu' for i in range(len(self.ANGLES))]) + # ... and every record points into that one series. + for index, case in enumerate(self.results.cases): + self.assertEqual(Path(case.vtk_path).name, f'firing-{index}.vtu') + self.assertTrue(os.path.isfile(case.vtk_path)) + self.assertFalse((Path(self.output_dir) / 'cases').exists(), + msg='single_jfh mode must not create per-case dirs') + + def test_cumulative_fields_are_in_the_last_firing_vtk(self): + last = Path(self.results.cases[-1].vtk_path) + head = last.read_bytes()[:8000].decode('utf-8', 'replace') + for field in ('strikes', 'cum_strikes', 'pressures', 'max_pressures', + 'shear_stress', 'max_shears', 'heat_flux_rate', + 'cum_heat_flux_load'): + self.assertIn(f'Name="{field}"', head) + + def test_sweep_envelope_bounds_every_pose(self): + envelope = self.study.envelope + self.assertIn('plate', envelope) + plate = envelope['plate'] + + peak_pressure = max(case.max_pressure for case in self.results.cases) + peak_shear = max(case.max_shear_stress + for case in self.results.cases) + self.assertAlmostEqual(plate['max_pressure'], peak_pressure, places=9) + self.assertAlmostEqual(plate['max_shear_stress'], peak_shear, + places=9) + + # Coverage: each pose struck at most the whole plate, and the sweep + # as a whole struck at least as much as any single pose did. + self.assertGreaterEqual( + plate['unique_struck_faces'], + max(case.struck_faces for case in self.results.cases)) + self.assertGreaterEqual( + plate['swept_affected_area'], + max(case.affected_area for case in self.results.cases)) + self.assertAlmostEqual(plate['component_area'], + self.results.cases[0].component_area, + places=9) + self.assertEqual(self.results.provenance['sweep_envelope'], envelope) + + def test_optional_plots_work_from_the_single_history_results(self): + paths = self.study.plot() + self.assertTrue(paths) + names = {Path(path).name for path in paths} + self.assertIn('force_vs_angle.png', names) + for path in paths: + self.assertTrue(os.path.isfile(path)) + + +if __name__ == '__main__': + unittest.main() diff --git a/tests/mdao/mdao_unit_test_03.py b/tests/mdao/mdao_unit_test_03.py new file mode 100644 index 0000000..fd4b2a4 --- /dev/null +++ b/tests/mdao/mdao_unit_test_03.py @@ -0,0 +1,350 @@ +# ======================== +# PyRPOD: tests/mdao/mdao_unit_test_03.py +# ======================== +# Unit tests for prescribed firing generation and the exact meaning of +# n_firings (pyrpod.mdao.firing_plan, pyrpod.rpod.approach_maneuvers): +# +# * invalid firing counts (zero, negative, fractional, non-numeric) are +# rejected rather than coerced; +# * one requested firing produces exactly one JFH entry, and N requested +# firings produce exactly N -- verified by reading the written file back +# through the normal JetFiringHistory parser; +# * an explicit firing list whose length disagrees with n_firings is an +# error, never a silent truncation; +# * the generated pose convention reproduces the committed sweep-JFH +# generator of case/plume/plume_flat_plate_sweep to file precision; +# * compute_1d_approach honors an exact n_firings on the dynamics path. +# +# No case file I/O beyond reading the committed sweep JFH; generated JFH +# files are written to a temporary directory. +# +# Run: python -m pytest mdao/mdao_unit_test_03.py (from tests/) + +import tempfile +import unittest +from pathlib import Path + +import numpy as np +import pytest + +from pyrpod.mdao.firing_plan import ( + build_case_firings, + build_sweep_firings, + pose_for, + validate_n_firings, + write_jfh_file, +) +from pyrpod.mdao.study_config import StudyConfigError, SweepSpec, TargetSpec +from pyrpod.rpod import JetFiringHistory +from pyrpod.rpod.approach_maneuvers import ApproachInputs, compute_1d_approach + +CASE_DIR = '../case/plume/plume_flat_plate_sweep/' + +# The committed sweep: 19 approach angles x 5 stand-off distances, plate +# centered at the origin, normal +Z, tangent +X (see the case's +# jfh/generate_sweep_jfh.py, run with --alpha0-deg 0 --distance 0). +SWEEP_ANGLES = np.arange(-90.0, 90.0 + 1e-9, 10.0) +SWEEP_DISTANCES = [2.0, 4.0, 6.0, 8.0, 10.0] + + +def _target(reference_point=(0.0, 0.0, 0.0)): + return TargetSpec.from_mapping( + {'reference_point': list(reference_point), + 'normal': [0.0, 0.0, 1.0], + 'tangent': [1.0, 0.0, 0.0]}, + default_geometry_id='test_geometry.stl') + + +def _sweep(n_firings=1, firings=None, angles=(0.0,), distances=(4.0,), + mode='per_case'): + data = {'plate_angles_deg': list(angles), + 'source_distances': list(distances), + 'n_firings': n_firings, 'firing_duration_s': 1.0, + 'thrusters': [1], 'mode': mode} + if firings is not None: + data['firings'] = firings + return SweepSpec.from_mapping(data) + + +def _read_back(firings): + """Write firings to a temp JFH and parse it with JetFiringHistory.""" + with tempfile.TemporaryDirectory() as tmp: + path = Path(tmp) / 'jfh' / 'generated.A' + written = write_jfh_file(path, firings) + + jfh = JetFiringHistory.JetFiringHistory(CASE_DIR) + jfh.case_dir = str(tmp) + '/' + jfh.config.set('jfh', 'jfh', 'generated.A') + jfh.read_jfh() + return written, list(jfh.JFH) + + +class FiringCountValidation(unittest.TestCase): + """n_firings is an exact entry count, so only positive integers pass.""" + + def test_invalid_counts_are_rejected(self): + for value in (0, -1, -10, 1.5, 0.0, '2', None, True): + with self.subTest(value=value): + with pytest.raises(ValueError): + validate_n_firings(value) + + def test_valid_counts_are_returned_as_int(self): + for value in (1, 3, 95, 10.0): + with self.subTest(value=value): + self.assertEqual(validate_n_firings(value), int(value)) + + def test_zero_firings_rejected_by_the_sweep_specification(self): + with pytest.raises(StudyConfigError): + _sweep(n_firings=0) + + def test_writing_an_empty_jfh_is_refused(self): + with tempfile.TemporaryDirectory() as tmp: + with pytest.raises(StudyConfigError): + write_jfh_file(Path(tmp) / 'empty.A', []) + + +class GeneratedFiringCounts(unittest.TestCase): + """N requested firings produce exactly N JFH entries.""" + + def test_one_firing_produces_one_entry(self): + firings = build_case_firings(_sweep(n_firings=1), _target(), 0.0, 4.0) + self.assertEqual(len(firings), 1) + + written, entries = _read_back(firings) + self.assertEqual(written, 1) + self.assertEqual(len(entries), 1) + self.assertEqual(entries[0]['thrusters'], [1]) + + def test_n_firings_produce_exactly_n_entries(self): + for n_firings in (1, 2, 5, 17): + with self.subTest(n_firings=n_firings): + firings = build_case_firings(_sweep(n_firings=n_firings), + _target(), -30.0, 6.0) + self.assertEqual(len(firings), n_firings) + + written, entries = _read_back(firings) + self.assertEqual(written, n_firings) + self.assertEqual(len(entries), n_firings) + # Every entry is a complete, valid firing record. + for index, entry in enumerate(entries): + self.assertEqual(int(entry['nt']), index + 1) + self.assertEqual(np.asarray(entry['dcm']).shape, (3, 3)) + self.assertEqual(len(entry['xyz']), 3) + self.assertGreater(float(entry['t']), 0.0) + + def test_repeated_firings_share_the_pose_and_advance_in_time(self): + firings = build_case_firings(_sweep(n_firings=3), _target(), 20.0, 8.0) + for firing in firings[1:]: + np.testing.assert_allclose(firing.position, firings[0].position) + np.testing.assert_allclose(firing.dcm, firings[0].dcm) + self.assertEqual([f.start_time_s for f in firings], [0.0, 1.0, 2.0]) + + +class PrescribedFiringList(unittest.TestCase): + """Explicit firing sequences are honored, mismatches are not.""" + + def test_explicit_firings_are_used_verbatim(self): + explicit = [{'position': [0.0, 0.0, 3.0], + 'dcm': [[0.0, 0.0, 1.0], [0.0, 1.0, 0.0], + [-1.0, 0.0, 0.0]], + 'thrusters': [1], 'duration_s': 0.5}, + {'position': [1.0, 0.0, 3.0], + 'dcm': [[0.0, 0.0, 1.0], [0.0, 1.0, 0.0], + [-1.0, 0.0, 0.0]], + 'thrusters': [1], 'duration_s': 0.5}] + sweep = _sweep(n_firings=2, firings=explicit) + firings = build_case_firings(sweep, _target(), 0.0, 4.0) + + self.assertEqual(len(firings), 2) + np.testing.assert_allclose(firings[0].position, [0.0, 0.0, 3.0]) + np.testing.assert_allclose(firings[1].position, [1.0, 0.0, 3.0]) + self.assertEqual([f.duration_s for f in firings], [0.5, 0.5]) + + written, entries = _read_back(firings) + self.assertEqual((written, len(entries)), (2, 2)) + + def test_count_mismatch_is_an_error(self): + explicit = [{'position': [0.0, 0.0, 3.0], + 'dcm': np.eye(3).tolist()}] + with pytest.raises(StudyConfigError) as excinfo: + _sweep(n_firings=4, firings=explicit) + self.assertIn('n_firings', str(excinfo.value)) + + +class WholeSweepFiringSequence(unittest.TestCase): + """single_jfh mode: one history spanning every pose, exact length.""" + + ANGLES = (-30.0, 0.0, 30.0) + DISTANCES = (2.0, 6.0) + + def _sweep(self, n_firings=1, firings=None): + return _sweep(n_firings=n_firings, firings=firings, + angles=self.ANGLES, distances=self.DISTANCES, + mode='single_jfh') + + def test_pose_order_is_distance_major(self): + sweep = self._sweep() + self.assertEqual( + sweep.poses, + tuple((angle, distance) for distance in self.DISTANCES + for angle in self.ANGLES)) + self.assertEqual(sweep.total_firings, 6) + + def test_sequence_holds_exactly_poses_times_n_firings(self): + for n_firings in (1, 3): + with self.subTest(n_firings=n_firings): + sweep = self._sweep(n_firings=n_firings) + firings = build_sweep_firings(sweep, _target()) + expected = len(self.ANGLES) * len(self.DISTANCES) * n_firings + self.assertEqual(len(firings), expected) + self.assertEqual(sweep.total_firings, expected) + + written, entries = _read_back(firings) + self.assertEqual(written, expected) + self.assertEqual(len(entries), expected) + + def test_each_firing_is_tagged_with_the_pose_it_realizes(self): + firings = build_sweep_firings(self._sweep(n_firings=2), _target()) + tags = [(f.plate_angle_deg, f.source_distance) for f in firings] + # Two consecutive entries per pose, poses in distance-major order. + expected = [pose for pose in self._sweep().poses for _ in range(2)] + self.assertEqual(tags, expected) + self.assertEqual([f.pose_index for f in firings], + [i for i in range(6) for _ in range(2)]) + + def test_firing_times_run_continuously_across_the_sweep(self): + firings = build_sweep_firings(self._sweep(n_firings=2), _target()) + self.assertEqual([f.start_time_s for f in firings], + [float(i) for i in range(len(firings))]) + + def test_poses_match_the_per_case_engine_pose_by_pose(self): + sweep = self._sweep() + combined = build_sweep_firings(sweep, _target()) + for index, (angle, distance) in enumerate(sweep.poses): + per_case = build_case_firings(_sweep(angles=self.ANGLES, + distances=self.DISTANCES), + _target(), angle, distance) + np.testing.assert_allclose(combined[index].position, + per_case[0].position) + np.testing.assert_allclose(combined[index].dcm, per_case[0].dcm) + + def test_explicit_firings_are_sliced_across_the_poses(self): + # Six poses x one firing each: the supplied sequence IS the history. + explicit = [{'position': [0.0, 0.0, float(i + 2)], + 'dcm': np.eye(3).tolist()} for i in range(6)] + sweep = self._sweep(n_firings=1, firings=explicit) + firings = build_sweep_firings(sweep, _target()) + + self.assertEqual(len(firings), 6) + self.assertEqual([f.position[2] for f in firings], + [2.0, 3.0, 4.0, 5.0, 6.0, 7.0]) + # ... and each entry still carries the pose it stands for. + self.assertEqual([(f.plate_angle_deg, f.source_distance) + for f in firings], list(sweep.poses)) + + def test_explicit_firing_count_must_cover_every_pose(self): + explicit = [{'position': [0.0, 0.0, 4.0], 'dcm': np.eye(3).tolist()}] + with pytest.raises(StudyConfigError) as excinfo: + self._sweep(n_firings=1, firings=explicit) + self.assertIn('single_jfh', str(excinfo.value)) + + +class PoseConvention(unittest.TestCase): + """The generated poses match the committed sweep-JFH generator.""" + + def test_head_on_pose_is_on_the_target_normal(self): + position, dcm = pose_for(0.0, 4.0, [0.0, 0.0, 0.0], + [0.0, 0.0, 1.0], [1.0, 0.0, 0.0]) + np.testing.assert_allclose(position, [0.0, 0.0, 4.0], atol=1e-12) + # The DCM's first column is the thruster axis, aimed at the target. + np.testing.assert_allclose(dcm[:, 0], [0.0, 0.0, -1.0], atol=1e-12) + self.assertAlmostEqual(float(np.linalg.det(dcm)), 1.0, places=12) + + def test_positive_angle_rotates_toward_the_tangent(self): + position, _ = pose_for(90.0, 5.0, [0.0, 0.0, 0.0], + [0.0, 0.0, 1.0], [1.0, 0.0, 0.0]) + np.testing.assert_allclose(position, [5.0, 0.0, 0.0], atol=1e-12) + + def test_matches_the_committed_flat_plate_sweep_jfh(self): + jfh = JetFiringHistory.JetFiringHistory(CASE_DIR) + jfh.read_jfh() + self.assertEqual(len(jfh.JFH), + len(SWEEP_ANGLES) * len(SWEEP_DISTANCES)) + + for index, entry in enumerate(jfh.JFH): + distance = SWEEP_DISTANCES[index // len(SWEEP_ANGLES)] + angle = float(SWEEP_ANGLES[index % len(SWEEP_ANGLES)]) + position, dcm = pose_for(angle, distance, [0.0, 0.0, 0.0], + [0.0, 0.0, 1.0], [1.0, 0.0, 0.0]) + # The JFH file stores positions to 9 and DCMs to 6 significant + # digits, so file precision is the tolerance here. + np.testing.assert_allclose(entry['xyz'], position, atol=1e-7) + np.testing.assert_allclose(entry['dcm'], dcm, atol=1e-6) + + +class _StubFuelManager: + def calc_delta_v(self, dt, v_e, m_dot_sum, m_o): + return 0.1 + + +class _StubGrouping: + def calc_m_dot_sum(self, group): + return 0.01 + + def calc_v_e(self, group): + return 3000.0 + + +class _StubVehicle: + """Minimal stand-in for the attributes compute_1d_approach reads.""" + + mass = 1000.0 + rcs_groups = {'neg_x': ['T1']} + thruster_data = {'T1': {'type': ['R4D']}} + thruster_metrics = {'R4D': {'MIB': 0.4, 'F': 400.0}} + + +class DynamicsPathFiringCount(unittest.TestCase): + """The dynamics-driven approach honors an exact firing count too.""" + + def _run(self, n_firings): + return compute_1d_approach( + inputs=ApproachInputs(v_ida=0.0, v_o=5.0, r_o=20.0), + vv=_StubVehicle(), fuel_mgr=_StubFuelManager(), + grouping=_StubGrouping(), cant_rad=0.0, + dt_strategy={'multiplier': 1.0}, n_firings=n_firings) + + def test_requested_count_is_exact(self): + for n_firings in (1, 3, 10): + with self.subTest(n_firings=n_firings): + results = self._run(n_firings) + self.assertEqual(len(results['t']), n_firings) + self.assertEqual(len(results['x']), n_firings) + self.assertEqual(len(results['rot']), n_firings) + + def test_invalid_count_is_rejected(self): + for value in (0, -3, 2.5): + with self.subTest(value=value): + with pytest.raises(ValueError): + self._run(value) + + def test_unreachable_count_is_reported_not_silently_shortened(self): + # 5 m/s of delta-v at 0.1 m/s per firing completes in ~50 firings. + with pytest.raises(ValueError) as excinfo: + self._run(500) + self.assertIn('n_firings', str(excinfo.value)) + + def test_omitting_the_count_keeps_the_original_behavior(self): + results = compute_1d_approach( + inputs=ApproachInputs(v_ida=0.0, v_o=5.0, r_o=20.0), + vv=_StubVehicle(), fuel_mgr=_StubFuelManager(), + grouping=_StubGrouping(), cant_rad=0.0, + dt_strategy={'multiplier': 1.0}) + # Unbounded, the loop runs until the required delta-v is spent: 5.0 + # m/s at 0.1 m/s per firing (one extra iteration from the residual + # of the accumulated subtraction), plus the initial state entry. + self.assertEqual(len(results['t']), 52) + + +if __name__ == '__main__': + unittest.main() diff --git a/tests/mdao/mdao_unit_test_04.py b/tests/mdao/mdao_unit_test_04.py new file mode 100644 index 0000000..a3f8fdb --- /dev/null +++ b/tests/mdao/mdao_unit_test_04.py @@ -0,0 +1,370 @@ +# ======================== +# PyRPOD: tests/mdao/mdao_unit_test_04.py +# ======================== +# Unit tests for the surface-load integration of pyrpod.mdao.surface_loads, +# on small meshes whose loading has a known closed-form resultant: +# +# * force integration: uniform pressure on a flat plate gives -p*A*n_hat, +# and shear acts along the tangential projection of the flow direction; +# * moment about a user-defined reference point: M = (r_c - r_ref) x F, +# checked against a hand-computed value and against the plate's own +# linear pressure distribution; +# * center of pressure: recovered exactly for a linearly varying pressure +# field, consistent with the reported moment, and reported as unavailable +# (with a status naming the reason) for zero load and for near-total +# cancellation; +# * coefficients: computed when the normalization inputs are complete, +# omitted entirely when they are not -- no invented defaults; +# * component selection: a bounds-selected pair of components reproduces +# the whole-mesh resultant when summed. +# +# Everything here is constructed in memory; no case files are read. +# +# Run: python -m pytest mdao/mdao_unit_test_04.py (from tests/) + +import unittest + +import numpy as np +import pytest + +from pyrpod.mdao.study_config import ComponentSpec, Normalization +from pyrpod.mdao.surface_loads import ( + face_areas, + flow_directions, + integrate_component_loads, + select_component_faces, +) + + +def unit_square_plate(n=4, half=1.0, z=0.0): + """Triangulated square plate in the z-plane, unit normals along +Z. + + Returns (vectors, centroids, normals, areas) for an n x n grid of cells + split into two triangles each. + """ + edges = np.linspace(-half, half, n + 1) + faces = [] + for i in range(n): + for j in range(n): + x0, x1 = edges[i], edges[i + 1] + y0, y1 = edges[j], edges[j + 1] + # Counter-clockwise seen from +Z, so the unit normal is +Z. + faces.append([[x0, y0, z], [x1, y0, z], [x1, y1, z]]) + faces.append([[x0, y0, z], [x1, y1, z], [x0, y1, z]]) + vectors = np.asarray(faces, dtype=float) + centroids = vectors.mean(axis=1) + normals = np.tile([0.0, 0.0, 1.0], (len(vectors), 1)) + return vectors, centroids, normals, face_areas(vectors) + + +def integrate(centroids, normals, areas, pressures, shears=None, + heat=None, moment_ref=(0.0, 0.0, 0.0), + source=(0.0, 0.0, 5.0), normalization=None, strikes=None): + n_faces = len(centroids) + zeros = np.zeros(n_faces) + return integrate_component_loads( + component_name='plate', + face_indices=np.arange(n_faces), + centroids=centroids, unit_normals=normals, areas=areas, + pressures=np.asarray(pressures, dtype=float), + shear_stresses=zeros if shears is None + else np.asarray(shears, dtype=float), + heat_fluxes=zeros if heat is None else np.asarray(heat, dtype=float), + strikes=strikes, + moment_reference_point=moment_ref, + source_position=source, + normalization=normalization) + + +class ForceIntegration(unittest.TestCase): + + def test_uniform_pressure_on_a_flat_plate(self): + vectors, centroids, normals, areas = unit_square_plate(n=4, half=1.0) + pressure = 250.0 + loads = integrate(centroids, normals, areas, + np.full(len(centroids), pressure)) + + total_area = 4.0 # 2 m x 2 m plate + self.assertAlmostEqual(loads.total_area, total_area, places=12) + # Pressure pushes into the surface: along -n_hat. + np.testing.assert_allclose(loads.force, + [0.0, 0.0, -pressure * total_area], + atol=1e-9) + np.testing.assert_allclose(loads.pressure_force, loads.force, + atol=1e-12) + np.testing.assert_allclose(loads.shear_force, np.zeros(3), atol=1e-12) + self.assertAlmostEqual(loads.force_magnitude, pressure * total_area, + places=9) + + def test_shear_acts_along_the_tangential_flow_direction(self): + vectors, centroids, normals, areas = unit_square_plate(n=2, half=1.0) + shear = 3.0 + # Flow arriving at 45 deg in the X-Z plane: the tangential + # projection on the plate is +X for every face. + flow = np.tile([np.sqrt(0.5), 0.0, -np.sqrt(0.5)], (len(centroids), 1)) + loads = integrate_component_loads( + component_name='plate', face_indices=np.arange(len(centroids)), + centroids=centroids, unit_normals=normals, areas=areas, + pressures=np.zeros(len(centroids)), + shear_stresses=np.full(len(centroids), shear), + heat_fluxes=np.zeros(len(centroids)), strikes=None, + moment_reference_point=[0.0, 0.0, 0.0], flow_unit_vectors=flow) + + np.testing.assert_allclose(loads.shear_force, + [shear * 4.0, 0.0, 0.0], atol=1e-9) + np.testing.assert_allclose(loads.pressure_force, np.zeros(3), + atol=1e-12) + + def test_affected_area_counts_struck_faces_only(self): + vectors, centroids, normals, areas = unit_square_plate(n=2, half=1.0) + strikes = np.zeros(len(centroids)) + strikes[:4] = 1.0 + loads = integrate(centroids, normals, areas, + np.full(len(centroids), 10.0), strikes=strikes) + + self.assertEqual(loads.n_struck_faces, 4) + self.assertAlmostEqual(loads.affected_area, float(np.sum(areas[:4])), + places=12) + self.assertAlmostEqual(loads.total_area, 4.0, places=12) + + +class MomentIntegration(unittest.TestCase): + + def test_moment_about_a_user_defined_point(self): + vectors, centroids, normals, areas = unit_square_plate(n=4, half=1.0) + pressure = 100.0 + reference = np.array([2.0, -1.0, 3.0]) + loads = integrate(centroids, normals, areas, + np.full(len(centroids), pressure), + moment_ref=reference) + + # Uniform load: the resultant acts at the plate centroid (origin). + force = np.array([0.0, 0.0, -pressure * 4.0]) + expected = np.cross(np.zeros(3) - reference, force) + np.testing.assert_allclose(loads.moment, expected, atol=1e-9) + np.testing.assert_allclose(loads.pressure_moment, expected, atol=1e-9) + np.testing.assert_allclose(loads.shear_moment, np.zeros(3), + atol=1e-12) + + def test_moment_reference_point_is_reported_back(self): + vectors, centroids, normals, areas = unit_square_plate(n=2, half=1.0) + loads = integrate(centroids, normals, areas, + np.full(len(centroids), 5.0), + moment_ref=(0.5, 0.25, 0.0)) + np.testing.assert_allclose(loads.moment_reference_point, + [0.5, 0.25, 0.0]) + + def test_zero_reference_moment_when_taken_at_the_load_center(self): + vectors, centroids, normals, areas = unit_square_plate(n=6, half=1.0) + loads = integrate(centroids, normals, areas, + np.full(len(centroids), 42.0), + moment_ref=(0.0, 0.0, 0.0)) + np.testing.assert_allclose(loads.moment, np.zeros(3), atol=1e-9) + + +class CenterOfPressure(unittest.TestCase): + + def test_recovered_for_a_linear_pressure_distribution(self): + # p(x) = p0 * (1 + x) on [-1, 1] has its area-weighted center at + # x = integral(x*p)/integral(p) = (2/3) / 2 = 1/3. + vectors, centroids, normals, areas = unit_square_plate(n=40, half=1.0) + pressures = 100.0 * (1.0 + centroids[:, 0]) + loads = integrate(centroids, normals, areas, pressures) + + self.assertEqual(loads.center_of_pressure_status, 'ok') + cop = loads.center_of_pressure + self.assertIsNotNone(cop) + # Face-centroid sampling of a continuous field is a midpoint rule, + # and the cell-diagonal split is not mirror-invariant, so both the + # recovered x and the residual y carry an O(h^2) = (2/40)^2 ~ 2e-4 + # discretization error. + self.assertAlmostEqual(float(cop[0]), 1.0 / 3.0, delta=5e-4) + self.assertAlmostEqual(float(cop[1]), 0.0, delta=5e-4) + self.assertAlmostEqual(float(cop[2]), 0.0, places=9) + # The auxiliary pressure-weighted centroid agrees for a planar + # component under unidirectional pressure. + np.testing.assert_allclose(loads.pressure_weighted_centroid, cop, + atol=1e-6) + + def test_is_consistent_with_the_reported_moment(self): + vectors, centroids, normals, areas = unit_square_plate(n=20, half=1.0) + pressures = 50.0 * (2.0 + centroids[:, 0] + 0.5 * centroids[:, 1]) + loads = integrate(centroids, normals, areas, pressures, + moment_ref=(0.3, -0.2, 0.1)) + + arm = loads.center_of_pressure - loads.moment_reference_point + np.testing.assert_allclose(np.cross(arm, loads.force), loads.moment, + atol=1e-6) + + def test_zero_load_returns_no_center_of_pressure(self): + vectors, centroids, normals, areas = unit_square_plate(n=2, half=1.0) + loads = integrate(centroids, normals, areas, + np.zeros(len(centroids))) + + np.testing.assert_allclose(loads.force, np.zeros(3), atol=0.0) + np.testing.assert_allclose(loads.moment, np.zeros(3), atol=0.0) + self.assertIsNone(loads.center_of_pressure) + self.assertEqual(loads.center_of_pressure_status, 'zero_load') + self.assertIsNone(loads.pressure_weighted_centroid) + self.assertEqual(loads.max_pressure, 0.0) + self.assertEqual(loads.affected_area, 0.0) + + def test_cancelling_load_is_reported_as_ill_conditioned(self): + # Two parallel plates with opposite normals and equal pressure: the + # face forces are large but the resultant cancels, so no line of + # action exists. A naive F x M / |F|^2 would explode here. + vectors, centroids, normals, areas = unit_square_plate(n=4, half=1.0) + vectors_b, centroids_b, normals_b, areas_b = unit_square_plate( + n=4, half=1.0, z=1.0) + centroids = np.vstack([centroids, centroids_b]) + normals = np.vstack([normals, -normals_b]) + areas = np.concatenate([areas, areas_b]) + loads = integrate(centroids, normals, areas, + np.full(len(centroids), 75.0)) + + self.assertLess(loads.force_magnitude, 1e-9) + self.assertIsNone(loads.center_of_pressure) + self.assertEqual(loads.center_of_pressure_status, 'ill_conditioned') + + def test_residual_couple_is_reported(self): + vectors, centroids, normals, areas = unit_square_plate(n=8, half=1.0) + pressures = 20.0 * (1.0 + centroids[:, 0]) + loads = integrate(centroids, normals, areas, pressures, + moment_ref=(0.0, 0.0, 0.0)) + # Pure pressure on a planar component: the moment is perpendicular + # to the force, so nothing is left over. + self.assertLess(loads.residual_couple, 1e-9) + + +class Coefficients(unittest.TestCase): + + def _loads(self, normalization): + vectors, centroids, normals, areas = unit_square_plate(n=4, half=1.0) + pressures = np.full(len(centroids), 10.0) + shears = np.full(len(centroids), 1.0) + heat = np.full(len(centroids), 500.0) + return integrate(centroids, normals, areas, pressures, shears=shears, + heat=heat, moment_ref=(1.0, 0.0, 0.0), + normalization=normalization) + + def test_computed_with_complete_normalization_inputs(self): + normalization = Normalization(reference_area=4.0, + reference_length=2.0, + dynamic_pressure=100.0, + reference_heat_flux=1000.0) + loads = self._loads(normalization) + + self.assertTrue(loads.has_coefficients) + self.assertAlmostEqual(loads.coefficients['CF'], + loads.force_magnitude / (100.0 * 4.0), + places=12) + self.assertAlmostEqual(loads.coefficients['CM'], + loads.moment_magnitude / (100.0 * 4.0 * 2.0), + places=12) + self.assertAlmostEqual(loads.coefficients['Cp_max'], 10.0 / 100.0, + places=12) + self.assertAlmostEqual(loads.coefficients['Cf_max'], 1.0 / 100.0, + places=12) + self.assertAlmostEqual(loads.coefficients['Cq_max'], 500.0 / 1000.0, + places=12) + + def test_omitted_entirely_when_no_normalization_is_supplied(self): + loads = self._loads(None) + self.assertFalse(loads.has_coefficients) + self.assertEqual(loads.coefficients, {}) + + def test_partial_normalization_yields_only_the_valid_coefficients(self): + # Dynamic pressure but no reference area: surface-load coefficients + # are available, force and moment coefficients are not. + loads = self._loads(Normalization(dynamic_pressure=100.0)) + self.assertIn('Cp_max', loads.coefficients) + self.assertNotIn('CF', loads.coefficients) + self.assertNotIn('CM', loads.coefficients) + + # Force inputs but no reference length: no moment coefficient. + loads = self._loads(Normalization(reference_area=4.0, + dynamic_pressure=100.0)) + self.assertIn('CF', loads.coefficients) + self.assertNotIn('CM', loads.coefficients) + self.assertNotIn('Cq_max', loads.coefficients) + + def test_invalid_normalization_values_are_rejected(self): + from pyrpod.mdao.study_config import StudyConfigError + + for payload in ({'reference_area': 0.0}, {'dynamic_pressure': -1.0}): + with self.subTest(payload=payload): + with pytest.raises(StudyConfigError): + Normalization.from_mapping(payload) + + +class ComponentSelection(unittest.TestCase): + + def test_components_partition_the_resultant(self): + vectors, centroids, normals, areas = unit_square_plate(n=6, half=1.0) + pressures = 30.0 * (2.0 + centroids[:, 0]) + + left = ComponentSpec.from_mapping( + {'name': 'left', 'bounds': {'min': [-1.0, -1.0, -0.1], + 'max': [0.0, 1.0, 0.1]}}) + right = ComponentSpec.from_mapping( + {'name': 'right', 'bounds': {'min': [0.0, -1.0, -0.1], + 'max': [1.0, 1.0, 0.1]}}) + whole = integrate(centroids, normals, areas, pressures) + + partial = np.zeros(3) + for component in (left, right): + indices = select_component_faces(component, centroids, + len(centroids)) + loads = integrate_component_loads( + component_name=component.name, face_indices=indices, + centroids=centroids, unit_normals=normals, areas=areas, + pressures=pressures, + shear_stresses=np.zeros(len(centroids)), + heat_fluxes=np.zeros(len(centroids)), strikes=None, + moment_reference_point=[0.0, 0.0, 0.0], + source_position=[0.0, 0.0, 5.0]) + partial = partial + loads.force + self.assertLess(loads.n_faces, whole.n_faces) + + np.testing.assert_allclose(partial, whole.force, atol=1e-9) + + def test_explicit_face_indices(self): + vectors, centroids, normals, areas = unit_square_plate(n=2, half=1.0) + component = ComponentSpec.from_mapping( + {'name': 'strip', 'face_indices': [0, 1, 2]}) + indices = select_component_faces(component, centroids, len(centroids)) + np.testing.assert_array_equal(indices, [0, 1, 2]) + + def test_out_of_range_and_empty_selections_are_rejected(self): + vectors, centroids, normals, areas = unit_square_plate(n=2, half=1.0) + with pytest.raises(ValueError): + select_component_faces( + ComponentSpec.from_mapping({'name': 'bad', + 'face_indices': [999]}), + centroids, len(centroids)) + with pytest.raises(ValueError): + select_component_faces( + ComponentSpec.from_mapping( + {'name': 'empty', + 'bounds': {'min': [10.0, 10.0, 10.0], + 'max': [11.0, 11.0, 11.0]}}), + centroids, len(centroids)) + + +class FlowDirections(unittest.TestCase): + + def test_unit_vectors_point_from_the_source_to_each_face(self): + centroids = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0]]) + directions = flow_directions(centroids, [0.0, 0.0, 2.0]) + np.testing.assert_allclose(directions[0], [0.0, 0.0, -1.0], atol=1e-12) + np.testing.assert_allclose(np.linalg.norm(directions, axis=1), + [1.0, 1.0], atol=1e-12) + + def test_face_at_the_source_yields_a_zero_direction(self): + directions = flow_directions(np.array([[0.0, 0.0, 2.0]]), + [0.0, 0.0, 2.0]) + np.testing.assert_allclose(directions[0], np.zeros(3)) + + +if __name__ == '__main__': + unittest.main() diff --git a/tests/mdao/mdao_unit_test_05.py b/tests/mdao/mdao_unit_test_05.py new file mode 100644 index 0000000..3b88da6 --- /dev/null +++ b/tests/mdao/mdao_unit_test_05.py @@ -0,0 +1,351 @@ +# ======================== +# PyRPOD: tests/mdao/mdao_unit_test_05.py +# ======================== +# Unit tests for the YAML study configuration (pyrpod.mdao.study_config): +# +# * the committed flat-plate example configurations parse, and every field +# the study engine relies on comes back with the expected value; +# * paths are resolved relative to the configuration file, and the case +# must be a real PyRPOD case (config.ini present); +# * validation rejects what must not be guessed: an unsupported plume +# model, a bad firing count, a firing list that disagrees with +# n_firings, non-orthogonal target axes, duplicate component names and +# non-positive normalization values; +# * incomplete normalization inputs disable the corresponding +# coefficients instead of inventing defaults; +# * backward compatibility: the study layer never touches the case's own +# config.ini, which keeps parsing exactly as before. +# +# Run: python -m pytest mdao/mdao_unit_test_05.py (from tests/) + +import configparser +import copy +import os +import tempfile +import unittest +from pathlib import Path + +import numpy as np +import pytest +import yaml + +from pyrpod.mdao.study_config import ( + Normalization, + StudyConfig, + StudyConfigError, + SUPPORTED_PLUME_MODEL, +) + +_TESTS_DIR = Path(__file__).resolve().parents[1] +CASE_DIR = _TESTS_DIR.parent / 'case' / 'plume' / 'plume_flat_plate_sweep' +BASELINE_YAML = CASE_DIR / 'study' / 'flat_plate_baseline.yaml' +SWEEP_YAML = CASE_DIR / 'study' / 'flat_plate_sweep.yaml' +SINGLE_JFH_YAML = CASE_DIR / 'study' / 'flat_plate_sweep_single_jfh.yaml' + + +def baseline_mapping(): + return yaml.safe_load(BASELINE_YAML.read_text(encoding='utf-8')) + + +def from_mapping(data): + return StudyConfig.from_mapping(data, source_path=str(BASELINE_YAML)) + + +class CommittedExampleConfigurations(unittest.TestCase): + + def test_baseline_configuration_parses(self): + config = StudyConfig.from_yaml(BASELINE_YAML) + + self.assertEqual(config.study_name, 'cai2016_flat_plate_baseline') + self.assertEqual(config.plume_model, SUPPORTED_PLUME_MODEL) + self.assertTrue(os.path.isdir(config.case_dir)) + self.assertTrue(config.case_dir.endswith(os.sep)) + self.assertTrue(os.path.isfile(os.path.join(config.case_dir, + 'config.ini'))) + self.assertEqual(config.source_path, str(BASELINE_YAML)) + + self.assertEqual(config.sweep.plate_angles_deg, (0.0,)) + self.assertEqual(config.sweep.source_distances, (4.0,)) + self.assertEqual(config.sweep.n_firings, 1) + self.assertEqual(config.sweep.thrusters, (1,)) + self.assertEqual(config.n_cases, 1) + + np.testing.assert_allclose(config.target.reference_point, + [0.0, 0.0, 0.0]) + np.testing.assert_allclose(config.target.normal, [0.0, 0.0, 1.0]) + np.testing.assert_allclose(config.target.tangent, [1.0, 0.0, 0.0]) + self.assertEqual([c.name for c in config.target.components], ['plate']) + + np.testing.assert_allclose(config.loads.moment_reference_point, + [0.0, 0.0, 0.0]) + self.assertTrue(config.loads.normalization.has_moment_inputs) + self.assertTrue(config.output.write_vtk) + self.assertFalse(config.output.write_plots) + self.assertEqual(config.units['force'], 'N') + + def test_sweep_configuration_enumerates_every_combination(self): + config = StudyConfig.from_yaml(SWEEP_YAML) + self.assertEqual(len(config.sweep.plate_angles_deg), 19) + self.assertEqual(len(config.sweep.source_distances), 5) + self.assertEqual(config.n_cases, 95) + self.assertTrue(config.output.write_plots) + + def test_provenance_records_the_configuration_source(self): + config = StudyConfig.from_yaml(BASELINE_YAML) + provenance = config.provenance() + self.assertEqual(provenance['config_path'], str(BASELINE_YAML)) + self.assertEqual(provenance['plume_model'], SUPPORTED_PLUME_MODEL) + self.assertIn('units', provenance) + + def test_output_directory_can_be_overridden(self): + config = StudyConfig.from_yaml(BASELINE_YAML) + with tempfile.TemporaryDirectory() as tmp: + moved = config.with_output_dir(tmp) + self.assertEqual(moved.output_dir, os.path.abspath(tmp)) + # The original is untouched (the dataclass is frozen). + self.assertNotEqual(config.output_dir, moved.output_dir) + + +class PathResolution(unittest.TestCase): + + def test_relative_paths_resolve_against_the_configuration_file(self): + config = StudyConfig.from_yaml(BASELINE_YAML) + self.assertEqual(Path(config.case_dir).resolve(), CASE_DIR.resolve()) + + def test_missing_case_directory_is_rejected(self): + data = baseline_mapping() + data['study']['case_dir'] = 'no/such/case' + with pytest.raises(StudyConfigError) as excinfo: + from_mapping(data) + self.assertIn('case_dir', str(excinfo.value)) + + def test_directory_without_config_ini_is_rejected(self): + data = baseline_mapping() + with tempfile.TemporaryDirectory() as tmp: + data['study']['case_dir'] = tmp + with pytest.raises(StudyConfigError) as excinfo: + from_mapping(data) + self.assertIn('config.ini', str(excinfo.value)) + + def test_missing_configuration_file_is_reported(self): + with pytest.raises(StudyConfigError): + StudyConfig.from_yaml(CASE_DIR / 'study' / 'does_not_exist.yaml') + + +class Validation(unittest.TestCase): + + def test_study_name_is_required(self): + data = baseline_mapping() + del data['study']['name'] + with pytest.raises(StudyConfigError): + from_mapping(data) + + def test_only_the_supported_plume_model_is_accepted(self): + data = baseline_mapping() + data['plume_model']['name'] = 'CollisionlessGasKinetics' + with pytest.raises(StudyConfigError) as excinfo: + from_mapping(data) + self.assertIn(SUPPORTED_PLUME_MODEL, str(excinfo.value)) + + def test_invalid_firing_counts_are_rejected(self): + for value in (0, -2, 1.5, 'many'): + with self.subTest(value=value): + data = baseline_mapping() + data['sweep']['n_firings'] = value + with pytest.raises(StudyConfigError): + from_mapping(data) + + def test_explicit_firings_must_match_n_firings(self): + data = baseline_mapping() + data['sweep']['n_firings'] = 3 + data['sweep']['firings'] = [ + {'position': [0.0, 0.0, 4.0], 'dcm': np.eye(3).tolist()}] + with pytest.raises(StudyConfigError) as excinfo: + from_mapping(data) + self.assertIn('n_firings', str(excinfo.value)) + + def test_explicit_firings_are_parsed_when_the_count_agrees(self): + data = baseline_mapping() + data['sweep']['n_firings'] = 2 + data['sweep']['firings'] = [ + {'position': [0.0, 0.0, 4.0], 'dcm': np.eye(3).tolist(), + 'thrusters': [1], 'duration_s': 0.25}, + {'position': [0.5, 0.0, 4.0], 'dcm': np.eye(3).tolist()}] + config = from_mapping(data) + self.assertEqual(len(config.sweep.firings), 2) + self.assertEqual(config.sweep.firings[0].duration_s, 0.25) + # The second firing inherits the sweep defaults. + self.assertEqual(config.sweep.firings[1].duration_s, 1.0) + self.assertEqual(config.sweep.firings[1].thrusters, (1,)) + + def test_malformed_firing_pose_is_rejected(self): + data = baseline_mapping() + data['sweep']['n_firings'] = 1 + data['sweep']['firings'] = [ + {'position': [0.0, 0.0, 4.0], 'dcm': [[1.0, 0.0], [0.0, 1.0]]}] + with pytest.raises(StudyConfigError): + from_mapping(data) + + def test_target_axes_must_be_orthogonal(self): + data = baseline_mapping() + data['target']['tangent'] = [0.0, 0.0, 1.0] + with pytest.raises(StudyConfigError) as excinfo: + from_mapping(data) + self.assertIn('orthogonal', str(excinfo.value)) + + def test_empty_sweep_axes_are_rejected(self): + for key in ('plate_angles_deg', 'source_distances'): + with self.subTest(key=key): + data = baseline_mapping() + data['sweep'][key] = [] + with pytest.raises(StudyConfigError): + from_mapping(data) + + def test_non_positive_distances_are_rejected(self): + data = baseline_mapping() + data['sweep']['source_distances'] = [4.0, 0.0] + with pytest.raises(StudyConfigError): + from_mapping(data) + + def test_duplicate_component_names_are_rejected(self): + data = baseline_mapping() + data['target']['components'] = [{'name': 'plate'}, {'name': 'plate'}] + with pytest.raises(StudyConfigError): + from_mapping(data) + + def test_component_bounds_are_validated(self): + data = baseline_mapping() + data['target']['components'] = [ + {'name': 'plate', + 'bounds': {'min': [1.0, 1.0, 1.0], 'max': [0.0, 0.0, 0.0]}}] + with pytest.raises(StudyConfigError): + from_mapping(data) + + +class SweepMode(unittest.TestCase): + """per_case (default) vs single_jfh, and what n_firings means in each.""" + + def test_default_mode_is_per_case(self): + self.assertEqual(StudyConfig.from_yaml(BASELINE_YAML).sweep.mode, + 'per_case') + + def test_single_jfh_configuration_parses(self): + config = StudyConfig.from_yaml(SINGLE_JFH_YAML) + self.assertEqual(config.sweep.mode, 'single_jfh') + self.assertEqual(config.n_cases, 95) + self.assertEqual(config.sweep.n_firings, 1) + # n_firings is per pose; the one history holds poses x n_firings. + self.assertEqual(config.sweep.total_firings, 95) + + def test_unknown_mode_is_rejected(self): + data = baseline_mapping() + data['sweep']['mode'] = 'one_big_jfh' + with pytest.raises(StudyConfigError) as excinfo: + from_mapping(data) + self.assertIn('sweep.mode', str(excinfo.value)) + + def test_total_firings_scales_with_firings_per_pose(self): + data = baseline_mapping() + data['sweep']['mode'] = 'single_jfh' + data['sweep']['plate_angles_deg'] = [-30.0, 0.0, 30.0] + data['sweep']['source_distances'] = [2.0, 4.0] + data['sweep']['n_firings'] = 4 + sweep = from_mapping(data).sweep + self.assertEqual(len(sweep.poses), 6) + self.assertEqual(sweep.total_firings, 24) + + def test_poses_are_enumerated_distance_major(self): + data = baseline_mapping() + data['sweep']['plate_angles_deg'] = [-10.0, 10.0] + data['sweep']['source_distances'] = [2.0, 4.0] + self.assertEqual(from_mapping(data).sweep.poses, + ((-10.0, 2.0), (10.0, 2.0), + (-10.0, 4.0), (10.0, 4.0))) + + def test_explicit_firing_count_is_checked_per_mode(self): + one_firing = [{'position': [0.0, 0.0, 4.0], + 'dcm': np.eye(3).tolist()}] + + # per_case: one firing per case is enough for any number of poses. + data = baseline_mapping() + data['sweep']['plate_angles_deg'] = [-10.0, 10.0] + data['sweep']['n_firings'] = 1 + data['sweep']['firings'] = one_firing + self.assertEqual(len(from_mapping(data).sweep.firings), 1) + + # single_jfh: the sequence must cover every pose. + data['sweep']['mode'] = 'single_jfh' + with pytest.raises(StudyConfigError) as excinfo: + from_mapping(data) + self.assertIn('single_jfh', str(excinfo.value)) + + data['sweep']['firings'] = one_firing * 2 + self.assertEqual(len(from_mapping(data).sweep.firings), 2) + + +class CoefficientNormalization(unittest.TestCase): + + def test_complete_inputs_enable_force_and_moment_coefficients(self): + normalization = StudyConfig.from_yaml( + BASELINE_YAML).loads.normalization + self.assertTrue(normalization.has_force_inputs) + self.assertTrue(normalization.has_moment_inputs) + self.assertEqual(normalization.reference_area, 64.0) + self.assertEqual(normalization.reference_length, 4.0) + + def test_absent_normalization_leaves_every_input_unset(self): + data = baseline_mapping() + del data['loads']['normalization'] + normalization = from_mapping(data).loads.normalization + self.assertFalse(normalization.has_force_inputs) + self.assertFalse(normalization.has_moment_inputs) + self.assertEqual(normalization.to_dict(), + {'reference_area': None, 'reference_length': None, + 'dynamic_pressure': None, + 'reference_heat_flux': None}) + + def test_partial_normalization_disables_only_what_it_must(self): + data = baseline_mapping() + del data['loads']['normalization']['reference_length'] + normalization = from_mapping(data).loads.normalization + self.assertTrue(normalization.has_force_inputs) + self.assertFalse(normalization.has_moment_inputs) + + def test_non_positive_values_are_rejected(self): + with pytest.raises(StudyConfigError): + Normalization.from_mapping({'dynamic_pressure': 0.0}) + + +class BackwardCompatibility(unittest.TestCase): + """The study layer sits on top of the case; it changes nothing below.""" + + def test_case_config_ini_still_parses_unchanged(self): + config = configparser.ConfigParser() + config.read(str(CASE_DIR / 'config.ini')) + self.assertEqual(config['pm']['kinetics'], 'Simplified') + self.assertEqual(config['tv']['stl'], 'flat_plate_transformed.stl') + self.assertEqual(config['jfh']['jfh'], 'jfh_flat_plate_sweep.A') + + def test_geometry_id_defaults_to_the_case_target_stl(self): + data = baseline_mapping() + del data['target']['geometry_id'] + self.assertEqual(from_mapping(data).target.geometry_id, + 'flat_plate_transformed.stl') + + def test_defaults_apply_to_a_minimal_configuration(self): + minimal = {'study': {'name': 'minimal', + 'case_dir': str(CASE_DIR)}, + 'sweep': {'plate_angles_deg': [0.0], + 'source_distances': [4.0]}} + config = StudyConfig.from_mapping(copy.deepcopy(minimal), + source_path=str(BASELINE_YAML)) + self.assertEqual(config.sweep.n_firings, 1) + self.assertEqual(config.plume_model, SUPPORTED_PLUME_MODEL) + self.assertEqual([c.name for c in config.target.components], + ['target']) + self.assertFalse(config.loads.normalization.has_force_inputs) + self.assertTrue(config.output.write_vtk) + + +if __name__ == '__main__': + unittest.main() diff --git a/tests/mdao/mdao_unit_test_06.py b/tests/mdao/mdao_unit_test_06.py new file mode 100644 index 0000000..99aa838 --- /dev/null +++ b/tests/mdao/mdao_unit_test_06.py @@ -0,0 +1,324 @@ +# ======================== +# PyRPOD: tests/mdao/mdao_unit_test_06.py +# ======================== +# Unit tests for the generic external-reference comparison +# (pyrpod.mdao.reference_data) and the structured result schema +# (pyrpod.mdao.study_results): +# +# * the metrics themselves -- absolute error, relative error, normalized +# RMSE, peak error, integrated-load error and center-of-pressure +# displacement -- against hand-computed values, including the cases +# where a metric is undefined and must return None instead of infinity; +# * loading reference datasets from CSV, JSON and YAML with the SAME +# comparison behavior, since the interface must not know (or care) +# whether the data came from DSMC, an analytical solution, an experiment +# or another code; +# * matching by case keys, and the reporting of unmatched cases and of +# quantities the candidate does not provide -- nothing is fabricated; +# * the result schema round-tripping to CSV rows and JSON. +# +# Run: python -m pytest mdao/mdao_unit_test_06.py (from tests/) + +import json +import math +import tempfile +import unittest +from pathlib import Path + +import numpy as np +import pytest +import yaml + +from pyrpod.mdao.reference_data import ( + ReferenceDataset, + ReferenceRecord, + absolute_error, + center_of_pressure_displacement, + compare_case, + compare_results, + integrated_load_error, + load_reference_dataset, + normalized_rmse, + peak_error, + relative_error, +) +from pyrpod.mdao.study_results import CaseResult, StudyResults +from pyrpod.mdao.surface_loads import ComponentLoads + + +def make_result(case_id='case000', component='plate', force=(0.0, 0.0, -10.0), + moment=(0.0, 2.0, 0.0), cop=(0.1, 0.0, 0.0), + plate_angle_deg=0.0, source_distance=4.0, + coefficients=None): + loads = ComponentLoads( + component=component, n_faces=8, n_struck_faces=8, total_area=4.0, + affected_area=4.0, + pressure_force=np.asarray(force, dtype=float), + shear_force=np.zeros(3), + force=np.asarray(force, dtype=float), + force_magnitude=float(np.linalg.norm(force)), + moment_reference_point=np.zeros(3), + pressure_moment=np.asarray(moment, dtype=float), + shear_moment=np.zeros(3), + moment=np.asarray(moment, dtype=float), + moment_magnitude=float(np.linalg.norm(moment)), + center_of_pressure=(None if cop is None + else np.asarray(cop, dtype=float)), + center_of_pressure_status='zero_load' if cop is None else 'ok', + residual_couple=0.0, + pressure_weighted_centroid=None if cop is None + else np.asarray(cop, dtype=float), + max_pressure=12.0, max_shear_stress=1.5, max_heat_flux=300.0, + total_heat_load=1200.0, + coefficients=dict(coefficients or {})) + return CaseResult.from_loads( + loads, study_name='unit', case_id=case_id, firing_id=1, + geometry_id='plate.stl', mesh_faces=8, + coordinate_system='case global frame', units={'force': 'N'}, + plume_source_position=[0.0, 0.0, 4.0], + plume_source_orientation=list(np.eye(3).ravel()), + target_normal=[0.0, 0.0, 1.0], target_tangent=[1.0, 0.0, 0.0], + target_reference_point=[0.0, 0.0, 0.0], + plate_angle_deg=plate_angle_deg, source_distance=source_distance, + firing_duration_s=1.0, thrusters=[1], + plume_model='SimplifiedGasKinetics', plume_model_parameters={}) + + +class Metrics(unittest.TestCase): + + def test_absolute_error_of_scalars_and_vectors(self): + self.assertAlmostEqual(absolute_error(5.0, 4.0), 1.0) + self.assertAlmostEqual(absolute_error([3.0, 4.0, 0.0], + [0.0, 0.0, 0.0]), 5.0) + + def test_relative_error_is_none_for_a_zero_reference(self): + self.assertAlmostEqual(relative_error(11.0, 10.0), 0.1) + self.assertIsNone(relative_error(1.0, 0.0)) + + def test_normalized_rmse_uses_the_reference_range(self): + reference = np.array([0.0, 1.0, 2.0, 3.0]) + candidate = reference + 0.3 + # RMSE 0.3 over a range of 3.0 + self.assertAlmostEqual(normalized_rmse(candidate, reference), 0.1) + self.assertAlmostEqual( + normalized_rmse(candidate, reference, norm='mean'), + 0.3 / 1.5) + + def test_normalized_rmse_is_none_when_undefined(self): + constant = np.array([2.0, 2.0, 2.0]) + self.assertIsNone(normalized_rmse(constant + 0.1, constant)) + self.assertIsNone(normalized_rmse([1.0, 2.0], [1.0, 2.0, 3.0])) + + def test_peak_error(self): + absolute, relative = peak_error([1.0, 4.0, 2.0], [1.0, 5.0, 2.0]) + self.assertAlmostEqual(absolute, 1.0) + self.assertAlmostEqual(relative, 0.2) + + def test_integrated_load_error_compares_vectors_not_magnitudes(self): + # Same magnitude, opposite direction: a magnitude-only comparison + # would report zero error. + absolute, relative = integrated_load_error([0.0, 0.0, 10.0], + [0.0, 0.0, -10.0]) + self.assertAlmostEqual(absolute, 20.0) + self.assertAlmostEqual(relative, 2.0) + + def test_center_of_pressure_displacement(self): + distance, normalized = center_of_pressure_displacement( + [0.3, 0.0, 0.0], [0.0, 0.4, 0.0], reference_length=2.0) + self.assertAlmostEqual(distance, 0.5) + self.assertAlmostEqual(normalized, 0.25) + + distance, normalized = center_of_pressure_displacement( + [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]) + self.assertAlmostEqual(distance, 0.0) + self.assertIsNone(normalized) + + +class DatasetLoading(unittest.TestCase): + """The same records, expressed in any supported format, compare alike.""" + + QUANTITIES = {'force': [0.0, 0.0, -9.0], 'max_pressure': 11.0} + + def _csv(self, directory): + path = Path(directory) / 'reference.csv' + path.write_text( + 'case_id,component,plate_angle_deg,source_distance,' + 'force_x,force_y,force_z,max_pressure\n' + 'case000,plate,0.0,4.0,0.0,0.0,-9.0,11.0\n', + encoding='utf-8') + return path + + def _json(self, directory): + path = Path(directory) / 'reference.json' + path.write_text(json.dumps({ + 'label': 'independent solver', + 'source': 'run-42', + 'units': {'force': 'N'}, + 'records': [{'key': {'case_id': 'case000', 'component': 'plate'}, + 'quantities': self.QUANTITIES}]}), encoding='utf-8') + return path + + def _yaml(self, directory): + path = Path(directory) / 'reference.yaml' + path.write_text(yaml.safe_dump({ + 'label': 'experiment', + 'records': [{'key': {'case_id': 'case000'}, + 'quantities': self.QUANTITIES}]}), encoding='utf-8') + return path + + def test_every_format_yields_the_same_comparison(self): + result = make_result() + with tempfile.TemporaryDirectory() as tmp: + for builder in (self._csv, self._json, self._yaml): + with self.subTest(builder=builder.__name__): + dataset = load_reference_dataset(builder(tmp)) + self.assertEqual(len(dataset), 1) + report = compare_results([result], dataset) + by_name = {c.quantity: c for c in report.comparisons} + self.assertEqual(set(by_name), {'force', 'max_pressure'}) + self.assertAlmostEqual(by_name['force'].absolute_error, + 1.0) + self.assertAlmostEqual( + by_name['max_pressure'].relative_error, + 1.0 / 11.0) + + def test_csv_vector_columns_are_assembled(self): + with tempfile.TemporaryDirectory() as tmp: + dataset = load_reference_dataset(self._csv(tmp)) + record = dataset.records[0] + self.assertEqual(record.quantities['force'], [0.0, 0.0, -9.0]) + self.assertEqual(record.key['case_id'], 'case000') + self.assertAlmostEqual(record.key['plate_angle_deg'], 0.0) + + def test_label_can_be_overridden_and_is_only_a_label(self): + with tempfile.TemporaryDirectory() as tmp: + dataset = load_reference_dataset(self._json(tmp), label='DSMC') + self.assertEqual(dataset.label, 'DSMC') + report = compare_results([make_result()], dataset) + # The metrics do not depend on where the data came from. + self.assertEqual(report.label, 'DSMC') + self.assertAlmostEqual( + report.max_relative_error('max_pressure'), 1.0 / 11.0) + + def test_unknown_format_and_missing_file_are_reported(self): + with tempfile.TemporaryDirectory() as tmp: + bad = Path(tmp) / 'reference.parquet' + bad.write_text('', encoding='utf-8') + with pytest.raises(ValueError): + load_reference_dataset(bad) + with pytest.raises(FileNotFoundError): + load_reference_dataset(Path(tmp) / 'absent.csv') + + +class Matching(unittest.TestCase): + + def test_records_match_on_swept_parameters(self): + dataset = ReferenceDataset(records=[ + ReferenceRecord(key={'plate_angle_deg': 30.0, + 'source_distance': 4.0}, + quantities={'max_pressure': 10.0})]) + matching = make_result(plate_angle_deg=30.0, source_distance=4.0) + other = make_result(plate_angle_deg=-30.0, source_distance=4.0) + + self.assertIsNotNone(dataset.match(matching)) + self.assertIsNone(dataset.match(other)) + + def test_unmatched_cases_are_listed_not_compared(self): + dataset = ReferenceDataset(records=[ + ReferenceRecord(key={'case_id': 'case000'}, + quantities={'max_pressure': 10.0})]) + report = compare_results( + [make_result(case_id='case000'), make_result(case_id='case001')], + dataset) + + self.assertEqual(len(report.comparisons), 1) + self.assertEqual(report.unmatched_cases, ['case001/plate']) + + def test_quantity_absent_from_the_candidate_is_flagged(self): + result = make_result() # no coefficients supplied + record = ReferenceRecord(key={'case_id': 'case000'}, + quantities={'CF': 0.5}) + comparison = compare_case(result, record)[0] + + self.assertEqual(comparison.status, 'missing_candidate') + self.assertIsNone(comparison.absolute_error) + self.assertEqual(comparison.reference, 0.5) + + def test_coefficients_are_compared_when_available(self): + result = make_result(coefficients={'CF': 0.4}) + record = ReferenceRecord(key={'case_id': 'case000'}, + quantities={'CF': 0.5}) + comparison = compare_case(result, record)[0] + + self.assertEqual(comparison.status, 'compared') + self.assertAlmostEqual(comparison.absolute_error, 0.1) + self.assertAlmostEqual(comparison.relative_error, 0.2) + + def test_shape_mismatch_is_reported_rather_than_coerced(self): + record = ReferenceRecord(key={'case_id': 'case000'}, + quantities={'force': [1.0, 2.0]}) + comparison = compare_case(make_result(), record)[0] + self.assertEqual(comparison.status, 'shape_mismatch') + + def test_center_of_pressure_displacement_is_recorded(self): + record = ReferenceRecord( + key={'case_id': 'case000'}, + quantities={'center_of_pressure': [0.0, 0.0, 0.0]}) + comparison = compare_case(make_result(cop=(0.1, 0.0, 0.0)), record, + reference_length=2.0)[0] + self.assertAlmostEqual(comparison.displacement, 0.1) + + def test_empty_dataset_produces_an_empty_report(self): + report = compare_results([make_result()], ReferenceDataset()) + self.assertEqual(len(report), 0) + self.assertEqual(report.unmatched_cases, ['case000/plate']) + with pytest.raises(ValueError): + report.write_csv(Path(tempfile.gettempdir()) / 'empty.csv') + + +class ResultSchema(unittest.TestCase): + + def test_row_expands_vectors_and_coefficients(self): + row = make_result(coefficients={'CF': 0.25}).to_row() + self.assertAlmostEqual(row['force_z'], -10.0) + self.assertAlmostEqual(row['moment_y'], 2.0) + self.assertAlmostEqual(row['center_of_pressure_x'], 0.1) + self.assertAlmostEqual(row['coeff_CF'], 0.25) + self.assertTrue(row['coefficients_available']) + + def test_missing_center_of_pressure_leaves_empty_columns(self): + row = make_result(cop=None).to_row() + self.assertEqual(row['center_of_pressure_x'], '') + self.assertEqual(row['center_of_pressure_status'], 'zero_load') + + def test_results_write_csv_and_json(self): + results = StudyResults(study_name='unit', output_dir='', + provenance={'plume_model': + 'SimplifiedGasKinetics'}, + cases=[make_result(case_id='case000'), + make_result(case_id='case001')]) + with tempfile.TemporaryDirectory() as tmp: + csv_path = results.write_csv(Path(tmp) / 'results.csv') + json_path = results.write_metadata(Path(tmp) / 'metadata.json') + + lines = Path(csv_path).read_text(encoding='utf-8').splitlines() + self.assertEqual(len(lines), 3) # header + two rows + self.assertIn('force_z', lines[0]) + + document = json.loads(Path(json_path).read_text(encoding='utf-8')) + self.assertEqual(document['n_cases'], 2) + self.assertEqual(document['provenance']['plume_model'], + 'SimplifiedGasKinetics') + self.assertEqual(document['cases'][0]['case_id'], 'case000') + self.assertIn('units', document['cases'][0]) + + def test_quantity_lookup_covers_the_comparable_fields(self): + result = make_result(coefficients={'CF': 0.25}) + self.assertEqual(result.quantity('force'), [0.0, 0.0, -10.0]) + self.assertTrue(math.isclose(result.quantity('max_heat_flux'), 300.0)) + self.assertAlmostEqual(result.quantity('CF'), 0.25) + self.assertIsNone(result.quantity('not_a_quantity')) + + +if __name__ == '__main__': + unittest.main() diff --git a/tests/test_manifest.yaml b/tests/test_manifest.yaml index f6b220a..be9a7c0 100644 --- a/tests/test_manifest.yaml +++ b/tests/test_manifest.yaml @@ -90,6 +90,75 @@ tests: manual_command: null collection_ignore_reason: null + - path: tests/mdao/mdao_unit_test_03.py + description: >- + Unit tests for prescribed firing generation and the exact meaning of + n_firings: invalid counts are rejected, N requested firings produce + exactly N JFH entries (read back through JetFiringHistory), an explicit + firing list that disagrees with n_firings is an error, a whole-sweep + sequence holds exactly poses x n_firings pose-tagged entries, the + generated pose convention reproduces the committed flat-plate sweep + JFH, and the dynamics-driven approach honors an exact count. + subsystem: mdao + category: unit + execution_mode: automated + development_status: implemented + collection_status: collected + reference: null + manual_command: null + collection_ignore_reason: null + + - path: tests/mdao/mdao_unit_test_04.py + description: >- + Unit tests for per-component surface-load integration on meshes with + analytically known loading: pressure and shear force integration, + moments about a user-defined reference point, center of pressure + (recovered, moment-consistent, and unavailable for zero-load and + near-cancellation cases), coefficient calculation and omission, and + component face selection. + subsystem: mdao + category: unit + execution_mode: automated + development_status: implemented + collection_status: collected + reference: null + manual_command: null + collection_ignore_reason: null + + - path: tests/mdao/mdao_unit_test_05.py + description: >- + Unit tests for the YAML study configuration: the committed flat-plate + examples parse, paths resolve against the configuration file, + validation rejects an unsupported plume model, an unknown sweep mode, + invalid firing counts, mismatched firing lists and bad geometry or + normalization inputs, the per-pose n_firings semantics hold in both + sweep modes, and the case's own config.ini keeps parsing unchanged. + subsystem: mdao + category: unit + execution_mode: automated + development_status: implemented + collection_status: collected + reference: null + manual_command: null + collection_ignore_reason: null + + - path: tests/mdao/mdao_unit_test_06.py + description: >- + Unit tests for the generic external-reference comparison (absolute and + relative error, normalized RMSE, peak error, integrated-load error, + center-of-pressure displacement; CSV / JSON / YAML datasets compared + identically regardless of origin; unmatched cases and missing + quantities reported rather than fabricated) and for the structured + result schema's CSV and JSON serialization. + subsystem: mdao + category: unit + execution_mode: automated + development_status: implemented + collection_status: collected + reference: null + manual_command: null + collection_ignore_reason: null + - path: tests/mdao/mdao_integration_test_01.py description: >- Placeholder MDAO integration test; the test body returns immediately @@ -103,6 +172,43 @@ tests: manual_command: null collection_ignore_reason: null + - path: tests/mdao/mdao_integration_test_02.py + description: >- + Baseline flat-plate case run end to end through the package-level + TradeStudy.from_config API: one case with exactly one JFH entry, the + standard per-face strike VTK at the advertised path, CSV and JSON + summaries with full provenance, physically consistent head-on loads, + available coefficients, and agreement of the integrated normal load + with the independent Cai 2016 exact reference. + subsystem: mdao + category: integration + execution_mode: automated + development_status: implemented + collection_status: collected + reference: >- + Cai 2016, Aerospace 3(4), 43, doi:10.3390/aerospace3040043, Eq. 15 + manual_command: null + collection_ignore_reason: null + + - path: tests/mdao/mdao_integration_test_03.py + description: >- + Multi-angle / multi-distance flat-plate sweep through the trade-study + API (a reduced subset of the committed sweep configuration): one result + per angle-distance case, mirror symmetry in +/- angle, load decay with + distance, center-of-pressure travel, agreement with the independent Cai + 2016 reference through the generic comparison interface, per-case VTK + isolation, optional trend plots, and the same machinery running on the + cylinder target with coefficients correctly unavailable. + subsystem: mdao + category: integration + execution_mode: automated + development_status: implemented + collection_status: collected + reference: >- + Cai 2016, Aerospace 3(4), 43, doi:10.3390/aerospace3040043, Eq. 15 + manual_command: null + collection_ignore_reason: null + - path: tests/mdao/mdao_verification_test_01.py description: >- OpenMDAO axial-positioning optimizer: minimizes the maximum heat-flux @@ -198,6 +304,25 @@ tests: manual_command: null collection_ignore_reason: null + - path: tests/mdao/mdao_integration_test_04.py + description: >- + The single-history sweep engine (sweep.mode: single_jfh): engine + dispatch from the configuration, exactly one Jet Firing History holding + poses x n_firings entries, one result record per firing keyed to its + pose, per-firing equivalence with the per-case engine, a single + results/strikes VTK series, the per-component sweep envelope (absent + rather than fabricated without the full pipeline), and the same + independent Cai 2016 reference comparison. + subsystem: mdao + category: integration + execution_mode: automated + development_status: implemented + collection_status: collected + reference: >- + Cai 2016, Aerospace 3(4), 43, doi:10.3390/aerospace3040043, Eq. 15 + manual_command: null + collection_ignore_reason: null + # ---------------------------------------------------------------- mission - path: tests/mission/mission_unit_test_01.py description: >- From 7e4422f2ec3785c7c1d0fe9c3d54e4563de7c657 Mon Sep 17 00:00:00 2001 From: andytorrestb Date: Mon, 27 Jul 2026 22:58:29 -0500 Subject: [PATCH 14/14] docs: document the plume-validation study docs/plume_validation_study.md covers the repaired TradeStudy architecture and the repairs made to the legacy sweeps, the package-level API, the two sweep decompositions and when to use each, the YAML configuration schema, the exact meaning of n_firings and prescribed firing lists, the force / moment / center-of-pressure definitions (including why a unique 3D center of pressure needs a stated convention and how degenerate cases are reported), coefficient normalization requirements, the result schema, VTK outputs, optional plots, the external reference-data format and how DSMC results plug into it later, how to add another validation geometry, and the known limitations. The root README gains a short pointer with the worked flat-plate example. Co-Authored-By: Claude Opus 5 --- README.md | 34 ++ docs/plume_validation_study.md | 622 +++++++++++++++++++++++++++++++++ 2 files changed, 656 insertions(+) create mode 100644 docs/plume_validation_study.md diff --git a/README.md b/README.md index 3b32ed2..098ad18 100644 --- a/README.md +++ b/README.md @@ -47,6 +47,40 @@ PyRPOD utilizies scientific libraries such as NumPy, SciPy, Matplotlib, and SymP python -m pytest tests/rpod # or just point pytest at a directory/file directly ``` +## Plume Validation Trade Studies + +Prescribed plume/target validation sweeps (fixed target, engineer-placed +firing poses, swept approach angle and source distance) run through a +package-level API on top of any existing case: + +```python +from pyrpod.mdao.TradeStudy import TradeStudy + +study = TradeStudy.from_config( + 'case/plume/plume_flat_plate_sweep/study/flat_plate_baseline.yaml') +results = study.run() # per-component force, moment, center of + # pressure, peak loads, VTK paths, CSV + JSON +``` + +A YAML study configuration adds the sweep, the exact JFH firing count, the +moment reference point and the coefficient normalization on top of the case's +own `config.ini`, which keeps owning every asset. `sweep.mode` picks how the +sweep is decomposed — `per_case` gives every angle-distance combination its +own Jet Firing History, `single_jfh` runs the whole sweep as one history with +a single strike series and a sweep-wide load envelope. Results can be compared +against independently generated reference data (DSMC, analytical, +experimental — the interface does not care which) through a generic +comparison layer. + +See [docs/plume_validation_study.md](docs/plume_validation_study.md) for the +architecture, the configuration schema, the exact meaning of `n_firings`, the +force / moment / center-of-pressure definitions, the result schema, the +reference-data format, and the known limitations. Worked examples: +[`flat_plate_baseline.yaml`](case/plume/plume_flat_plate_sweep/study/flat_plate_baseline.yaml), +[`flat_plate_sweep.yaml`](case/plume/plume_flat_plate_sweep/study/flat_plate_sweep.yaml) +and its single-history form +[`flat_plate_sweep_single_jfh.yaml`](case/plume/plume_flat_plate_sweep/study/flat_plate_sweep_single_jfh.yaml). + ## Test Reporting Test information lives in two places, and the split is deliberate: diff --git a/docs/plume_validation_study.md b/docs/plume_validation_study.md new file mode 100644 index 0000000..de48ec4 --- /dev/null +++ b/docs/plume_validation_study.md @@ -0,0 +1,622 @@ +# Prescribed plume-validation trade studies + +This document covers the repaired `TradeStudy` architecture and the +package-level API for **prescribed** plume/target validation sweeps: studies +in which the firing poses are placed by the engineer (or generated from a +swept approach angle and source distance about a stationary target) rather +than flown from vehicle dynamics. + +- [Quick start](#quick-start) +- [Architecture](#architecture) +- [Sweep modes: one JFH per case, or one for the sweep](#sweep-modes-one-jfh-per-case-or-one-for-the-sweep) +- [YAML configuration](#yaml-configuration) +- [`n_firings` and prescribed firings](#n_firings-and-prescribed-firings) +- [Integrated loads](#integrated-loads) + - [Force](#force) + - [Moment](#moment) + - [Center of pressure](#center-of-pressure) + - [Thermal quantities](#thermal-quantities) + - [Coefficients](#coefficients) +- [Result schema](#result-schema) +- [VTK outputs](#vtk-outputs) +- [Optional plots](#optional-plots) +- [External reference data](#external-reference-data) +- [Adding another validation geometry](#adding-another-validation-geometry) +- [Known limitations](#known-limitations) + +--- + +## Quick start + +A complete flat-plate example ships with the repository. It runs the Cai 2016 +Section-4 conditions (argon round jet, `D = 1 m`, `S0 = 2.0`, `T0 = 200 K`, +`Tw = 300 K`, fully diffuse) against the 8 m x 8 m plate of +`case/plume/plume_flat_plate_sweep`, with the plume source head-on at +`L = 4D`: + +```python +from pyrpod.mdao.TradeStudy import TradeStudy + +study = TradeStudy.from_config( + 'case/plume/plume_flat_plate_sweep/study/flat_plate_baseline.yaml') +results = study.run() + +case = results.cases[0] +print(case.force) # [~0, ~0, -2.763] N (plume pushes -Z) +print(case.moment) # ~0 about the plate center +print(case.center_of_pressure) # ~[0, 0, 0] +print(case.max_pressure) # 0.340 Pa +print(case.max_heat_flux) # 50.1 W/m^2 +print(case.coefficients['CF']) # 0.0391 +print(case.vtk_path) # .../cases/case000_.../results/strikes/firing-0.vtu +print(results.summary_csv_path) # .../case_results.csv +``` + +The multi-angle / multi-distance sweep (19 angles x 5 distances) is the same +call against `flat_plate_sweep.yaml`: + +```python +study = TradeStudy.from_config( + 'case/plume/plume_flat_plate_sweep/study/flat_plate_sweep.yaml') +results = study.run() +study.plot() # optional trend figures +``` + +Outputs land under the configured `output_dir` (by default +`/results/studies//`, which is gitignored): + +``` +results/studies/flat_plate_sweep/ + jfh/case000_alpha0p0_d4.A one prescribed JFH per case + cases/case000_alpha0p0_d4/results/strikes/firing-0.vtu + sweep_results.csv one row per case x component x firing + sweep_metadata.json provenance + nested case records + plots/ optional trend figures + reference_comparison.csv written when reference data is configured +``` + +To run a committed configuration into a scratch location (this is what the +automated tests do), pass `output_dir`: + +```python +TradeStudy.from_config(path, output_dir='/tmp/my_study').run() +``` + +--- + +## Architecture + +`TradeStudy` is a thin façade. The work lives in small modules with one +responsibility each, all under `pyrpod/mdao/`: + +| Module | Responsibility | +| --- | --- | +| `study_config.py` | Parse and validate the YAML study configuration | +| `firing_plan.py` | Build prescribed firings; write Jet Firing Histories | +| `plume_validation.py` | **Engine**: one JFH per case (`mode: per_case`) | +| `parameter_sweep.py` | **Engine**: one JFH for the sweep (`mode: single_jfh`) | +| `study_runtime.py` | Plumbing both engines share (assets, geometry, strikes, records) | +| `surface_loads.py` | Integrate per-face fields into component loads | +| `study_results.py` | Result schema; CSV + JSON output | +| `reference_data.py` | Generic external-reference comparison and metrics | +| `study_plots.py` | Optional sweep trend figures | + +Everything domain-specific is delegated to the existing PyRPOD objects — +`JetFiringHistory`, `TargetVehicle`, `VisitingVehicle`, `MissionEnvironment` +and `PlumeStrikeEstimationStudy` — so a study inherits the pipeline's +validation, logging, plume physics and VTK conventions instead of +reimplementing them. The single plume model is `SimplifiedGasKinetics`; the +configuration records it explicitly and rejects anything else. No plume-model +registry is introduced. + +`TradeStudy.from_config` picks the engine from `sweep.mode` (see the next +section); both do the same seven things, differing only in how many Jet +Firing Histories the sweep is decomposed into: + +1. build the prescribed firings; +2. write the Jet Firing History and read it back with the normal + `JetFiringHistory` parser; +3. run the plume-strike calculation; +4. integrate per-component loads; +5. export the per-face VTK fields; +6. record a structured result; +7. optionally compare against external reference data. + +### Repairs to the legacy sweeps + +The dynamics-driven sweeps (`run_axial_overshoot_sweep`, +`run_surface_cant_sweep`, `run_multi_var_sweep`) were broken before any +physics ran. Fixed on this branch: + +- `init_trade_study` constructed `PlumeStrikeEstimationStudy.RPOD(case_dir)` + — a class that does not exist — and passed a case directory where the study + takes a `MissionEnvironment`. It now builds the environment and the real + `PlumeStrikeEstimationStudy`. +- The sweeps called `jfh_plume_strikes(trade_study=True)`, which that + method's `(parallel, workers)` signature rejects with `TypeError`. +- `print_mission_report` read impingement maxima from study attributes that + are never set; it now summarizes them from the per-firing arrays + `jfh_plume_strikes` returns (the old attribute path still works for + existing callers). +- `interpret_mission_report` read `self.rpod.config`, which does not exist; + the configuration lives on the study's `MissionEnvironment`. + +### Execution paths + +With VTK output enabled (the default) each case runs the full +`PlumeStrikeEstimationStudy.jfh_plume_strikes()` pipeline, with the target's +output root temporarily redirected to that case's folder so sweep cases +cannot overwrite one another's artifacts. With VTK output disabled the study +calls the pipeline's own per-firing core, `compute_plume_strikes()` — the +same function `jfh_plume_strikes` calls per firing, so the numbers are +identical — which keeps large sweeps and automated tests fast and +artifact-free. This applies to both engines. + +--- + +## Sweep modes: one JFH per case, or one for the sweep + +The same configuration can be decomposed two ways, chosen with `sweep.mode`. +Per-firing numbers are **identical** either way — an automated test asserts +it pose by pose — so the choice is about what you want to look at, not about +accuracy. + +| | `per_case` (default) | `single_jfh` | +| --- | --- | --- | +| Engine | `PlumeValidationStudy` | `ParameterSweepStudy` | +| Jet Firing Histories | one per angle-distance case | **one**, spanning every pose | +| `case_id` / `firing_id` | one case per pose, firings 1..`n_firings` | one case, firings 1..`poses x n_firings` | +| Strike runs | one per case | one, over the whole history | +| VTK layout | `cases//results/strikes/firing-.vtu` | `results/strikes/firing-.vtu` — one series, scrubbable in ParaView | +| Cumulative fields | restart each case | accumulate across the sweep | +| Sweep envelope | n/a | recorded per component | + +Use `per_case` when each pose is an independent experiment (a validation +matrix compared pose-by-pose against reference data): nothing one pose does +can affect another's artifacts. Use `single_jfh` when the poses form one +sequence and you want the pipeline's cumulative fields to mean something: + +- `max_pressures` / `max_shears` — the worst load each face saw **anywhere** + in the sweep; +- `cum_strikes` — a coverage map over the whole sweep; +- `cum_heat_flux_load` — accumulated heat-flux dose. + +Those per-face fields live in the last firing's VTK file; their per-component +summary is exposed as `study.validation_study.envelope` and written into the +metadata document under `provenance.sweep_envelope`: + +```python +study = TradeStudy.from_config('.../flat_plate_sweep_single_jfh.yaml') +results = study.run() +study.validation_study.envelope['plate'] +# {'max_pressure': 0.340, 'max_shear_stress': 0.120, +# 'max_heat_flux_load': 126.95, 'total_strike_events': 31104.0, +# 'unique_struck_faces': 10368, 'swept_affected_area': 64.0, +# 'component_area': 64.0} +``` + +The envelope needs the pipeline's cumulative arrays, which only the full +strike path produces: with `output.vtk.enabled: false` it is reported as +empty rather than reconstructed from the per-firing records. + +`n_firings` keeps one meaning in both modes — **entries contributed by each +pose** — so a configuration switches modes by changing one line. In +`single_jfh` mode the single history is required to hold exactly +`len(poses) x n_firings` entries, checked after writing and again after +reading the file back. + +--- + +## YAML configuration + +A study configuration is a **layer on top of an existing PyRPOD case**. The +case's `config.ini` keeps owning the STL, TCF, TDF, plume model and gating +geometry, so every existing case and public API is unaffected; the YAML adds +only what a trade study needs and an INI cannot express. + +```yaml +study: + name: cai2016_flat_plate_baseline # required; identifies every result row + description: ... + case_dir: .. # required; relative to THIS file + output_dir: ../results/studies/flat_plate_baseline + +thruster: + id: T1 # optional; validated against the case TCF + +plume_model: + name: SimplifiedGasKinetics # the only accepted value + parameters: # recorded for provenance only + gas: argon + speed_ratio_S0: 2.0 + +target: + geometry_id: flat_plate_transformed.stl # defaults to the case's [tv] stl + reference_point: [0.0, 0.0, 0.0] # sweep is built about this point + normal: [0.0, 0.0, 1.0] # toward the plume-source side + tangent: [1.0, 0.0, 0.0] # sweep plane's in-plane axis + components: + - name: plate + selector: all # or face_indices: [...] / bounds: {...} + +sweep: + mode: per_case # per_case (default) | single_jfh + plate_angles_deg: [0.0] # 0 = head-on along `normal` + source_distances: [4.0] # from `reference_point` + n_firings: 1 # EXACT JFH entries per pose + firing_duration_s: 1.0 + thrusters: [1] + # firings: [...] # optional explicit poses; see below + +loads: + moment_reference_point: [0.0, 0.0, 0.0] + normalization: # optional; omit for no coefficients + reference_area: 64.0 + reference_length: 4.0 + dynamic_pressure: 1.1044652197738332 + reference_heat_flux: 637.3520362956127 + +output: + vtk: {enabled: true} + summary: {csv: case_results.csv, metadata: study_metadata.json} + plots: {enabled: false} + +reference: + path: null # optional external reference data + label: null + +metadata: + coordinate_system: case global frame + units: {} # overrides the SI defaults +``` + +Validation is strict and specific: a missing case directory, a case without a +`config.ini`, an unsupported plume model, an empty or non-positive sweep +axis, non-orthogonal target axes, duplicate component names, a non-positive +normalization value, or a firing list whose length disagrees with +`n_firings` each raise `StudyConfigError` naming the offending key. + +Angles and distances are geometry, not physics: `normal` and `tangent` define +the plane the source is swept in, so a curved target simply supplies the axes +its sweep should use (see [Adding another validation +geometry](#adding-another-validation-geometry)). + +--- + +## `n_firings` and prescribed firings + +**`n_firings` is the exact number of Jet Firing History entries contributed +by each swept pose.** It is validated before anything is generated (positive +integer only — zero, negative, fractional and non-numeric values are +rejected), the generated sequence is asserted to match it, and the JFH read +back from disk is asserted to report the same count. In `per_case` mode that +is the exact length of every case's history; in `single_jfh` mode the one +shared history must hold exactly `len(poses) * n_firings` entries. + +Two ways to supply the poses: + +- **Generated** (the usual case): the pose comes from the swept + `plate_angles_deg` / `source_distances` entry and is repeated for + `n_firings` successive firing intervals — one JFH entry per firing, each of + `firing_duration_s`. +- **Explicitly prescribed**: list the poses under `sweep.firings`. Their + count must agree with `n_firings` (in `single_jfh` mode, with + `len(poses) * n_firings`, the sequence being assigned to the poses in sweep + order); a mismatch is an error, never a silent truncation or extension. + +```yaml +sweep: + n_firings: 2 + firings: + - position: [0.0, 0.0, 4.0] + dcm: [[0, 0, 1], [0, 1, 0], [-1, 0, 0]] # first COLUMN = thruster axis + thrusters: [1] + duration_s: 0.5 + - position: [0.5, 0.0, 4.0] + dcm: [[0, 0, 1], [0, 1, 0], [-1, 0, 0]] +``` + +The generated pose convention, for a target reference point `C` with outward +normal `n_hat` and in-plane tangent `t_hat`: + +``` +d_hat(alpha) = cos(alpha) * n_hat + sin(alpha) * t_hat +source position = C + L * d_hat(alpha) +thruster axis = -d_hat(alpha) (aimed at C) +``` + +`alpha = 0` is head-on. The JFH DCM carries the thruster axis as its first +**column**, which is what the strike pipeline reads as the plume normal. This +reproduces the committed sweep-JFH generators +(`case/plume/plume_flat_plate_sweep/jfh/generate_sweep_jfh.py`) to file +precision — a regression test pins all 95 poses of that file. + +The dynamics-driven path was corrected to the same meaning: +`PlumeStrikeEstimationStudy.print_jfh_1d_approach_n_fire(..., n_firings=N)` +previously accepted `N` and ignored it. It now validates the count, stops the +approach simulation at exactly `N` entries, and raises if the approach +completes in fewer firings than requested. + +--- + +## Integrated loads + +`pyrpod/mdao/surface_loads.py` turns the strike pipeline's per-face arrays +into component-level resultants. Sign and direction conventions follow the +pipeline's own face-selection rule: a face is struck only when its unit +normal `n_hat` opposes the plume direction, i.e. points back toward the +source. + +### Force + +``` +F_p,i = -p_i * A_i * n_hat_i (pressure, into the surface) +t_hat_i = normalize(u_hat_i - (u_hat_i . n_in_i) * n_in_i), n_in_i = -n_hat_i +F_s,i = tau_i * A_i * t_hat_i (shear, along the tangential flow) +``` + +`u_hat_i` is the local flow direction — the unit vector from the plume source +to the face centroid, since the collisionless flow is radial. Faces at exactly +normal incidence have no tangential direction and contribute no shear force. +Reported: `pressure_force`, `shear_force`, `force`, `force_magnitude`. + +### Moment + +``` +M_ref = sum (r_i - r_ref) x F_i +``` + +about the user-defined `loads.moment_reference_point`. The pressure and shear +contributions are reported separately (`pressure_moment`, `shear_moment`) +alongside the combined `moment` and `moment_magnitude`. + +### Center of pressure + +A unique three-dimensional center of pressure does not exist in general: the +resultant of a distributed load is a force **plus a couple**, and only the +moment component perpendicular to the force can be represented by shifting +the force's line of action. The definition used here is stated explicitly: + +``` +r_cop = r_ref + (F x M_ref) / |F|^2 +``` + +— the point **on the line of action of the resultant force that lies closest +to the moment reference point**. The force-parallel moment component is +reported separately as `residual_couple`; no choice of `r_cop` can remove it. + +Degenerate cases are handled explicitly and never return a misleading value: + +| Situation | `center_of_pressure` | `center_of_pressure_status` | +| --- | --- | --- | +| Non-zero resultant | the point above | `ok` | +| Nothing loaded | `None` | `zero_load` | +| Resultant is cancellation noise (below `1e-6` of the summed face-force magnitudes, e.g. a closed target loaded symmetrically) | `None` | `ill_conditioned` | + +`pressure_weighted_centroid` — `sum(p_i A_i r_i) / sum(p_i A_i)` — is +reported as an auxiliary, always-defined location; it coincides with the +classical center of pressure for a planar component under unidirectional +pressure. + +### Thermal quantities + +Instantaneous and peak quantities only (time-integrated heat dose is +deliberately out of scope on this branch): `max_heat_flux` (peak per-face +heat flux, W/m^2) and `total_heat_load` (`sum q_i A_i` over the component, W), +plus `max_pressure`, `max_shear_stress`, `affected_area` and `struck_faces`. +The per-face heat-flux load the pipeline already accumulates stays in the VTK +output. + +### Coefficients + +Coefficients are computed **only** when the configuration supplies every +input a given coefficient needs. Nothing is defaulted or invented. + +| Coefficients | Required inputs | +| --- | --- | +| `Cp_max`, `Cf_max` | `dynamic_pressure` | +| `Cq_max` | `reference_heat_flux` | +| `CF`, `CFx`, `CFy`, `CFz` | `dynamic_pressure`, `reference_area` | +| `CM`, `CMx`, `CMy`, `CMz` | + `reference_length` | + +With no `normalization` block at all, `coefficients` is empty and +`coefficients_available` is `False` — which is what the cylinder example +does, since no reference values exist for it yet. + +--- + +## Result schema + +`StudyResults` holds one `CaseResult` per **case x component x firing**, each +carrying enough metadata to reproduce and later compare the calculation. In +`per_case` mode there is one `case_id` per pose; in `single_jfh` mode every +record shares the study's single `case_id` and is distinguished by +`firing_id`, but both carry the same `plate_angle_deg` / `source_distance`, +so downstream code (CSV, plots, reference comparison) is mode-agnostic: + +- **identity** — `study_name`, `case_id`, `component`, `firing_id` +- **geometry** — `geometry_id`, `mesh_faces`, `component_faces`, + `component_area`, `coordinate_system`, `units` +- **pose and sweep** — `plume_source_position`, `plume_source_orientation` + (9 DCM values), `target_normal`, `target_tangent`, + `target_reference_point`, `plate_angle_deg`, `source_distance`, + `firing_duration_s`, `thrusters` +- **model** — `plume_model`, `plume_model_parameters` +- **loads** — `pressure_force`, `shear_force`, `force`, `force_magnitude`, + `moment_reference_point`, `pressure_moment`, `shear_moment`, `moment`, + `moment_magnitude`, `center_of_pressure`, `center_of_pressure_status`, + `residual_couple`, `pressure_weighted_centroid` +- **surface fields** — `max_pressure`, `max_shear_stress`, `max_heat_flux`, + `total_heat_load`, `affected_area`, `struck_faces` +- **coefficients** — `coefficients`, `coefficients_available` +- **artifacts and provenance** — `vtk_path`, `jfh_path`, `config_path`, + `case_dir`, `code_version` (git commit when available), `generated_at` + +Two machine-readable artifacts are written, in formats the repository already +uses (no Parquet, no new dependency): + +- **CSV** (`StudyResults.write_csv`) — one flat row per record; vectors are + expanded to `_x/_y/_z` columns and each coefficient gets its own + `coeff_` column, so the file is directly plottable; +- **JSON** (`StudyResults.write_metadata`) — `schema`, study-level + `provenance` (including the plume model, mesh size, component list, code + version and known limitations) and the nested per-case records. This is the + document an externally generated dataset is later transformed into for + comparison. + +--- + +## VTK outputs + +Per-face data are never reduced away. With `output.vtk.enabled` (the default) +the standard pipeline writer produces one `.vtu` per firing, laid out +according to the sweep mode: + +- `per_case`: `/cases//results/strikes/firing-.vtu` +- `single_jfh`: `/results/strikes/firing-.vtu` — one numbered + series over the whole sweep, which ParaView opens as a time sequence + +Both carry the pipeline's own cell fields: + +``` +strikes, cum_strikes, pressures, max_pressures, shear_stress, max_shears, +heat_flux_rate, heat_flux_load, cum_heat_flux_load +``` + +Every case result advertises its own `vtk_path`, and that path is also in the +CSV and JSON summaries, so artifacts are discoverable from the returned +results. Set `output.vtk.enabled: false` for a fast, artifact-free run (the +numbers are identical; only the files are skipped). + +--- + +## Optional plots + +Plot generation is entirely optional — no automated test requires graphical +output — and matplotlib's non-interactive `Agg` backend is pinned when +`pyrpod.mdao.study_plots` is imported (lazily, only when plots are asked +for). Enable with `output.plots.enabled: true` or call `study.plot()`: + +``` +force_vs_angle.png moment_vs_angle.png heat_flux_vs_angle.png +force_vs_distance.png moment_vs_distance.png center_of_pressure.png +reference_comparison.png (when a comparison report exists) +``` + +--- + +## External reference data + +`pyrpod/mdao/reference_data.py` compares study results against independently +generated data **without knowing where it came from**. A reference record is +just named quantities attached to matching keys; DSMC, an analytical +solution, an experiment and another code are all handled identically, and no +producer-specific importer exists. + +Supported formats: + +- **CSV** — one row per record. `case_id`, `component`, `plate_angle_deg` and + `source_distance` are recognized as keys; every other numeric column is a + quantity, and `_x/_y/_z` triplets are folded into vectors. This is + exactly the layout `StudyResults.write_csv` emits, so an external producer + can be transformed into it column-for-column. +- **JSON / YAML**: + +```yaml +label: DSMC run 12 +source: /runs/dsmc-12 +units: {force: N} +records: + - key: {plate_angle_deg: -40.0, source_distance: 2.0, component: plate} + quantities: + force: [1.117, 0.0, -2.376] + max_pressure: 0.34 + center_of_pressure: [0.031, 0.0, 0.015] +``` + +```python +from pyrpod.mdao.reference_data import load_reference_dataset + +report = study.compare(load_reference_dataset('dsmc_run_12.json')) +print(report.max_relative_error('force')) +report.write_csv('comparison.csv') +``` + +Metrics are applied only where they are mathematically meaningful: absolute +error, relative error (`None` when the reference is zero, never infinity), +normalized RMSE (by reference range or mean), peak-value error, +integrated-load error (compared as vectors, so a load of the right magnitude +pointing the wrong way is an error) and center-of-pressure displacement. A +quantity the reference supplies but the result does not is reported as +`missing_candidate`; a case with no matching record is listed in +`unmatched_cases`. Nothing is fabricated or defaulted. + +For the flat-plate tests the independent reference is the exact Cai 2016 +solution (`pyrpod/plume/CaiImpingement2016.py`, Eq. 15 quadrature), +converted to dimensional loads with the case's own normalization — PyRPOD's +own output is never used as the validation baseline. Measured agreement over +the sweep: integrated normal load within ~1.5% mean / 2.2% max, component +heat load within ~15%, matching the documented accuracy of the Maxwellian +engineering chain against the exact collisionless solution. + +### Adding DSMC results later + +Nothing DSMC-specific is needed. Transform the external results into either +supported format, keyed by `plate_angle_deg` / `source_distance` / +`component` (or `case_id`, which the study metadata records), then call +`study.compare(...)` or point `reference.path` at the file. Because the +result schema records the mesh, coordinate system, units, source pose, model +and configuration provenance, the transformation can be done from the stored +metadata without rerunning PyRPOD. + +--- + +## Adding another validation geometry + +1. Create (or reuse) a PyRPOD case directory with its `config.ini`, target + STL, TCF and TDF — the study layer adds nothing here. +2. Add a study YAML next to it under `study/`, setting `target.reference_point` + to the geometry's reference (plate center, cylinder centroid, ...) and + `target.normal` / `target.tangent` to the plane the source should be swept + in. +3. Break the target into components if the loads should be reported + separately, using `face_indices` or a `bounds` box on face centroids. +4. Supply `loads.normalization` only if valid reference values exist for that + geometry; otherwise leave it out and the coefficients are reported as + unavailable. +5. Point `reference.path` at external data whenever it exists. + +`case/plume/plume_cylinder_sweep/study/cylinder_baseline.yaml` is a worked +example on a curved, closed target. Its results are a pipeline smoke case, +not a validated physical answer — no cylinder reference data exists yet — and +it deliberately supplies no normalization inputs. + +--- + +## Known limitations + +- **No shadowing or occlusion.** Face selection is the existing pipeline + behavior: a face is struck when it lies inside the plume wedge/radius and + its normal faces the source. Plume shadowing, occlusion by intervening + geometry, self-shadowing of concave or closed targets and back-facing + surfaces are **not** modeled, so a face hidden behind other geometry still + receives load. This branch changes nothing about that; results on closed + targets must be read with it in mind. +- **Edge-on poses are degenerate.** At exactly `+/-90 deg` the facing test's + dot product is zero in exact arithmetic, so float32 STL normals decide + strike membership arbitrarily. The plate-averaged load there is ~zero, but + individual near-nozzle faces can carry large grazing values. The committed + sweep keeps those poses; the symmetry checks exclude them. +- **Cumulative fields mean different things per mode.** In `single_jfh` mode + they are a genuine sweep envelope; in `per_case` mode they only restate the + case they belong to, since each case's history starts fresh. Read + `max_pressures` / `cum_strikes` in a per-case VTK accordingly. +- **Coefficients require explicit normalization.** By design: nothing is + inferred from the geometry. +- **One thruster, one plume model.** The study workflow prescribes a single + firing source and `SimplifiedGasKinetics`. Multi-thruster and multi-group + behavior is untouched elsewhere in PyRPOD but is not exercised here. +- **Cylinder validation is not quantitative yet.** The architecture and the + comparison interface are ready for it; the reference data is not. +- **Heat loading is instantaneous/peak only.** Time-integrated heat dose and + impulse are deliberately out of scope on this branch.