diff --git a/experiments/REDUCED_SURFACE_STUDY_2026-09-18.md b/experiments/REDUCED_SURFACE_STUDY_2026-09-18.md new file mode 100644 index 0000000..3a6daee --- /dev/null +++ b/experiments/REDUCED_SURFACE_STUDY_2026-09-18.md @@ -0,0 +1,64 @@ +# One-hour reduced smooth surface study — 2026-09-18 + +Authorized wall-clock window: 19:30:15–20:30:15 UTC. Branch: +research/reduced-smooth-surface-20260918, parent 6ceda61. + +Keep Q3 displacement/surface, Q4 reference geometry, Q2 DG density/potential, +Q3 enthalpy and RT2 gravity. No central-density constraint. n=1.5,J=0. +No production API, operator or extension change. Single MPI rank. + +Frozen affine embedding q=q0+B a. The base q0 is never projected. Construct +real spherical harmonics through degree L at actual surface radial nodes; +orthonormalize the Q3 interpolants using q14 prepared physical surface area. +Thin QR of S B and R^-1 B supplies separate normalized state and dual test +maps. Solve Qf^T Jnormalized Qx p=-Qf^T Fnormalized. The projected full inverse +Qx^T Mnormalized^-1 Qf is an approximation, not an exact reduced inverse. +All other unknowns/equations remain active. Cold FGMRES, eta=.03, max200, +restart160, same gravity-Schur preconditioner (one sweep, relaxation.5). + +Priority: fresh full-space direction at the prior Q3 six-update endpoint, +then L4 there, then L2 and L8 (L6 optional); analytic-seed L4 only if time +permits. Stop launching FEM cases by20:20; reserve remaining time for independent +artifact audit and reporting. No convergence claim from frozen steps. Cases, +builds, and full numerical audits run serially. A library-only Q3 warm rebuild +may overlap editing experiment-only headers; no production source changes. + +Raw independent trial alphas 1,1/4,1/16,1/1024. Preserve original and reduced +residuals, separate linear gates, full/reduced slopes, excluded shape forces, +and exact restoration checks. Observational merit criterion is +norm(Ftrial)<=(1-1e-4 alpha)*norm(Fbase), evaluated separately in both spaces. +No trial accepted, no SOC, no positivity repair. q14/q18 physical measurements +at base and the first/largest prepared trial only; all trial residuals saved. +Direction physical observer q10/q14 reused from the preceding audited study. + +Acceptance of an observation requires algebra/transpose/full-rank control tests, +true reduced gate classification, finite vectors, source/input/output hashes, +and exact saved-state restoration. Reduced success is not full success. +Virial1e-5, interior physical errors, density sign, and omitted force remain +separate measurements. No selected norm reduction can certify a physical root. + +## Implementation + +Experiment-only headers: +- `reduced_surface_basis.hpp`: scalar radial-amplitude FE trace basis and + quadrature-defined physical surface mass metric, explicit rank/orthogonality + checks, and raw data for independent verification. +- `reduced_surface_operators.hpp`: two separate QR maps, reduced operator, + projected flexible preconditioner and deterministic algebra checks. +- `reduced_surface_diagnostics.hpp`: frozen solve, both original/reduced gates, + excluded-force decomposition, raw trials and restoration invariants. +- `audit_reduced_surface.py`: independent artifact algebra and corruption tests. + +The diagnostic executable adds experiment-only flags --reduced-surface-probe +and --surface-degree (0 means full-space control); no production-facing API is +modified. Degrees are runtime experimental choices while operator/inverse +compatibility and the inherited preconditioner application contract remain +compile-time constrained. Actual mesh dimensions and normalization factors are +necessarily checked at runtime. Existing physical observers are reused. + +Evidence root: `../../diagnostics/reduced_surface_2026-09-18` relative to this +checkout. The preceding study remains on its two committed research branches. + +## Results + +Pending the bounded serial campaign and its independent artifact audit. diff --git a/experiments/audit_reduced_surface.py b/experiments/audit_reduced_surface.py new file mode 100644 index 0000000..c8c6c6b --- /dev/null +++ b/experiments/audit_reduced_surface.py @@ -0,0 +1,517 @@ +#!/usr/bin/env python3 +"""Audit saved reduced-surface algebra; a reduced gate is not a full solve. + +No FEM action, preparation, or physical quadrature is independently repeated. +The dual residual restriction is checked against the physical basis and saved +normalization, rather than assumed to be the state injection's transpose. +""" +import argparse +import json +import math +from pathlib import Path +import sys + +from audit_shape_field_weak import (Audit, EPS, RTOL, checked_einsum as einsum, + digest, finite, integer, layout, near, norm, + np, rows, table, vector, vector_near) +from analyze_shape_surface_modes import harmonics + +P = "reduced_surface_" +MAPS = ("basis_physical", "state_basis", "residual_basis", "state_transform", "residual_transform") +VECTORS = ("base", "physical_base", "residual", "direction", "physical_direction", "action", "defect", + "state_factors", "residual_factors", "reduced_residual", "reduced_direction", "reduced_action", + "reduced_defect", "excluded_residual", "excluded_defect", "restored_residual", "restored_action") +ALPHAS = (1., .25, .0625, 1 / 1024) + + +def matrix(path): + with Path(path).open() as stream: + dimensions = [int(x) for x in stream.readline().split()] + values = np.loadtxt(stream, ndmin=1) + if len(dimensions) != 2 or min(dimensions) <= 0 or values.ndim != 1 or values.size != math.prod(dimensions) or not np.isfinite(values).all(): + raise ValueError(f"Invalid finite row-major matrix: {path}") + return values.reshape(dimensions) + + +def close_array(audit, label, actual, expected, scale=0.): + audit.check(label, vector_near(np.asarray(actual).ravel(), np.asarray(expected).ravel(), scale)) + + +class Coordinates: + def __init__(self, summary, maps): + self.n = integer(summary["full_size"]) + self.s = integer(summary["shape_size"]) + self.k = integer(summary["mode_count"]) + self.vo = integer(summary["shape_value_offset"]) + self.ro = integer(summary["shape_residual_offset"]) + self.identity = summary["identity_control"] == 1 + self.maps = maps + if not (0 < self.k <= self.s <= self.n and 0 <= self.vo <= self.n - self.s and 0 <= self.ro <= self.n - self.s): + raise ValueError("Invalid coordinate dimensions") + + def restrict(self, values, residual=True): + if self.identity: + return values.copy() + offset = self.ro if residual else self.vo + basis = self.maps["residual_basis" if residual else "state_basis"] + middle = einsum("ij,i->j", basis, values[offset:offset + self.s]) + return np.concatenate((values[:offset], middle, values[offset + self.s:])) + + def inject(self, values, residual=False): + if self.identity: + return values.copy() + offset = self.ro if residual else self.vo + basis = self.maps["residual_basis" if residual else "state_basis"] + middle = einsum("ij,j->i", basis, values[offset:offset + self.k]) + return np.concatenate((values[:offset], middle, values[offset + self.k:])) + + def excluded(self, values): + return values - self.inject(self.restrict(values), residual=True) + + def validate(self, factors, residual_factors, audit): + if self.identity: + audit.check("identity:full_mode_count", self.k == self.s) + return + b, px, qf, tx, tf = (self.maps[k] for k in MAPS) + audit.check("maps:dimensions", b.shape == px.shape == qf.shape == (self.s, self.k) and tx.shape == tf.shape == (self.k, self.k)) + sx = factors[self.vo:self.vo + self.s] + rf = residual_factors[self.ro:self.ro + self.s] + close_array(audit, "maps:state_QR", sx[:, None] * b, einsum("ij,jk->ik", px, tx)) + close_array(audit, "maps:dual_QR", b / rf[:, None], einsum("ij,jk->ik", qf, tf)) + for name, basis, transform in (("state", px, tx), ("residual", qf, tf)): + close_array(audit, "maps:" + name + "_orthonormal", einsum("ij,ik->jk", basis, basis), np.eye(self.k)) + close_array(audit, "maps:" + name + "_triangular", np.tril(transform, -1), np.zeros((self.k, self.k)), norm(transform.ravel())) + audit.check("maps:" + name + "_nonsingular", np.all(np.abs(np.diag(transform)) > 64 * EPS * max(norm(transform.ravel()), 1e-300))) + + +def basis_package(folder, coords, physical_base, required, audit): + first = folder / (P + "basis_raw.txt") + if not first.exists(): + raise ValueError("Missing physical harmonic basis package") + def get(name): + path = folder / (P + name + ".txt") + required.add(path.name) + return matrix(path) + raw, gram, orthogram, transform, samples = (get(k) for k in + ("basis_raw", "gram_raw", "gram_orthonormal", "coefficient_transform", "basis_samples_raw")) + weights_path = folder / (P + "basis_weights.txt") + required.add(weights_path.name) + weights = vector(weights_path) + audit.check("basis:positive_quadrature", len(weights) == len(samples) and np.all(weights > 0)) + close_array(audit, "basis:quadrature_gram", gram, einsum("ni,n,nj->ij", samples, weights, samples)) + close_array(audit, "basis:physical_from_raw", coords.maps["basis_physical"], einsum("ij,jk->ik", raw, transform)) + actual = einsum("ij,ik,kl->jl", transform, gram, transform) + close_array(audit, "basis:transformed_gram", orthogram, actual) + close_array(audit, "basis:mass_orthonormal", actual, np.eye(coords.k)) + node_path, mode_path, meta_path = (folder / (P + n + ".csv") for n in ("basis_nodes", "basis_modes", "basis_metadata")) + required.update(p.name for p in (node_path, mode_path, meta_path)) + nodes, modes = rows(node_path), rows(mode_path) + audit.check("basis:node_indices", [integer(r["parameter"]) for r in nodes] == list(range(coords.s))) + node_base = np.array([finite(r["base_parameter"]) for r in nodes]) + reference_radius = np.array([finite(r["reference_radius"]) for r in nodes]) + physical_radius = np.array([finite(r["physical_radius"]) for r in nodes]) + close_array(audit, "basis:frozen_base_parameters", node_base, physical_base[coords.vo:coords.vo + coords.s]) + close_array(audit, "basis:physical_radius", physical_radius, reference_radius + node_base) + audit.check("basis:positive_radii", np.all(reference_radius > 0) and np.all(physical_radius > 0)) + points = np.array([[finite(r["unit_radial_" + axis]) for axis in "xyz"] for r in nodes]) + lengths = np.sqrt(einsum("ij,ij->i", points, points)) + close_array(audit, "basis:unit_directions", lengths, np.ones(coords.s)) + maximum_degree = max(integer(r["degree"]) for r in modes) + expected, labels = harmonics(points / lengths[:, None], maximum_degree) + audit.check("basis:mode_indices", [integer(r["column"]) for r in modes] == list(range(coords.k))) + audit.check("basis:mode_order", [(integer(r["degree"]), integer(r["m"])) for r in modes] == [(l["degree"], l["m"]) for l in labels]) + close_array(audit, "basis:independent_harmonics", raw, expected) + text_keys = {"metric_scope", "scalar_finite_element_collection", "displacement_finite_element_collection"} + metadata = {k: r["value"] if k in text_keys else finite(r["value"]) + for k, r in table(meta_path, lambda r: r["metric"]).items()} + for key, value in {"maximum_degree": maximum_degree, "modes": coords.k, "surface_parameters": coords.s, + "quadrature": 14, "surface_quadrature_points": len(weights), "accepted_buffers_unchanged": 1, + "prepared_state_unchanged": 1, "normalization_factors_unchanged": 1, + "vector_displacement_metric_claimed": 0, "base_projection_performed": 0, + "trial_preparations": 0, "radial_trace_validation_actions": 1}.items(): + audit.check("basis:metadata:" + key, metadata[key] == value) + audit.check("basis:metric_scope", metadata["metric_scope"] == "scalar_radial_amplitude_FE_trace_mass_on_current_physical_surface") + audit.check("basis:reported_error_bounds", 0 <= metadata["maximum_gram_orthogonality_error"] <= 1e-10 and + 0 <= metadata["maximum_sample_reconstruction_error"] <= 1e-10 and 0 <= metadata["adjoint_relative_error"] <= 1e-12) + for j, row in enumerate(modes): + pivot, relative = finite(row["pivot_norm"]), finite(row["relative_pivot"]) + audit.check(f"basis:independent_mode:{j}", pivot > 0 and relative > metadata["rank_relative_pivot_threshold"] > 0) + audit.close(f"basis:pivot_transform:{j}", pivot * transform[j, j], 1.) + audit.close(f"basis:relative_pivot:{j}", relative, pivot / math.sqrt(gram[j, j])) + audit.close("basis:surface_area", metadata["physical_surface_area"], math.fsum(weights)) + return {"status": "nodal_basis_and_saved_quadrature_gram_checked", "metadata": metadata, + "scope": "Saved FE trace samples and weights are checked algebraically; physical quadrature is not rerun."} + + +def provenance(folder, required, audit): + run_path = folder.parent / (folder.name + "_run.json") + source_path = folder.parent / "provenance.json" + if not run_path.exists() and not source_path.exists(): + return {"status": "unavailable", "bound": False} + if not run_path.exists() or not source_path.exists(): + raise ValueError("Incomplete campaign provenance") + run, source = json.loads(run_path.read_text()), json.loads(source_path.read_text()) + audit.check("provenance:required_hashes", required <= set(run["output_hashes"])) + for name, expected in run["output_hashes"].items(): + audit.check("provenance:output:" + name, digest(folder / name) == expected) + for name, expected in run.get("input_state_hashes", {}).items(): + audit.check("provenance:input:" + name, digest(name) == expected) + command = run.get("command", run.get("cmd")) + if command is None: + raise ValueError("Missing recorded command") + case = next(c for c in source["cases"] if c["id"] == folder.name) + expected = [str(folder.parent / "plateau_diagnostic"), *case["arguments"], "--mesh", source["mesh"], "--output", str(folder)] + audit.check("provenance:command", command == expected) + audit.check("provenance:binary", digest(command[0]) == source["binary_sha256"]) + audit.check("provenance:mesh", digest(source["mesh"]) == source["mesh_sha256"]) + if "--probe-state" in command: + path = Path(command[command.index("--probe-state") + 1]) + if not path.is_absolute(): + path = Path(source["source"]) / path + audit.check("provenance:base_input_bound", str(path.resolve()) in run.get("input_state_hashes", {})) + close_array(audit, "provenance:physical_base_input", vector(path), vector(folder / (P + "physical_base.txt"))) + audit.check("provenance:successful_probe_exit", run["exit_code"] == 0) + return {"status": "runner_hashes_checked", "bound": True, "source_commit": source["commit"], + "run_sha256": digest(run_path), "provenance_sha256": digest(source_path)} + + +def inspect(folder): + folder = Path(folder).resolve() + audit, required, observations = Audit(), set(), {} + try: + def get(name): + path = folder / (P + name + ".txt") + required.add(path.name) + return vector(path) + summary_path = folder / (P + "summary.csv") + trials_path = folder / (P + "trials.csv") + layout_path = folder / (P + "layout.csv") + required.update(p.name for p in (summary_path, trials_path, layout_path)) + summary = {k: finite(r["value"]) for k, r in table(summary_path, lambda r: r["metric"]).items()} + audit.check("identity:boolean", summary["identity_control"] in (0, 1)) + maps = {} + for name in MAPS: + path = folder / (P + name + ".txt") + if summary["identity_control"]: + audit.check("identity:matrix_absent:" + name, not path.exists()) + else: + required.add(path.name) + maps[name] = matrix(path) + coordinates = Coordinates(summary, maps) + n, s, k, vo, ro = coordinates.n, coordinates.s, coordinates.k, coordinates.vo, coordinates.ro + degree = integer(summary["maximum_degree"]) + audit.check("prescribed_basis_size", (degree == 0 and k == s) if coordinates.identity else (degree in (2, 4, 6, 8) and k == (degree + 1) ** 2)) + audit.check("reduced_size", summary["reduced_size"] == n - s + k) + audit.check("schema_version", summary["schema_version"] == 1) + audit.check("reported_algebra_selfcheck", summary["algebra_self_check_passed"] == 1) + blocks = layout(layout_path, n) + audit.check("layout:value_shape", ("surface_deformation", vo, s) in blocks["value"]) + audit.check("layout:residual_shape", ("surface_shape_balance", ro, s) in blocks["residual"]) + data = {name: get(name) for name in VECTORS} + for name, values in data.items(): + audit.check("size:" + name, len(values) == (n - s + k if name.startswith("reduced_") else n)) + sx, rf = data["state_factors"], data["residual_factors"] + audit.check("normalization:positive", np.all(sx > 0) and np.all(rf > 0)) + coordinates.validate(sx, rf, audit) + for key, value in {"accepted_buffers_bitwise_identical": 1, "normalization_factors_bitwise_identical": 1, + "direction_buffers_bitwise_identical": 1, "restored_residual_error": 0, + "restored_action_error": 0, "probe_exception": 0, "accepted_commits": 0}.items(): + audit.check("invariant:" + key, summary[key] == value) + audit.check("restoration:residual_exact", np.array_equal(data["residual"], data["restored_residual"])) + audit.check("restoration:action_exact", np.array_equal(data["action"], data["restored_action"])) + close_array(audit, "normalization:base", data["physical_base"], data["base"] / sx) + close_array(audit, "normalization:direction", data["physical_direction"], data["direction"] / sx) + close_array(audit, "injection:direction", data["direction"], coordinates.inject(data["reduced_direction"])) + for name in ("residual", "action", "defect"): + # Forming R(F+Jp) and RF+RJp rounds differently when F and Jp + # cancel. The allowance must scale with the operands, not only + # with the potentially zero defect being checked. + scale = norm(data["residual"]) + norm(data["action"]) if name == "defect" else 0. + close_array(audit, "restriction:" + name, data["reduced_" + name], coordinates.restrict(data[name]), scale) + close_array(audit, "defect:full", data["defect"], data["residual"] + data["action"], norm(data["residual"]) + norm(data["action"])) + close_array(audit, "defect:reduced", data["reduced_defect"], data["reduced_residual"] + data["reduced_action"], + norm(data["reduced_residual"]) + norm(data["reduced_action"])) + for name in ("residual", "defect"): + excluded = coordinates.excluded(data[name]) + close_array(audit, "excluded:" + name, data["excluded_" + name], excluded, norm(data[name])) + audit.close("pythagoras:" + name, norm(data[name]) ** 2, + norm(data["reduced_" + name]) ** 2 + norm(excluded) ** 2, norm(data[name]) ** 2) + eta = summary["linear_relative_tolerance"] + absolute = summary["linear_absolute_tolerance"] + audit.check("linear:controls", 0 < eta < 1 and absolute >= 0 and integer(summary["linear_maximum_iterations"]) > 0) + gates, slopes = {}, {} + for label, prefix in (("full", ""), ("reduced", "reduced_")): + fn, dn = norm(data[prefix + "residual"]), norm(data[prefix + "defect"]) + audit.check(label + ":nonzero_rhs", fn > 0) + passed = dn <= max(absolute, eta * fn) + audit.check(label + ":reported_gate", summary[label + "_verified_gate"] in (0, 1) and bool(summary[label + "_verified_gate"]) == passed) + audit.close(label + ":relative_defect", summary[label + "_relative_defect"], dn / fn) + gates[label] = {"verified_gate": passed, "residual_norm": fn, "defect_norm": dn, "relative_defect": dn / fn} + slopes[label] = float(einsum("i,i->", data[prefix + "residual"], data[prefix + "action"])) + audit.close(label + ":reported_norm", summary[label + "_base_norm"], fn) + audit.close(label + ":reported_slope", summary[label + "_merit_slope"], slopes[label], fn * norm(data[prefix + "action"])) + audit.close("reported_direction_norm", summary["direction_norm"], norm(data["direction"])) + audit.close("reported_excluded_base", summary["excluded_base_norm"], norm(data["excluded_residual"]), norm(data["residual"])) + audit.close("linear:reported_true_norm", summary["linear_reported_true_norm"], norm(data["reduced_defect"]), norm(data["reduced_residual"])) + audit.check("linear:iteration_budget", 0 <= integer(summary["linear_iterations"]) <= integer(summary["linear_maximum_iterations"])) + trials = rows(trials_path) + audit.check("trial:complete_grid", len(trials) == len(ALPHAS)) + trial_observations = [] + for index, row in enumerate(trials): + alpha = finite(row["alpha"]) + audit.check(f"trial:{index}:alpha", index < len(ALPHAS) and alpha == ALPHAS[index]) + audit.check(f"trial:{index}:sample", integer(row["sample"]) == index) + prepared = integer(row["prepared"]) + audit.check(f"trial:{index}:prepared_flag", prepared in (0, 1)) + result = {"sample": index, "alpha": alpha, "prepared": bool(prepared), "reason": integer(row["reason"]), "stage": integer(row["stage"])} + candidate_path = folder / (P + f"candidate_a{index}.txt") + if prepared or candidate_path.exists(): + candidate = get(f"candidate_a{index}") + close_array(audit, f"trial:{index}:candidate", candidate, data["base"] + alpha * data["direction"]) + if prepared: + audit.check(f"trial:{index}:prepared_status", result["reason"] == result["stage"] == -1) + physical_candidate = get(f"physical_candidate_a{index}") + close_array(audit, f"trial:{index}:physical_candidate", physical_candidate, candidate / sx) + full, reduced = get(f"residual_a{index}"), get(f"reduced_residual_a{index}") + close_array(audit, f"trial:{index}:restriction", reduced, coordinates.restrict(full)) + excluded = coordinates.excluded(full) + full_norm, reduced_norm, excluded_norm = norm(full), norm(reduced), norm(excluded) + remainder = full - data["residual"] - alpha * data["action"] + reduced_remainder = reduced - data["reduced_residual"] - alpha * data["reduced_action"] + for label, actual in (("full_norm", full_norm), ("reduced_norm", reduced_norm), ("excluded_shape_norm", excluded_norm), + ("full_remainder_norm", norm(remainder)), ("reduced_remainder_norm", norm(reduced_remainder))): + audit.close(f"trial:{index}:" + label, finite(row[label]), actual, norm(data["residual"])) + audit.close(f"trial:{index}:pythagoras", full_norm ** 2, reduced_norm ** 2 + excluded_norm ** 2, full_norm ** 2) + for label, actual in (("full", full_norm), ("reduced", reduced_norm)): + accepted = actual <= (1 - 1e-4 * alpha) * gates[label]["residual_norm"] + audit.check(f"trial:{index}:" + label + "_accept", integer(row[label + "_merit_pass"]) in (0, 1) and bool(integer(row[label + "_merit_pass"])) == accepted) + result[label + "_merit_pass"] = accepted + result.update(full_norm=full_norm, reduced_norm=reduced_norm, excluded_shape_norm=excluded_norm, + full_remainder_norm=norm(remainder), reduced_remainder_norm=norm(reduced_remainder)) + else: + audit.check(f"trial:{index}:failure_status", result["reason"] >= 0 and result["stage"] >= 0) + for name in (f"residual_a{index}", f"reduced_residual_a{index}"): + audit.check(f"trial:{index}:unavailable:" + name, not (folder / (P + name + ".txt")).exists()) + for label in ("full_norm", "reduced_norm", "excluded_shape_norm", "full_remainder_norm", "reduced_remainder_norm"): + audit.check(f"trial:{index}:unavailable:" + label, math.isnan(float(row[label]))) + audit.check(f"trial:{index}:failure_not_accepted", integer(row["full_merit_pass"]) == integer(row["reduced_merit_pass"]) == 0) + trial_observations.append(result) + package = basis_package(folder, coordinates, data["physical_base"], required, audit) if not coordinates.identity else {"status": "identity_control"} + if not coordinates.identity: + audit.check("basis:degree_matches_summary", package["metadata"]["maximum_degree"] == degree) + binding = provenance(folder, required, audit) + observations = {"summary": summary, "linear": gates, "merit_directional_derivatives": slopes, + "classification": "full_and_reduced_gates_pass" if gates["full"]["verified_gate"] and gates["reduced"]["verified_gate"] else + "reduced_gate_only_not_full_Newton" if gates["reduced"]["verified_gate"] else "finite_reduced_solve_failed_gate", + "base_excluded_shape_norm": norm(data["excluded_residual"]), "trials": trial_observations, + "basis": package, "provenance": binding} + except (ValueError, OSError, KeyError, TypeError, IndexError, StopIteration, np.linalg.LinAlgError) as error: + audit.check("complete_readable_consistent_artifacts:" + str(error), False) + return {"directory": str(folder), "artifact_pass": not audit.failures, "checks_count": audit.count, + "failures": audit.failures, "required_artifacts": sorted(required), "observations": observations if not audit.failures else {}} + + +def self_test(): + count = 0 + def checked(condition): + nonlocal count + assert condition + count += 1 + # Deliberately different value/residual factors and shape offsets catch the + # tempting but invalid reuse of the state basis as a residual test basis. + b = np.array([[1., .4], [1., -.1], [1., -.3]]) + sx = np.array([2., 3., 5., 7., 11., 13., 17.]) + rf = np.array([.2, .3, .7, 1.1, 1.3, 1.7, 1.9]) + summary = dict(full_size=7, shape_size=3, mode_count=2, shape_value_offset=1, shape_residual_offset=3, identity_control=0) + px, tx = np.linalg.qr(sx[1:4, None] * b) + qf, tf = np.linalg.qr(b / rf[3:6, None]) + maps = dict(zip(MAPS, (b, px, qf, tx, tf))) + c = Coordinates(summary, maps) + a = Audit(); c.validate(sx, rf, a); checked(not a.failures) + z = np.arange(1., 7.) + f = np.arange(-3., 4.) + checked(vector_near(c.restrict(c.inject(z), False), z)) + checked(vector_near(c.restrict(c.inject(z, True)), z)) + checked(near(norm(f) ** 2, norm(c.restrict(f)) ** 2 + norm(c.excluded(f)) ** 2)) + physical_dual = einsum("ij,i->j", b, f[3:6] / rf[3:6]) + checked(vector_near(einsum("ji,j->i", tf, c.restrict(f)[3:5]), physical_dual)) + bad = Coordinates(summary, {**maps, "residual_basis": px}) + a = Audit(); bad.validate(sx, rf, a); checked("maps:dual_QR" in a.failures) + a = Audit(); c.validate(sx, rf[::-1], a); checked(bool(a.failures)) + checked(norm(c.restrict(c.excluded(f))) < 1e-14) + # Full-rank bases may rotate the coordinates. Verify the recovered full + # direction, not equality of the rotated coefficient arrays. + bfull = np.array([[1., .4, .2], [1., -.1, .7], [1., -.3, -.6]]) + px, tx = np.linalg.qr(sx[1:4, None] * bfull) + qf, tf = np.linalg.qr(bfull / rf[3:6, None]) + full = Coordinates({**summary, "mode_count": 3}, dict(zip(MAPS, (bfull, px, qf, tx, tf)))) + op = np.diag(np.arange(2., 9.)) + .03 * np.ones((7, 7)) + columns = [] + for i in range(7): + columns.append(full.restrict(einsum("ij,j->i", op, full.inject(np.eye(7)[:, i])))) + reduced_op = np.column_stack(columns) + result = full.inject(np.linalg.solve(reduced_op, -full.restrict(f))) + checked(vector_near(result, np.linalg.solve(op, -f))) + identity = Coordinates({**summary, "mode_count": 3, "identity_control": 1}, {}) + checked(np.array_equal(identity.inject(f), f) and np.array_equal(identity.restrict(f), f)) + checked(np.array_equal(identity.excluded(f), np.zeros(7))) + # A reduced zero defect can coexist with a substantial discarded full one. + discarded = c.excluded(f) + checked(norm(discarded) > .1 and norm(c.restrict(discarded)) < 1e-14) + import csv + import tempfile + def fixture(folder, identity=False): + folder.mkdir() + def save(name, array): + array = np.asarray(array) + with (folder / (P + name + ".txt")).open("w") as stream: + stream.write((f"{len(array)}" if array.ndim == 1 else f"{array.shape[0]} {array.shape[1]}") + "\n") + np.savetxt(stream, array.ravel(), fmt="%.17g") + def csv_save(name, columns, records): + with (folder / (P + name + ".csv")).open("w", newline="") as stream: + writer = csv.writer(stream); writer.writerow(columns); writer.writerows(records) + n, s, k, vo, ro = 18, 12, 9, 1, 4 + points = [] + for i in range(s): + z = 1 - 2 * (i + .5) / s + phi = i * math.pi * (3 - math.sqrt(5)) + points.append((math.sqrt(1 - z * z) * math.cos(phi), math.sqrt(1 - z * z) * math.sin(phi), z)) + points = np.asarray(points) + raw, modes = harmonics(points, 2) + weights = np.ones(s) + gram = einsum("ni,n,nj->ij", raw, weights, raw) + transform = np.linalg.inv(np.linalg.cholesky(gram).T) + b = einsum("ij,jk->ik", raw, transform) + sx, rf = np.linspace(1., 3., n), np.linspace(.4, 2., n) + px, tx = np.linalg.qr(sx[vo:vo + s, None] * b) + qf, tf = np.linalg.qr(b / rf[ro:ro + s, None]) + maps = dict(zip(MAPS, (b, px, qf, tx, tf))) + summary = dict(schema_version=1, full_size=n, shape_size=s, mode_count=s if identity else k, + shape_value_offset=vo, shape_residual_offset=ro, identity_control=int(identity), + maximum_degree=0 if identity else 2, reduced_size=n if identity else n - s + k, + linear_relative_tolerance=.03, linear_absolute_tolerance=0., linear_maximum_iterations=200, + linear_restart=160, linear_status=0, linear_iterations=7, linear_seconds=.01, + algebra_self_check_passed=1, accepted_commits=0, accepted_buffers_bitwise_identical=1, + normalization_factors_bitwise_identical=1, direction_buffers_bitwise_identical=1, + restored_residual_error=0, restored_action_error=0, preparations=8, probe_exception=0) + c = Coordinates(summary, {} if identity else maps) + reduced_direction = np.sin(np.arange(summary["reduced_size"]) + 1.) + direction = c.inject(reduced_direction) + action = direction.copy() + residual = -action + c.excluded(np.cos(np.arange(n))) + base = np.linspace(.1, .3, n) + data = dict(base=base, physical_base=base / sx, residual=residual, direction=direction, + physical_direction=direction / sx, action=action, defect=residual + action, + state_factors=sx, residual_factors=rf, reduced_direction=reduced_direction, + restored_residual=residual, restored_action=action) + for name in ("residual", "action", "defect"): + data["reduced_" + name] = c.restrict(data[name]) + for name in ("residual", "defect"): + data["excluded_" + name] = c.excluded(data[name]) + for label, prefix in (("full", ""), ("reduced", "reduced_")): + fn, dn = norm(data[prefix + "residual"]), norm(data[prefix + "defect"]) + summary[label + "_base_norm"] = fn + summary[label + "_relative_defect"] = dn / fn + summary[label + "_verified_gate"] = int(dn <= .03 * fn) + summary[label + "_merit_slope"] = float(einsum("i,i->", data[prefix + "residual"], data[prefix + "action"])) + summary["linear_reported_true_norm"] = norm(data["reduced_defect"]) + summary["direction_norm"] = norm(direction) + summary["excluded_base_norm"] = norm(data["excluded_residual"]) + for name, values in data.items(): + save(name, values) + csv_save("summary", ["metric", "value"], summary.items()) + csv_save("layout", ["kind", "block", "offset", "size"], + [("value", "density", 0, vo), ("value", "surface_deformation", vo, s), ("value", "tail", vo + s, n - vo - s), + ("residual", "head", 0, ro), ("residual", "surface_shape_balance", ro, s), ("residual", "tail", ro + s, n - ro - s)]) + trials = [] + for index, alpha in enumerate(ALPHAS): + save(f"candidate_a{index}", base + alpha * direction) + if index == 0: + trials.append((index, alpha, 0, 0, 0, "nan", "nan", "nan", 0, 0, "nan", "nan")) + continue + f = residual + alpha * action + .025 * alpha ** 2 * residual + fr = c.restrict(f) + remainder = f - residual - alpha * action + save(f"physical_candidate_a{index}", (base + alpha * direction) / sx) + save(f"residual_a{index}", f); save(f"reduced_residual_a{index}", fr) + trials.append((index, alpha, 1, -1, -1, norm(f), norm(fr), norm(c.excluded(f)), + int(norm(f) <= (1 - 1e-4 * alpha) * norm(residual)), + int(norm(fr) <= (1 - 1e-4 * alpha) * norm(data["reduced_residual"])), + norm(remainder), norm(c.restrict(remainder)))) + csv_save("trials", ["sample", "alpha", "prepared", "reason", "stage", "full_norm", "reduced_norm", "excluded_shape_norm", + "full_merit_pass", "reduced_merit_pass", "full_remainder_norm", "reduced_remainder_norm"], trials) + if not identity: + for name, values in maps.items(): + save(name, values) + for name, values in {"basis_raw": raw, "gram_raw": gram, "coefficient_transform": transform, + "gram_orthonormal": einsum("ij,ik,kl->jl", transform, gram, transform), + "basis_samples_raw": raw, "basis_weights": weights}.items(): + save(name, values) + csv_save("basis_nodes", ["parameter", "unit_radial_x", "unit_radial_y", "unit_radial_z", "reference_radius", "physical_radius", "base_parameter"], + [(i, *point, 1., 1 + data["physical_base"][vo + i], data["physical_base"][vo + i]) for i, point in enumerate(points)]) + csv_save("basis_modes", ["column", "degree", "m", "pivot_norm", "relative_pivot"], + [(i, m["degree"], m["m"], 1 / transform[i, i], 1 / (transform[i, i] * math.sqrt(gram[i, i]))) for i, m in enumerate(modes)]) + metadata = dict(maximum_degree=2, modes=k, surface_parameters=s, quadrature=14, surface_quadrature_points=s, + physical_surface_area=s, accepted_buffers_unchanged=1, prepared_state_unchanged=1, normalization_factors_unchanged=1, + vector_displacement_metric_claimed=0, base_projection_performed=0, trial_preparations=0, + radial_trace_validation_actions=1, metric_scope="scalar_radial_amplitude_FE_trace_mass_on_current_physical_surface", + maximum_gram_orthogonality_error=1e-15, maximum_sample_reconstruction_error=0., + adjoint_relative_error=0., rank_relative_pivot_threshold=1e-12) + csv_save("basis_metadata", ["metric", "value"], metadata.items()) + return save + with tempfile.TemporaryDirectory() as tmp: + folder = Path(tmp) / "reduced" + save = fixture(folder) + result = inspect(folder) + assert result["artifact_pass"], result["failures"] + checked(result["artifact_pass"]) + checked(result["observations"]["classification"] == "reduced_gate_only_not_full_Newton") + for name in ("residual_basis", "candidate_a1", "excluded_residual", "gram_raw", "state_factors"): + path = folder / (P + name + ".txt") + original = path.read_bytes() + values = matrix(path) if name in ("residual_basis", "gram_raw") else vector(path) + values.flat[0] += .02 + save(name, values) + checked(not inspect(folder)["artifact_pass"]) + path.write_bytes(original) + summary_path = folder / (P + "summary.csv") + original = summary_path.read_text() + summary_path.write_text(original.replace("full_verified_gate,0", "full_verified_gate,1")) + checked(not inspect(folder)["artifact_pass"]) + summary_path.write_text(original) + save("residual_a0", np.ones(18)) + checked(not inspect(folder)["artifact_pass"]) + (folder / (P + "residual_a0.txt")).unlink() + identity_folder = Path(tmp) / "identity" + fixture(identity_folder, True) + result = inspect(identity_folder) + checked(result["artifact_pass"] and result["observations"]["classification"] == "full_and_reduced_gates_pass") + print(f"{count} tiny algebra, end-to-end artifact and corruption checks passed") + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("directories", nargs="*") + parser.add_argument("--output", type=Path) + parser.add_argument("--self-test", action="store_true") + args = parser.parse_args() + if args.self_test: + self_test() + return 0 + if not args.directories or args.output is None: + parser.error("Provide probe directories and a fresh output JSON") + output = args.output.resolve() + if output.exists() or any(output.is_relative_to(Path(d).resolve()) for d in args.directories): + parser.error("Use a fresh output outside immutable probe directories") + cases = [inspect(d) for d in args.directories] + result = {"schema_version": 1, "passed": all(c["artifact_pass"] for c in cases), "checker_sha256": digest(__file__), + "tolerances": {"relative": RTOL, "cancellation_epsilon_multiplier": 2048}, "cases": cases, + "scope": "Independent saved-vector and map algebra only. Reduced success does not certify the unrestricted equations, FEM evaluation, or physical validity."} + output.parent.mkdir(parents=True, exist_ok=True) + output.write_text(json.dumps(result, indent=2, allow_nan=False) + "\n") + print(f"{'PASS' if result['passed'] else 'FAIL'}: {len(cases)} cases; {output}") + return 0 if result["passed"] else 1 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/experiments/plateau_diagnostic.cpp b/experiments/plateau_diagnostic.cpp index 935d124..3f5a42f 100644 --- a/experiments/plateau_diagnostic.cpp +++ b/experiments/plateau_diagnostic.cpp @@ -31,6 +31,7 @@ import mean_field; #include "plateau_positivity_strategy_diagnostics.hpp" #include "plateau_polynomial_certificate.hpp" #include "shape_field_weak_diagnostics.hpp" +#include "reduced_surface_diagnostics.hpp" using namespace mean_field; namespace { @@ -42,6 +43,8 @@ struct Options { bool strategyProbe = false, positivityProbe = false, positivityCertificate = false; bool preconditionerProbe = false, nonlinearStepProbe = false; bool weakShapeProbe = false, weakShapeFinal = false; + bool reducedSurfaceProbe = false, reducedSurfaceAlgebraCheck = false; + int surfaceDegree = 4; bool coupledTrialProbe = false, coupledSeedControl = false; bool dampedCoupledProbe = false; std::string coupledSearch = "near"; @@ -115,6 +118,12 @@ Options Parse(int argc, char **argv) { o.preconditionerProbe = true; else if (a == "--nonlinear-step-probe") o.nonlinearStepProbe = true; + else if (a == "--reduced-surface-probe") + o.reducedSurfaceProbe = true; + else if (a == "--reduced-surface-algebra-check") + o.reducedSurfaceAlgebraCheck = true; + else if (a == "--surface-degree") + o.surfaceDegree = std::stoi(v()); else if (a == "--weak-shape-probe") o.weakShapeProbe = true; else if (a == "--weak-shape-final") @@ -177,6 +186,18 @@ Options Parse(int argc, char **argv) { } if (o.output.empty() || o.n < 1 || o.n >= 5 || o.maxNewton < 0 || o.maxLinear < 1 || o.restart < 1) throw std::invalid_argument("Invalid options"); + if (o.reducedSurfaceProbe || o.reducedSurfaceAlgebraCheck) { + if ((o.reducedSurfaceProbe && o.reducedSurfaceAlgebraCheck) || o.maxNewton != 0 || + o.strategy != "baseline" || o.probes || o.pressureProbe || o.shapeCurvature || + o.surfaceTrace || o.strategyProbe || o.positivityProbe || o.positivityCertificate || + o.preconditionerProbe || o.nonlinearStepProbe || o.coupledTrialProbe || o.dampedCoupledProbe || + o.weakShapeProbe || o.weakShapeFinal || !o.probeDirection.empty() || + !o.replayCandidate.empty() || !o.captureIterations.empty() || !o.noCentral || + (o.surfaceDegree != 0 && o.surfaceDegree != 2 && o.surfaceDegree != 4 && + o.surfaceDegree != 6 && o.surfaceDegree != 8) || + !std::isfinite(o.linearTolerance) || o.linearTolerance <= 0.) + throw std::invalid_argument("Reduced-surface study requires a standalone unbordered frozen probe; degree 0 is the full control"); + } if (o.weakShapeProbe || o.weakShapeFinal) { if ((o.weakShapeProbe && o.weakShapeFinal) || o.probes || o.pressureProbe || o.shapeCurvature || o.surfaceTrace || o.strategyProbe || o.positivityProbe || o.positivityCertificate || @@ -234,8 +255,8 @@ Options Parse(int argc, char **argv) { throw std::invalid_argument("Newton preconditioner variants are independent, gravity-schur, stellar-ldu, or stellar-ldu-gravity-schur"); if (o.preconditionerVariant != "independent") { if (o.probes || o.pressureProbe || o.shapeCurvature || o.surfaceTrace || o.strategyProbe || - o.positivityProbe || o.positivityCertificate || (!o.probeState.empty() && !o.nonlinearStepProbe && !o.coupledTrialProbe && !o.dampedCoupledProbe && !o.weakShapeProbe) || - (!o.probeDirection.empty() && !o.coupledTrialProbe && !o.dampedCoupledProbe && !o.weakShapeProbe) || + o.positivityProbe || o.positivityCertificate || (!o.probeState.empty() && !o.nonlinearStepProbe && !o.coupledTrialProbe && !o.dampedCoupledProbe && !o.weakShapeProbe && !o.reducedSurfaceProbe) || + (!o.probeDirection.empty() && !o.coupledTrialProbe && !o.dampedCoupledProbe && !o.weakShapeProbe && !o.reducedSurfaceProbe) || !o.replayCandidate.empty()) throw std::invalid_argument("A custom Newton preconditioner cannot be combined with frozen or certificate probes"); if (o.strategy != "baseline" && o.strategy != "positivity" && o.strategy != "curvature" && @@ -529,7 +550,7 @@ int RunContext(Context &context, const Options &o, const Scales &scales) { << "\nrestart=" << o.restart << "\nprobe_restart=" << o.probeRestart << "\nfrozen_probe_backend=" << (o.preconditionerProbe ? "direct_mfem_fgmres_preconditioner_benchmark" : (o.coupledTrialProbe || o.dampedCoupledProbe) ? "supplied_outer_direction_verified_with_full_J_and_auxiliary_fgmres" : - (o.nonlinearStepProbe || ((o.weakShapeProbe || o.weakShapeFinal) && o.probeDirection.empty())) + (o.reducedSurfaceProbe || o.nonlinearStepProbe || ((o.weakShapeProbe || o.weakShapeFinal) && o.probeDirection.empty())) ? "prepared_fgmres_cold_with_independent_true_gate" : o.positivityCertificate ? "polynomial_certificate_no_linear_solve" : o.probeDirection.empty() ? "direct_mfem_fgmres_cold" : "supplied_direction_no_linear_solve") @@ -539,6 +560,7 @@ int RunContext(Context &context, const Options &o, const Scales &scales) { << "\ngravity_iterations=" << o.gravityIterations << "\ngravity_restart=" << o.gravityRestart << "\npreconditioner_probe=" << o.preconditionerProbe << "\nnonlinear_step_probe=" << o.nonlinearStepProbe + << "\nreduced_surface_probe=" << o.reducedSurfaceProbe << "\nsurface_degree=" << o.surfaceDegree << "\nweak_shape_probe=" << o.weakShapeProbe << "\nweak_shape_final=" << o.weakShapeFinal << "\ndamped_coupled_probe=" << o.dampedCoupledProbe << "\ncoupled_search=" << o.coupledSearch @@ -564,6 +586,23 @@ int RunContext(Context &context, const Options &o, const Scales &scales) { << "weak_shape_operator=full_monolithic\n" << "weak_shape_trial_policy=raw_prepare_only\n"; metadata.flush(); + if (o.reducedSurfaceProbe) { + metadata << "reduced_surface_input=" << (o.probeState.empty() ? "lane_emden_seed" : "physical_snapshot") << '\n' + << "reduced_surface_trial_policy=raw_prepare_only_no_commit\n" + << "probe_state=" << o.probeState << '\n'; + metadata.flush(); + solver::detail::StellarEquilibriumContextDiagnostics::WithState(context, [&](auto &s, const fem::FEM &f) { + mfem::Vector base(s.acceptedNormalizedState); + if (!o.probeState.empty()) { + const auto physical = LoadVector(o.probeState, s.AcceptedPhysicalState().Size()); + s.normalizedOperator->NormalizeState(physical, base); + } + experiment::plateau_diagnostics::ProbeReducedSurface(s, f, base, o.output, o.surfaceDegree, o.n, scales.K, { + .linear = {.relativeTolerance = o.linearTolerance, .absoluteTolerance = 0., + .maximumIterations = o.maxLinear}, .restartLength = o.restart}); + }); + return 0; + } if (o.weakShapeProbe) { metadata << "weak_shape_input=" << (o.probeState.empty() ? "lane_emden_seed" : "physical_snapshot") << '\n' << "probe_state=" << o.probeState << "\nprobe_direction=" << o.probeDirection << '\n'; @@ -1088,7 +1127,14 @@ int main(int argc, char **argv) { throw std::runtime_error("Single rank only"); mfem::Device device("cpu"); std::cout << std::setprecision(12); - result = Run(Parse(argc, argv)); + const auto options = Parse(argc, argv); + if (options.reducedSurfaceAlgebraCheck) { + const auto checks = experiment::plateau_diagnostics::RunReducedSurfaceAlgebraSelfCheck(MPI_COMM_WORLD); + std::cout << "Reduced surface algebra self-check passed: " << checks.checks + << " checks, maximum error " << checks.maximumError << std::endl; + } else { + result = Run(options); + } } catch (const std::exception &e) { std::cerr << "Diagnostic failure: " << e.what() << std::endl; result = 2; diff --git a/experiments/reduced_surface_basis.hpp b/experiments/reduced_surface_basis.hpp new file mode 100644 index 0000000..5a6de77 --- /dev/null +++ b/experiments/reduced_surface_basis.hpp @@ -0,0 +1,339 @@ +#pragma once + +// Include after import mean_field. This is an experiment-local coordinate +// construction, not a change to the stellar model or its deformation policy. +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include "plateau_coupled_coarse.hpp" + +namespace experiment::plateau_diagnostics { + +struct ReducedSurfaceBasis { + mfem::DenseMatrix physicalBasis; + mfem::DenseMatrix rawPhysicalBasis; + mfem::DenseMatrix rawGram; + mfem::DenseMatrix orthonormalGram; + mfem::DenseMatrix coefficientTransform; + int maximumDegree{}; + double physicalSurfaceArea{}; + double maximumGramOrthogonalityError{}; +}; + +namespace reduced_surface_basis_detail { +inline void Require(bool condition, const std::string &message) { + if (!condition) throw std::runtime_error("Reduced surface basis: " + message); +} + +inline bool Identical(const mfem::Vector &a, const mfem::Vector &b) { + return a.Size() == b.Size() && (a.Size() == 0 || + std::memcmp(a.GetData(), b.GetData(), sizeof(mfem::real_t) * a.Size()) == 0); +} + +inline std::ofstream Output(const std::filesystem::path &directory, std::string_view name, + std::string_view extension = ".txt") { + const auto path = directory / ("reduced_surface_" + std::string(name) + std::string(extension)); + Require(!std::filesystem::exists(path), "refusing to overwrite " + path.string()); + std::ofstream stream(path); + stream.exceptions(std::ios::badbit | std::ios::failbit); + stream << std::setprecision(std::numeric_limits::max_digits10); + return stream; +} + +inline void Save(const std::filesystem::path &directory, std::string_view name, + const mfem::DenseMatrix &matrix) { + auto stream = Output(directory, name); + stream << matrix.Height() << ' ' << matrix.Width() << '\n'; + for (int i = 0; i < matrix.Height(); ++i) + for (int j = 0; j < matrix.Width(); ++j) { + Require(std::isfinite(matrix(i, j)), "nonfinite matrix output"); + stream << matrix(i, j) << '\n'; + } +} + +// The normalized recurrence and signed-m convention match the existing +// analyze_shape_surface_modes.py analysis, including Condon--Shortley phase. +inline std::vector Harmonics(int degree, double x, double y, double z) { + std::array, 9> p{}; + const double phi = std::atan2(y, x), rxy = std::hypot(x, y); + p[0][0] = 1. / std::sqrt(4. * std::numbers::pi); + for (int m = 0; m <= degree; ++m) { + if (m > 0) p[m][m] = -std::sqrt((2. * m + 1.) / (2. * m)) * rxy * p[m - 1][m - 1]; + if (m < degree) p[m + 1][m] = std::sqrt(2. * m + 3.) * z * p[m][m]; + for (int l = m + 2; l <= degree; ++l) { + const double a = std::sqrt((4. * l * l - 1.) / (l * l - m * m)); + const double b = std::sqrt(((2. * l + 1.) * ((l - 1.) * (l - 1.) - m * m)) / + ((2. * l - 3.) * (l * l - m * m))); + p[l][m] = a * z * p[l - 1][m] - b * p[l - 2][m]; + } + } + std::vector values; + values.reserve((degree + 1) * (degree + 1)); + for (int l = 0; l <= degree; ++l) + for (int m = -l; m <= l; ++m) { + double value = p[l][std::abs(m)]; + if (m < 0) value *= std::sqrt(2.) * std::sin(-m * phi); + if (m > 0) value *= std::sqrt(2.) * std::cos(m * phi); + Require(std::isfinite(value), "nonfinite harmonic value"); + values.push_back(value); + } + return values; +} + +inline long double ColumnDot(const mfem::DenseMatrix &a, int i, + const mfem::DenseMatrix &b, int j) { + long double value = 0.; + for (int row = 0; row < a.Height(); ++row) + value += static_cast(a(row, i)) * b(row, j); + return value; +} + +struct Face { + int boundary{}; + mfem::Array scalarDofs; + mfem::Vector displacement; + mfem::DenseMatrix coefficients; +}; +} // namespace reduced_surface_basis_detail + +template +ReducedSurfaceBasis MakeReducedSurfaceBasis(State &state, const mean_field::fem::FEM &fem, + int maximumDegree, const std::filesystem::path &output) { + namespace d = reduced_surface_basis_detail; + using Schema = mean_field::utils::domain::CoreEnvelopeVacuumDomainSchema; + const auto start = std::chrono::steady_clock::now(); + d::Require(maximumDegree == 2 || maximumDegree == 4 || maximumDegree == 6 || maximumDegree == 8, + "maximum degree must be one of 2, 4, 6, 8"); + const auto &problem = state.Problem(); + int ranks = 0; + MPI_Comm_size(problem.GetCommunicator(), &ranks); + d::Require(ranks == 1, "this bounded experiment supports exactly one MPI rank"); + d::Require(state.normalizedOperator->IsPrepared() && problem.GetPhysicalOperator().IsPrepared(), + "a prepared frozen base is required"); + const auto generation = problem.GetPreparationGeneration(); + const mfem::Vector accepted(state.acceptedNormalizedState), acceptedResidual(state.acceptedNormalizedResidual), + acceptedPhysical(state.AcceptedPhysicalState()), trial(state.trialNormalizedState), + trialResidual(state.trialNormalizedResidual), candidatePhysical(state.candidatePhysicalState), + preparedPhysical(state.normalizedOperator->GetPhysicalState()), + preparedResidual(state.normalizedOperator->GetPhysicalResidual()), + stateFactors(state.normalizedOperator->GetNormalization().StateFactors()), + residualFactors(state.normalizedOperator->GetNormalization().ResidualFactors()); + const auto &physical = problem.GetPhysicalOperator(); + const mfem::Vector baseDisplacement(physical.GetGeneratedVolumeDisplacement()); + // This verifies the actual production radial trace, rather than relying + // only on its erased runtime descriptor. It does not prepare any state. + const auto radial = coupled_shape_detail::MakeDefaultRadialSurface(problem); + const int n = radial.parameterCount(), k = (maximumDegree + 1) * (maximumDegree + 1); + d::Require(n >= k && fem.mesh && fem.surfaceDeformationFes && fem.displacementFes, + "incomplete spaces or fewer surface parameters than requested modes"); + d::Require(!fem.surfaceDeformationFes->Nonconforming() && + !fem.surfaceDeformationFes->IsVariableOrder(), "conforming fixed-order scalar traces are required"); + const auto ¶meters = physical.GetSurfaceDeformationParameters(); + d::Require(parameters.Size() == n, "surface parameter ordering mismatch"); + ReducedSurfaceBasis result; + result.maximumDegree = maximumDegree; + result.rawPhysicalBasis.SetSize(n, k); + for (int i = 0; i < n; ++i) { + const auto values = d::Harmonics(maximumDegree, radial.radialDirection(i, 0), + radial.radialDirection(i, 1), radial.radialDirection(i, 2)); + for (int j = 0; j < k; ++j) result.rawPhysicalBasis(i, j) = values[j]; + d::Require(std::isfinite(parameters(i)) && radial.referenceRadius(i) + parameters(i) > 0., + "nonpositive or nonfinite nodal physical radius"); + } + + mfem::ParGridFunction displacement(fem.displacementFes.get()), scalar(fem.surfaceDeformationFes.get()); + displacement.SetFromTrueDofs(baseDisplacement); + std::vector faces; + constexpr int attribute = Schema::template boundary_attribute(); + constexpr int quadrature = 14; + int sampleCount = 0; + // Tagged internal interfaces are included by GetNBE; a Boundary-only face + // restriction would silently omit the stellar/vacuum interface here. + for (int b = 0; b < fem.surfaceDeformationFes->GetNBE(); ++b) { + if (fem.surfaceDeformationFes->GetBdrAttribute(b) != attribute) continue; + d::Face face; + face.boundary = b; + const auto &element = *fem.surfaceDeformationFes->GetBE(b); + mfem::Array vectorDofs; + d::Require(fem.surfaceDeformationFes->GetBdrElementVDofs(b, face.scalarDofs) == nullptr && + fem.displacementFes->GetBdrElementVDofs(b, vectorDofs) == nullptr, + "nontrivial boundary DOF transformations are unsupported"); + d::Require(face.scalarDofs.Size() == element.GetDof() && vectorDofs.Size() == 3 * element.GetDof(), + "the scalar and vector trace elements disagree"); + displacement.GetSubVector(vectorDofs, face.displacement); + face.coefficients.SetSize(element.GetDof(), k); + sampleCount += mfem::IntRules.Get(element.GetGeomType(), quadrature).GetNPoints(); + faces.push_back(std::move(face)); + } + d::Require(sampleCount > 0, "empty stellar surface quadrature"); + mfem::Vector mode(n), ambient(fem.surfaceDeformationFes->GetTrueVSize()), local; + for (int j = 0; j < k; ++j) { + for (int i = 0; i < n; ++i) mode(i) = result.rawPhysicalBasis(i, j); + radial.surfaceDofMap().scatter(mode, ambient); + scalar.SetFromTrueDofs(ambient); + for (auto &face : faces) { + scalar.GetSubVector(face.scalarDofs, local); + for (int i = 0; i < local.Size(); ++i) face.coefficients(i, j) = local(i); + } + } + mfem::DenseMatrix samples(sampleCount, k), weighted(sampleCount, k); + mfem::Vector weights(sampleCount), shape; + mfem::DenseMatrix derivative, currentJacobian(3, 2); + int sample = 0; + long double surfaceArea = 0.; + for (const auto &face : faces) { + const auto &element = *fem.surfaceDeformationFes->GetBE(face.boundary); + auto &transform = *fem.surfaceDeformationFes->GetBdrElementTransformation(face.boundary); + const auto &rule = mfem::IntRules.Get(element.GetGeomType(), quadrature); + const int ndof = element.GetDof(); + shape.SetSize(ndof); + derivative.SetSize(ndof, 2); + for (int q = 0; q < rule.GetNPoints(); ++q, ++sample) { + const auto &point = rule.IntPoint(q); + transform.SetIntPoint(&point); + element.CalcShape(point, shape); + element.CalcDShape(point, derivative); + currentJacobian = transform.Jacobian(); + // Different polynomial orders in the mesh and displacement are + // kept distinct: this is the exact trace Jacobian of X + u_h. + for (int c = 0; c < 3; ++c) + for (int a = 0; a < 2; ++a) + for (int i = 0; i < ndof; ++i) + currentJacobian(c, a) += face.displacement(i + c * ndof) * derivative(i, a); + const double weight = point.weight * currentJacobian.Weight(); + d::Require(std::isfinite(weight) && weight > 0., "invalid physical surface area weight"); + weights(sample) = weight; + surfaceArea += weight; + for (int j = 0; j < k; ++j) { + long double value = 0.; + for (int i = 0; i < ndof; ++i) value += static_cast(shape(i)) * face.coefficients(i, j); + samples(sample, j) = static_cast(value); + weighted(sample, j) = std::sqrt(weight) * samples(sample, j); + } + } + } + result.physicalSurfaceArea = static_cast(surfaceArea); + result.rawGram.SetSize(k, k); + for (int i = 0; i < k; ++i) + for (int j = 0; j < k; ++j) + result.rawGram(i, j) = static_cast(d::ColumnDot(weighted, i, weighted, j)); + + // Ordered, twice-reorthogonalized columns preserve nested harmonic spans. + // Rank loss rejects the requested space; it never silently changes K. + mfem::DenseMatrix orthogonal(weighted); + result.coefficientTransform.SetSize(k, k); + result.coefficientTransform = 0.; + std::vector pivots(k), relativePivots(k); + for (int j = 0; j < k; ++j) { + result.coefficientTransform(j, j) = 1.; + const double originalNorm = std::sqrt(static_cast(d::ColumnDot(weighted, j, weighted, j))); + for (int pass = 0; pass < 2; ++pass) + for (int i = 0; i < j; ++i) { + const double projection = static_cast(d::ColumnDot(orthogonal, i, orthogonal, j)); + for (int row = 0; row < sampleCount; ++row) + orthogonal(row, j) -= projection * orthogonal(row, i); + for (int row = 0; row < k; ++row) + result.coefficientTransform(row, j) -= projection * result.coefficientTransform(row, i); + } + pivots[j] = std::sqrt(static_cast(d::ColumnDot(orthogonal, j, orthogonal, j))); + relativePivots[j] = pivots[j] / originalNorm; + d::Require(std::isfinite(relativePivots[j]) && originalNorm > 0. && relativePivots[j] > 1e-12, + "the requested harmonic trace space is numerically rank deficient"); + for (int row = 0; row < sampleCount; ++row) orthogonal(row, j) /= pivots[j]; + for (int row = 0; row < k; ++row) result.coefficientTransform(row, j) /= pivots[j]; + } + result.physicalBasis.SetSize(n, k); + mfem::Mult(result.rawPhysicalBasis, result.coefficientTransform, result.physicalBasis); + mfem::DenseMatrix reconstructed(sampleCount, k); + mfem::Mult(weighted, result.coefficientTransform, reconstructed); + double reconstructionError = 0.; + for (int i = 0; i < sampleCount; ++i) + for (int j = 0; j < k; ++j) + reconstructionError = std::max(reconstructionError, std::abs(reconstructed(i, j) - orthogonal(i, j))); + result.orthonormalGram.SetSize(k, k); + for (int i = 0; i < k; ++i) + for (int j = 0; j < k; ++j) { + const double value = static_cast(d::ColumnDot(reconstructed, i, reconstructed, j)); + d::Require(std::isfinite(value), "nonfinite reconstructed mass Gram matrix"); + result.orthonormalGram(i, j) = value; + result.maximumGramOrthogonalityError = std::max(result.maximumGramOrthogonalityError, + std::abs(value - (i == j ? 1. : 0.))); + } + d::Require(result.maximumGramOrthogonalityError <= 1e-10 && reconstructionError <= 1e-10, + "physical mass orthonormalization failed its reconstruction check"); + mfem::Vector a(k), r(n), ba(n), btr(k); + for (int i = 0; i < k; ++i) a(i) = std::cos(.21 * (i + 1)); + for (int i = 0; i < n; ++i) r(i) = std::sin(.37 * (i + 1)); + result.physicalBasis.Mult(a, ba); + result.physicalBasis.MultTranspose(r, btr); + const double adjointError = std::abs((ba * r) - (a * btr)) / + std::max({1., ba.Norml2() * r.Norml2(), a.Norml2() * btr.Norml2()}); + d::Require(std::isfinite(adjointError) && adjointError <= 1e-12, "basis transpose identity failed"); + + const bool acceptedUnchanged = d::Identical(accepted, state.acceptedNormalizedState) && + d::Identical(acceptedResidual, state.acceptedNormalizedResidual) && + d::Identical(acceptedPhysical, state.AcceptedPhysicalState()); + const bool preparedUnchanged = generation == problem.GetPreparationGeneration() && + state.normalizedOperator->IsPrepared() && d::Identical(preparedPhysical, state.normalizedOperator->GetPhysicalState()) && + d::Identical(preparedResidual, state.normalizedOperator->GetPhysicalResidual()) && + d::Identical(baseDisplacement, physical.GetGeneratedVolumeDisplacement()) && + d::Identical(trial, state.trialNormalizedState) && d::Identical(trialResidual, state.trialNormalizedResidual) && + d::Identical(candidatePhysical, state.candidatePhysicalState); + const bool factorsUnchanged = d::Identical(stateFactors, state.normalizedOperator->GetNormalization().StateFactors()) && + d::Identical(residualFactors, state.normalizedOperator->GetNormalization().ResidualFactors()); + d::Require(acceptedUnchanged && preparedUnchanged && factorsUnchanged, "a solver state buffer changed"); + std::filesystem::create_directories(output); + d::Save(output, "basis_raw", result.rawPhysicalBasis); + d::Save(output, "basis_physical", result.physicalBasis); + d::Save(output, "gram_raw", result.rawGram); + d::Save(output, "gram_orthonormal", result.orthonormalGram); + d::Save(output, "coefficient_transform", result.coefficientTransform); + d::Save(output, "basis_samples_raw", samples); + auto weightOutput = d::Output(output, "basis_weights"); + weightOutput << sampleCount << '\n'; + for (int i = 0; i < sampleCount; ++i) weightOutput << weights(i) << '\n'; + auto nodes = d::Output(output, "basis_nodes", ".csv"); + nodes << "parameter,unit_radial_x,unit_radial_y,unit_radial_z,reference_radius,physical_radius,base_parameter\n"; + for (int i = 0; i < n; ++i) + nodes << i << ',' << radial.radialDirection(i, 0) << ',' << radial.radialDirection(i, 1) << ',' + << radial.radialDirection(i, 2) << ',' << radial.referenceRadius(i) << ',' + << radial.referenceRadius(i) + parameters(i) << ',' << parameters(i) << '\n'; + auto modes = d::Output(output, "basis_modes", ".csv"); + modes << "column,degree,m,pivot_norm,relative_pivot\n"; + int column = 0; + for (int l = 0; l <= maximumDegree; ++l) + for (int m = -l; m <= l; ++m, ++column) + modes << column << ',' << l << ',' << m << ',' << pivots[column] << ',' << relativePivots[column] << '\n'; + auto metadata = d::Output(output, "basis_metadata", ".csv"); + metadata << "metric,value\nmaximum_degree," << maximumDegree << "\nmodes," << k + << "\nsurface_parameters," << n << "\nquadrature," << quadrature + << "\nsurface_quadrature_points," << sampleCount << "\nphysical_surface_area," << result.physicalSurfaceArea + << "\nmaximum_gram_orthogonality_error," << result.maximumGramOrthogonalityError + << "\nmaximum_sample_reconstruction_error," << reconstructionError << "\nadjoint_relative_error," << adjointError + << "\naccepted_buffers_unchanged," << acceptedUnchanged << "\nprepared_state_unchanged," << preparedUnchanged + << "\nnormalization_factors_unchanged," << factorsUnchanged + << "\nmetric_scope,scalar_radial_amplitude_FE_trace_mass_on_current_physical_surface" + << "\nscalar_trace_order," << fem.surfaceDeformationFes->GetBE(faces.front().boundary)->GetOrder() + << "\nscalar_finite_element_collection," << fem.surfaceDeformationFes->FEColl()->Name() + << "\ndisplacement_finite_element_collection," << fem.displacementFes->FEColl()->Name() + << "\nstellar_surface_faces," << faces.size() + << "\nvector_displacement_metric_claimed,0\nbase_projection_performed,0\ntrial_preparations,0" + << "\nradial_trace_validation_actions,1\nrank_relative_pivot_threshold,1e-12" + << "\nconstruction_seconds," << std::chrono::duration(std::chrono::steady_clock::now() - start).count() << '\n'; + return result; +} +} // namespace experiment::plateau_diagnostics diff --git a/experiments/reduced_surface_diagnostics.hpp b/experiments/reduced_surface_diagnostics.hpp new file mode 100644 index 0000000..1625ba4 --- /dev/null +++ b/experiments/reduced_surface_diagnostics.hpp @@ -0,0 +1,270 @@ +#pragma once + +// A frozen, affine surface restriction tests approximation-space compatibility +// without changing the production problem or committing a filtered model. +#include "plateau_nonlinear_step_diagnostics.hpp" +#include "plateau_coupled_trial_physical.hpp" +#include "shape_field_physical_diagnostics.hpp" +#include "reduced_surface_basis.hpp" +#include "reduced_surface_operators.hpp" +#include + +namespace experiment::plateau_diagnostics { +namespace reduced_surface_probe_detail { +inline void SaveMatrix(const std::filesystem::path &path, const mfem::DenseMatrix &matrix) { + auto file = nonlinear_step_detail::Open(path); + file << matrix.Height() << ' ' << matrix.Width() << '\n'; + for (int i = 0; i < matrix.Height(); ++i) + for (int j = 0; j < matrix.Width(); ++j) { + if (!std::isfinite(matrix(i, j))) throw std::runtime_error("Nonfinite reduced surface matrix."); + file << matrix(i, j) << '\n'; + } +} +} // namespace reduced_surface_probe_detail + +template +void ProbeReducedSurface(State &state, const mean_field::fem::FEM &fem, + const mfem::Vector &normalizedBase, const std::filesystem::path &output, + int maximumDegree, double n, double K, NonlinearStepProbeOptions options = {}) { + using namespace mean_field; + namespace d = nonlinear_step_detail; + namespace pd = reduced_surface_probe_detail; + const auto start = d::Clock::now(); + options.linear.Validate(); + int ranks = 0; + MPI_Comm_size(state.Problem().GetCommunicator(), &ranks); + const int size = state.Problem().StateSize(); + if (ranks != 1 || size != state.Problem().EquationSize() || normalizedBase.Size() != size || + !d::Finite(normalizedBase) || options.restartLength < 1) + throw std::invalid_argument("Reduced surface probes require one rank and finite square root data."); + const auto algebra = RunReducedSurfaceAlgebraSelfCheck(state.Problem().GetCommunicator()); + const bool identity = maximumDegree == 0; + const mfem::Vector base(normalizedBase); + const auto &normalization = state.normalizedOperator->GetNormalization(); + const mfem::Vector stateFactors(normalization.StateFactors()), residualFactors(normalization.ResidualFactors()); + const mfem::Vector savedAccepted(state.acceptedNormalizedState), savedResidual(state.acceptedNormalizedResidual), + savedPhysical(state.AcceptedPhysicalState()), savedCorrection(state.normalizedCorrection), + savedPhysicalCorrection(state.physicalCorrection), savedVolumeDirection(state.volumeDisplacementDirection); + const double savedMinimumJacobian = state.acceptedMinimumJacobianDeterminant; + const auto values = state.Problem().GetManifest().valueBlocks(); + const auto rows = state.Problem().GetManifest().residualBlocks(); + const auto shapeValue = shape_field_physical_detail::RequireBlock(values, "surface_deformation", size); + const auto shapeRow = shape_field_physical_detail::RequireBlock(rows, "surface_shape_balance", size); + if (shapeValue.size != shapeRow.size) throw std::invalid_argument("Shape dimensions disagree."); + std::filesystem::create_directories(output); + const auto file = [&](const std::string &name) { return output / ("reduced_surface_" + name); }; + auto summary = d::Open(file("summary.csv")); + summary << "metric,value\n"; + const auto metric = [&](const char *key, auto value) { summary << key << ',' << value << '\n'; }; + auto layout = d::Open(file("layout.csv")); + layout << "kind,block,offset,size\n"; + for (const auto &b : values) layout << "value," << b.stableId << ',' << b.offset << ',' << b.size << '\n'; + for (const auto &b : rows) layout << "residual," << b.stableId << ',' << b.offset << ',' << b.size << '\n'; + layout.flush(); + mfem::Vector residual(size), direction(size), action(size), physicalBase, physicalDirection; + bool haveBase = false, finalized = false; + int preparations = 0; + double restoreError = 0.; + const auto prepare = [&](const mfem::Vector &candidate) { + state.trialNormalizedState = candidate; + ++preparations; + return state.PrepareTrial(); + }; + const auto restore = [&]() { + if (!prepare(base)) throw std::runtime_error("Reduced surface probe cannot restore base."); + if (!haveBase) { residual = state.trialNormalizedResidual; haveBase = true; } + if (!d::Finite(state.trialNormalizedResidual)) throw std::runtime_error("Nonfinite base residual."); + restoreError = std::max(restoreError, d::Difference(residual, state.trialNormalizedResidual)); + }; + const auto finish = [&]() { + finalized = true; + restore(); + state.physicalInverse->Refresh(); + state.normalizedCorrection = savedCorrection; + state.physicalCorrection = savedPhysicalCorrection; + state.volumeDisplacementDirection = savedVolumeDirection; + const bool acceptedSame = d::Identical(savedAccepted, state.acceptedNormalizedState) && + d::Identical(savedResidual, state.acceptedNormalizedResidual) && + d::Identical(savedPhysical, state.AcceptedPhysicalState()) && + savedMinimumJacobian == state.acceptedMinimumJacobianDeterminant; + const bool factorsSame = d::Identical(stateFactors, normalization.StateFactors()) && + d::Identical(residualFactors, normalization.ResidualFactors()); + metric("accepted_buffers_bitwise_identical", acceptedSame); + metric("normalization_factors_bitwise_identical", factorsSame); + metric("direction_buffers_bitwise_identical", d::Identical(savedCorrection, state.normalizedCorrection) && + d::Identical(savedPhysicalCorrection, state.physicalCorrection) && + d::Identical(savedVolumeDirection, state.volumeDisplacementDirection)); + metric("restored_residual_error", restoreError); + d::Save(file("restored_residual.txt"), state.trialNormalizedResidual); + metric("preparations", preparations); + metric("total_probe_seconds", d::Seconds(start)); + summary.flush(); + if (!acceptedSame || !factorsSame || restoreError != 0.) + throw std::runtime_error("Reduced surface probe changed frozen invariants."); + }; + try { + restore(); + state.physicalInverse->Refresh(); + physicalBase = state.normalizedOperator->GetPhysicalState(); + metric("schema_version", 1); + metric("maximum_degree", maximumDegree); + metric("identity_control", identity); + metric("full_size", size); + metric("shape_size", shapeValue.size); + metric("shape_value_offset", shapeValue.offset); + metric("shape_residual_offset", shapeRow.offset); + metric("linear_relative_tolerance", options.linear.relativeTolerance); + metric("linear_absolute_tolerance", options.linear.absoluteTolerance); + metric("linear_maximum_iterations", options.linear.maximumIterations); + metric("linear_restart", options.restartLength); + metric("accepted_commits", 0); + metric("algebra_self_check_passed", 1); + metric("algebra_self_check_count", algebra.checks); + metric("algebra_self_check_maximum_error", algebra.maximumError); + + std::unique_ptr coordinates; + if (!identity) { + auto basis = MakeReducedSurfaceBasis(state, fem, maximumDegree, output); + coordinates = std::make_unique(basis.physicalBasis, + stateFactors, residualFactors, shapeValue.offset, shapeRow.offset, state.Problem().GetCommunicator()); + pd::SaveMatrix(file("state_basis.txt"), coordinates->StateBasis()); + pd::SaveMatrix(file("residual_basis.txt"), coordinates->ResidualBasis()); + pd::SaveMatrix(file("state_transform.txt"), coordinates->StateTransform()); + pd::SaveMatrix(file("residual_transform.txt"), coordinates->ResidualTransform()); + } + const int reducedSize = identity ? size : coordinates->ReducedSize(); + metric("mode_count", identity ? shapeValue.size : coordinates->ModeCount()); + metric("reduced_size", reducedSize); + const auto restrictResidual = [&](const mfem::Vector &full, mfem::Vector &reduced) { + if (identity) reduced = full; else coordinates->RestrictResidual(full, reduced); + }; + const auto excluded = [&](const mfem::Vector &full, mfem::Vector &result) { + mfem::Vector reduced, retained; + restrictResidual(full, reduced); + if (identity) retained = reduced; else coordinates->InjectResidual(reduced, retained); + result = full; result -= retained; + }; + mfem::Vector reducedResidual, reducedDirection(reducedSize), reducedAction, reducedDefect; + restrictResidual(residual, reducedResidual); + reducedDirection = 0.; + mfem::Vector rhs(reducedResidual); rhs *= -1.; + const auto solve = [&](auto &operation, auto &inverse) { + auto backend = solver::linear::prepareLinearBackend( + solver::linear::FGMRES({.restartLength = options.restartLength, .printLevel = -1}), + operation, inverse, state.Problem().GetCommunicator()); + const auto result = backend.Solve(rhs, reducedDirection, options.linear); + metric("linear_status", int(result.status)); + metric("linear_iterations", result.iterations); + metric("linear_seconds", result.solveSeconds); + metric("linear_reported_true_norm", result.trueResidualNorm); + }; + if (identity) solve(*state.normalizedOperator, *state.normalizedInverse); + else { + ReducedSurfaceOperator operation(*state.normalizedOperator, *coordinates); + ReducedSurfacePreconditioner inverse(*state.normalizedInverse, *coordinates); + solve(operation, inverse); + } + restore(); + if (identity) direction = reducedDirection; else coordinates->InjectState(reducedDirection, direction); + if (!d::Finite(direction) || d::Norm(direction) == 0.) throw std::runtime_error("Invalid reduced direction."); + state.normalizedOperator->Mult(direction, action); + if (!d::Finite(action)) throw std::runtime_error("Nonfinite full Jacobian action."); + restrictResidual(action, reducedAction); + mfem::Vector defect(action); defect += residual; + reducedDefect = reducedAction; reducedDefect += reducedResidual; + mfem::Vector excludedResidual, excludedDefect; + excluded(residual, excludedResidual); excluded(defect, excludedDefect); + const double fullBaseNorm = d::Norm(residual), reducedBaseNorm = d::Norm(reducedResidual); + for (double norm : {fullBaseNorm, reducedBaseNorm, d::Norm(defect), d::Norm(reducedDefect), + d::Norm(direction), d::Norm(excludedResidual), d::Norm(excludedDefect)}) + if (!std::isfinite(norm)) throw std::runtime_error("Nonfinite reduced diagnostic norm."); + if (!(fullBaseNorm > 0.) || !(reducedBaseNorm > 0.)) + throw std::runtime_error("This comparison requires nonzero full and retained residuals."); + metric("full_base_norm", fullBaseNorm); + metric("reduced_base_norm", reducedBaseNorm); + metric("excluded_base_norm", d::Norm(excludedResidual)); + metric("full_relative_defect", d::Ratio(d::Norm(defect), fullBaseNorm)); + metric("reduced_relative_defect", d::Ratio(d::Norm(reducedDefect), reducedBaseNorm)); + metric("full_verified_gate", d::Norm(defect) <= options.linear.ConvergenceThreshold(fullBaseNorm)); + metric("reduced_verified_gate", d::Norm(reducedDefect) <= options.linear.ConvergenceThreshold(reducedBaseNorm)); + metric("full_merit_slope", d::Dot(residual, action, 0, size)); + metric("reduced_merit_slope", d::Dot(reducedResidual, reducedAction, 0, reducedSize)); + metric("direction_norm", d::Norm(direction)); + state.normalizedOperator->DenormalizeState(direction, physicalDirection); + for (const auto &[name, vector] : std::initializer_list>{ + {"base", &base}, {"physical_base", &physicalBase}, {"residual", &residual}, + {"direction", &direction}, {"physical_direction", &physicalDirection}, {"action", &action}, + {"defect", &defect}, {"state_factors", &stateFactors}, {"residual_factors", &residualFactors}, + {"reduced_residual", &reducedResidual}, {"reduced_direction", &reducedDirection}, + {"reduced_action", &reducedAction}, {"reduced_defect", &reducedDefect}, + {"excluded_residual", &excludedResidual}, {"excluded_defect", &excludedDefect}}) + d::Save(file(std::string(name) + ".txt"), *vector); + summary.flush(); + std::cout << "Reduced surface degree " << maximumDegree << ": linear iterations complete; full relative defect " + << d::Ratio(d::Norm(defect), fullBaseNorm) << ", reduced " + << d::Ratio(d::Norm(reducedDefect), reducedBaseNorm) << std::endl; + MeasureShapeFieldPhysical(state, fem, physicalBase, physicalDirection, output); + auto physical = d::Open(file("physical.csv")); + physical << "sample,quadrature,metric,value\n"; + const auto measurePhysical = [&](int sample) { + coupled_physical_detail::PreparedSnapshotView view(state); + experiment::polytrope_validation::PhysicalState snapshot(view, fem); + for (int q : {14, 18}) + for (const auto &[key, value] : MeasurePhysical(snapshot, n, K, q)) + physical << sample << ',' << q << ',' << key << ',' << value << '\n'; + physical.flush(); + }; + measurePhysical(-1); + auto trials = d::Open(file("trials.csv")); + trials << "sample,alpha,prepared,reason,stage,full_norm,reduced_norm,excluded_shape_norm," + "full_merit_pass,reduced_merit_pass,full_remainder_norm,reduced_remainder_norm\n"; + bool measuredTrial = false; + int sample = 0; + for (double alpha : {1., .25, .0625, 1./1024.}) { + // The affine base is retained verbatim: no projection of q0 and no + // enthalpy repair, SOC, or accepted-state update belongs in this test. + mfem::Vector candidate(base); candidate.Add(alpha, direction); + const auto prepared = prepare(candidate); + d::Save(file("candidate_a" + std::to_string(sample) + ".txt"), candidate); + if (prepared) { + const mfem::Vector trial(state.trialNormalizedResidual); + if (!d::Finite(trial)) throw std::runtime_error("Nonfinite reduced trial residual."); + mfem::Vector reducedTrial, omitted, remainder(trial), reducedRemainder; + restrictResidual(trial, reducedTrial); excluded(trial, omitted); + remainder -= residual; remainder.Add(-alpha, action); + restrictResidual(remainder, reducedRemainder); + const double fullNorm = d::Norm(trial), reducedNorm = d::Norm(reducedTrial); + trials << sample << ',' << alpha << ",1,-1,-1," << fullNorm << ',' << reducedNorm << ',' + << d::Norm(omitted) << ',' << (fullNorm <= (1.-1e-4*alpha)*fullBaseNorm) << ',' + << (reducedNorm <= (1.-1e-4*alpha)*reducedBaseNorm) << ',' + << d::Norm(remainder) << ',' << d::Norm(reducedRemainder) << '\n'; + d::Save(file("residual_a" + std::to_string(sample) + ".txt"), trial); + d::Save(file("reduced_residual_a" + std::to_string(sample) + ".txt"), reducedTrial); + mfem::Vector physicalCandidate; + state.normalizedOperator->DenormalizeState(candidate, physicalCandidate); + d::Save(file("physical_candidate_a" + std::to_string(sample) + ".txt"), physicalCandidate); + if (!measuredTrial) { measurePhysical(sample); measuredTrial = true; } + } else { + trials << sample << ',' << alpha << ",0," << int(prepared.error().reason) << ',' + << int(prepared.error().stage) << ",nan,nan,nan,0,0,nan,nan\n"; + } + trials.flush(); + ++sample; + } + restore(); + mfem::Vector restoredAction(size); + state.normalizedOperator->Mult(direction, restoredAction); + d::Save(file("restored_action.txt"), restoredAction); + metric("restored_action_error", d::Difference(action, restoredAction)); + if (!d::Identical(action, restoredAction)) throw std::runtime_error("Restored Jacobian changed."); + finish(); + metric("probe_exception", 0); + summary.flush(); + } catch (...) { + const auto failure = std::current_exception(); + metric("probe_exception", 1); + if (!finalized) finish(); + std::rethrow_exception(failure); + } +} +} // namespace experiment::plateau_diagnostics diff --git a/experiments/reduced_surface_operators.hpp b/experiments/reduced_surface_operators.hpp new file mode 100644 index 0000000..98f71a9 --- /dev/null +++ b/experiments/reduced_surface_operators.hpp @@ -0,0 +1,527 @@ +#pragma once + +// Include after importing mean_field. These maps change the discrete shape +// space; they are not an invertible preconditioner of the original equations. +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +namespace experiment::plateau_diagnostics { + struct ReducedSurfaceCoordinateOptions final { + double relativeRankTolerance{1.e-12}; + double orthogonalityTolerance{5.e-12}; + }; + + struct ReducedSurfaceCoordinateReport final { + double stateOrthogonalityError{0.0}; + double residualOrthogonalityError{0.0}; + double stateReconstructionError{0.0}; + double residualReconstructionError{0.0}; + double stateMinimumIndependentFraction{1.0}; + double residualMinimumIndependentFraction{1.0}; + }; + + namespace reduced_surface_detail { + inline void RequireFinite(const mfem::Vector &value, const char *description) { + for (int i = 0; i < value.Size(); ++i) + if (!std::isfinite(value[i])) throw std::invalid_argument(description); + } + + inline double StableNorm(const mfem::Vector &value) { + double norm = 0.0; + for (int i = 0; i < value.Size(); ++i) norm = std::hypot(norm, value[i]); + return norm; + } + + inline double DotColumn(const mfem::DenseMatrix &matrix, int column, const mfem::Vector &value) { + long double result = 0.0L; + for (int i = 0; i < matrix.Height(); ++i) + result += static_cast(matrix(i, column)) * value[i]; + return static_cast(result); + } + + struct QRReport final { + double orthogonalityError{0.0}; + double reconstructionError{0.0}; + double minimumIndependentFraction{1.0}; + }; + + // Reorthogonalization preserves the advertised coordinate identities + // without a new dense linear-algebra dependency. Rejecting a dependent + // mode keeps the scientific degree/mode ordering explicit in artifacts. + inline QRReport ThinQR( + const mfem::DenseMatrix &input, + mfem::DenseMatrix &basis, + mfem::DenseMatrix &transform, + ReducedSurfaceCoordinateOptions options + ) { + const int rows = input.Height(), columns = input.Width(); + basis.SetSize(rows, columns); + transform.SetSize(columns, columns); + basis = 0.0; + transform = 0.0; + QRReport report; + mfem::Vector column(rows); + for (int j = 0; j < columns; ++j) { + for (int i = 0; i < rows; ++i) column[i] = input(i, j); + RequireFinite(column, "The reduced surface basis contains nonfinite entries."); + const double originalNorm = StableNorm(column); + if (!std::isfinite(originalNorm) || originalNorm <= 0.0) + throw std::invalid_argument("The reduced surface basis has a zero or overflowing column."); + for (int pass = 0; pass < 2; ++pass) { + for (int previous = 0; previous < j; ++previous) { + const double projection = DotColumn(basis, previous, column); + if (!std::isfinite(projection)) + throw std::invalid_argument("The reduced surface QR projection is nonfinite."); + transform(previous, j) += projection; + for (int i = 0; i < rows; ++i) column[i] -= projection * basis(i, previous); + } + } + const double independentNorm = StableNorm(column); + const double fraction = independentNorm / originalNorm; + if (!std::isfinite(fraction) || fraction <= options.relativeRankTolerance) + throw std::invalid_argument("The reduced surface basis is numerically rank deficient."); + report.minimumIndependentFraction = std::min(report.minimumIndependentFraction, fraction); + transform(j, j) = independentNorm; + for (int i = 0; i < rows; ++i) basis(i, j) = column[i] / independentNorm; + } + double inputNorm = 0.0, errorNorm = 0.0; + for (int j = 0; j < columns; ++j) { + for (int l = 0; l < columns; ++l) { + long double product = 0.0L; + for (int i = 0; i < rows; ++i) + product += static_cast(basis(i, j)) * basis(i, l); + report.orthogonalityError = std::max( + report.orthogonalityError, std::abs(static_cast(product) - (j == l ? 1.0 : 0.0)) + ); + } + for (int i = 0; i < rows; ++i) { + long double reconstructed = 0.0L; + for (int l = 0; l <= j; ++l) + reconstructed += static_cast(basis(i, l)) * transform(l, j); + inputNorm = std::hypot(inputNorm, input(i, j)); + errorNorm = std::hypot(errorNorm, static_cast(reconstructed) - input(i, j)); + } + } + report.reconstructionError = errorNorm / inputNorm; + if (!std::isfinite(report.reconstructionError) || + !std::isfinite(report.orthogonalityError) || + report.orthogonalityError > options.orthogonalityTolerance || + report.reconstructionError > options.orthogonalityTolerance) + throw std::invalid_argument("The reduced surface QR failed its coordinate consistency check."); + return report; + } + } // namespace reduced_surface_detail + + class ReducedSurfaceCoordinates final { + public: + ReducedSurfaceCoordinates( + const mfem::DenseMatrix &physicalBasis, + const mfem::Vector &stateFactors, + const mfem::Vector &residualFactors, + int shapeValueOffset, + int shapeResidualOffset, + MPI_Comm communicator, + ReducedSurfaceCoordinateOptions options = {} + ) + : m_physicalBasis(physicalBasis), m_stateFactors(stateFactors), m_residualFactors(residualFactors), + m_shapeValueOffset(shapeValueOffset), m_shapeResidualOffset(shapeResidualOffset) { + int ranks = 0; + if (MPI_Comm_size(communicator, &ranks) != MPI_SUCCESS || ranks != 1) + throw std::invalid_argument("The reduced surface experiment currently requires one MPI rank."); + if (!std::isfinite(options.relativeRankTolerance) || options.relativeRankTolerance <= 0.0 || + options.relativeRankTolerance >= 1.0 || !std::isfinite(options.orthogonalityTolerance) || + options.orthogonalityTolerance <= 0.0 || options.orthogonalityTolerance >= 1.0) + throw std::invalid_argument("Invalid reduced surface coordinate tolerances."); + const int size = stateFactors.Size(), shape = physicalBasis.Height(), modes = physicalBasis.Width(); + if (size != residualFactors.Size() || shape <= 0 || modes <= 0 || modes > shape || shape > size || + shapeValueOffset < 0 || shapeValueOffset > size - shape || + shapeResidualOffset < 0 || shapeResidualOffset > size - shape) + throw std::invalid_argument("Invalid reduced surface dimensions or distinct block offsets."); + for (int i = 0; i < size; ++i) + if (!std::isfinite(stateFactors[i]) || stateFactors[i] <= 0.0 || + !std::isfinite(residualFactors[i]) || residualFactors[i] <= 0.0) + throw std::invalid_argument("Reduced surface normalization factors must be finite and positive."); + + mfem::DenseMatrix stateImage(shape, modes), residualDualImage(shape, modes); + for (int j = 0; j < modes; ++j) { + for (int i = 0; i < shape; ++i) { + stateImage(i, j) = stateFactors[shapeValueOffset + i] * physicalBasis(i, j); + residualDualImage(i, j) = physicalBasis(i, j) / residualFactors[shapeResidualOffset + i]; + } + } + // State and residual scales act on different layouts and are not + // interchangeable. Q_r^T Fhat = T_r^{-T} B^T Fphysical preserves + // precisely the physical Galerkin equations with an invertible + // reduced residual scaling, even when Q_x differs from Q_r. + const auto state = reduced_surface_detail::ThinQR(stateImage, m_stateBasis, m_stateTransform, options); + const auto residual = reduced_surface_detail::ThinQR( + residualDualImage, m_residualBasis, m_residualTransform, options + ); + m_report = { + state.orthogonalityError, residual.orthogonalityError, + state.reconstructionError, residual.reconstructionError, + state.minimumIndependentFraction, residual.minimumIndependentFraction + }; + } + + [[nodiscard]] int FullSize() const noexcept { return m_stateFactors.Size(); } + [[nodiscard]] int ReducedSize() const noexcept { return FullSize() - ShapeSize() + ModeCount(); } + [[nodiscard]] int ShapeSize() const noexcept { return m_physicalBasis.Height(); } + [[nodiscard]] int ModeCount() const noexcept { return m_physicalBasis.Width(); } + [[nodiscard]] int ShapeValueOffset() const noexcept { return m_shapeValueOffset; } + [[nodiscard]] int ShapeResidualOffset() const noexcept { return m_shapeResidualOffset; } + [[nodiscard]] const mfem::DenseMatrix &PhysicalBasis() const noexcept { return m_physicalBasis; } + [[nodiscard]] const mfem::DenseMatrix &StateBasis() const noexcept { return m_stateBasis; } + [[nodiscard]] const mfem::DenseMatrix &ResidualBasis() const noexcept { return m_residualBasis; } + [[nodiscard]] const mfem::DenseMatrix &StateTransform() const noexcept { return m_stateTransform; } + [[nodiscard]] const mfem::DenseMatrix &ResidualTransform() const noexcept { return m_residualTransform; } + [[nodiscard]] const ReducedSurfaceCoordinateReport &Report() const noexcept { return m_report; } + + [[nodiscard]] bool MatchesNormalization( + const mfem::Vector &stateFactors, const mfem::Vector &residualFactors + ) const noexcept { + if (stateFactors.Size() != FullSize() || residualFactors.Size() != FullSize()) return false; + for (int i = 0; i < FullSize(); ++i) + if (stateFactors[i] != m_stateFactors[i] || residualFactors[i] != m_residualFactors[i]) return false; + return true; + } + + void InjectState(const mfem::Vector &reduced, mfem::Vector &full) const { + Inject(reduced, full, m_shapeValueOffset, m_stateBasis); + } + void RestrictState(const mfem::Vector &full, mfem::Vector &reduced) const { + Restrict(full, reduced, m_shapeValueOffset, m_stateBasis); + } + void RestrictResidual(const mfem::Vector &full, mfem::Vector &reduced) const { + Restrict(full, reduced, m_shapeResidualOffset, m_residualBasis); + } + void InjectResidual(const mfem::Vector &reduced, mfem::Vector &full) const { + Inject(reduced, full, m_shapeResidualOffset, m_residualBasis); + } + + private: + void Inject( + const mfem::Vector &reduced, mfem::Vector &full, int offset, const mfem::DenseMatrix &basis + ) const { + if (reduced.Size() != ReducedSize() || (&reduced == &full)) + throw std::invalid_argument("Reduced surface injection requires distinct correctly sized vectors."); + reduced_surface_detail::RequireFinite(reduced, "Cannot inject a nonfinite reduced surface vector."); + full.SetSize(FullSize()); + for (int i = 0; i < offset; ++i) full[i] = reduced[i]; + for (int i = 0; i < ShapeSize(); ++i) { + long double value = 0.0L; + for (int j = 0; j < ModeCount(); ++j) + value += static_cast(basis(i, j)) * reduced[offset + j]; + full[offset + i] = static_cast(value); + } + for (int i = offset + ShapeSize(); i < FullSize(); ++i) + full[i] = reduced[i - ShapeSize() + ModeCount()]; + reduced_surface_detail::RequireFinite(full, "Reduced surface injection overflowed."); + } + + void Restrict( + const mfem::Vector &full, mfem::Vector &reduced, int offset, const mfem::DenseMatrix &basis + ) const { + if (full.Size() != FullSize() || (&full == &reduced)) + throw std::invalid_argument("Reduced surface restriction requires distinct correctly sized vectors."); + reduced_surface_detail::RequireFinite(full, "Cannot restrict a nonfinite full surface vector."); + reduced.SetSize(ReducedSize()); + for (int i = 0; i < offset; ++i) reduced[i] = full[i]; + for (int j = 0; j < ModeCount(); ++j) { + long double value = 0.0L; + for (int i = 0; i < ShapeSize(); ++i) + value += static_cast(basis(i, j)) * full[offset + i]; + reduced[offset + j] = static_cast(value); + } + for (int i = offset + ShapeSize(); i < FullSize(); ++i) + reduced[i - ShapeSize() + ModeCount()] = full[i]; + reduced_surface_detail::RequireFinite(reduced, "Reduced surface restriction overflowed."); + } + + mfem::DenseMatrix m_physicalBasis, m_stateBasis, m_residualBasis, m_stateTransform, m_residualTransform; + mfem::Vector m_stateFactors, m_residualFactors; + int m_shapeValueOffset, m_shapeResidualOffset; + ReducedSurfaceCoordinateReport m_report; + }; + + template FullOperator> + class ReducedSurfaceOperator final : public mfem::Operator { + public: + ReducedSurfaceOperator(const FullOperator &operation, const ReducedSurfaceCoordinates &coordinates) + : mfem::Operator(coordinates.ReducedSize()), m_operation(&operation), m_coordinates(&coordinates), + m_fullDirection(coordinates.FullSize()), m_fullAction(coordinates.FullSize()) { + Refresh(); + } + + // Refresh does not prepare the physical problem. It adopts an already + // prepared generation while requiring the frozen coordinate scales. + void Refresh() { + m_current = false; + if (m_operation->Height() != m_coordinates->FullSize() || + m_operation->Width() != m_coordinates->FullSize()) + throw std::invalid_argument("Reduced surface operator needs the matching square full Jacobian."); + if constexpr (requires { m_operation->IsPrepared(); }) + if (!m_operation->IsPrepared()) throw std::invalid_argument("The full Jacobian is not prepared."); + if constexpr (requires { m_operation->GetNormalization(); }) { + const auto &normalization = m_operation->GetNormalization(); + if (!m_coordinates->MatchesNormalization(normalization.StateFactors(), normalization.ResidualFactors())) + throw std::invalid_argument("The frozen reduced surface normalization has changed."); + } + if constexpr (requires { m_operation->GetProblem().GetPreparationGeneration(); }) + m_generation = m_operation->GetProblem().GetPreparationGeneration(); + m_current = true; + } + + [[nodiscard]] bool IsCurrent() const { + if (!m_current || m_operation->Height() != m_coordinates->FullSize() || + m_operation->Width() != m_coordinates->FullSize()) return false; + if constexpr (requires { m_operation->IsPrepared(); }) + if (!m_operation->IsPrepared()) return false; + if constexpr (requires { m_operation->GetProblem().GetPreparationGeneration(); }) + if (m_generation != m_operation->GetProblem().GetPreparationGeneration()) return false; + if constexpr (requires { m_operation->GetNormalization(); }) { + const auto &normalization = m_operation->GetNormalization(); + if (!m_coordinates->MatchesNormalization(normalization.StateFactors(), normalization.ResidualFactors())) + return false; + } + return true; + } + + void Mult(const mfem::Vector &direction, mfem::Vector &action) const override { + if (!IsCurrent()) throw std::invalid_argument("The reduced surface Jacobian is stale."); + m_coordinates->InjectState(direction, m_fullDirection); + m_operation->Mult(m_fullDirection, m_fullAction); + m_coordinates->RestrictResidual(m_fullAction, action); + } + + [[nodiscard]] const FullOperator &GetFullOperator() const noexcept { return *m_operation; } + + private: + const FullOperator *m_operation; + const ReducedSurfaceCoordinates *m_coordinates; + mutable mfem::Vector m_fullDirection, m_fullAction; + std::uint64_t m_generation{0}; + bool m_current{false}; + }; + + template + requires std::derived_from && + mean_field::solver::LinearPreconditionerApplicationContractAvailable + class ReducedSurfacePreconditioner final : public mfem::Solver { + public: + // Projecting an inverse is an approximation, not an inverse identity. + // In particular a flexible underlying application remains flexible. + static constexpr auto applicationContract = + mean_field::solver::linearPreconditionerApplicationContract; + + ReducedSurfacePreconditioner(FullInverse &inverse, const ReducedSurfaceCoordinates &coordinates) + : mfem::Solver(coordinates.ReducedSize()), m_inverse(&inverse), m_coordinates(&coordinates), + m_fullResidual(coordinates.FullSize()), m_fullCorrection(coordinates.FullSize()) { + if (!IsCurrent()) throw std::invalid_argument("The full reduced-surface inverse is not current."); + iterative_mode = false; + } + + void SetOperator(const mfem::Operator &operation) override { + if (operation.Height() != Height() || operation.Width() != Width()) + throw std::invalid_argument("The projected inverse received an incompatible reduced operator."); + // The physical inverse is already bound to the authoritative full + // operator. Passing this reduced operator into it would be invalid. + if (!IsCurrent()) throw std::invalid_argument("The projected inverse is stale."); + } + + [[nodiscard]] bool IsCurrent() const { + if (m_inverse->Height() != m_coordinates->FullSize() || + m_inverse->Width() != m_coordinates->FullSize() || m_inverse->iterative_mode) return false; + if constexpr (requires { m_inverse->IsCurrent(); }) + if (!m_inverse->IsCurrent()) return false; + return true; + } + + void Mult(const mfem::Vector &residual, mfem::Vector &correction) const override { + if (!IsCurrent()) throw std::invalid_argument("The projected inverse is stale."); + m_coordinates->InjectResidual(residual, m_fullResidual); + m_fullCorrection = 0.0; + m_inverse->Mult(m_fullResidual, m_fullCorrection); + m_coordinates->RestrictState(m_fullCorrection, correction); + } + + private: + FullInverse *m_inverse; + const ReducedSurfaceCoordinates *m_coordinates; + mutable mfem::Vector m_fullResidual, m_fullCorrection; + }; + + struct ReducedSurfaceAlgebraSelfCheckReport final { + int checks{0}; + double maximumError{0.0}; + bool rejectedDependentBasis{false}; + bool rejectedNonfiniteBasis{false}; + bool rejectedStaleOperator{false}; + }; + + namespace reduced_surface_detail { + class SelfCheckOperator final : public mfem::Operator { + public: + explicit SelfCheckOperator(int size) : mfem::Operator(size), matrix(size) { + matrix = 0.0; + for (int i = 0; i < size; ++i) { + matrix(i, i) = 2.0 + i; + for (int j = 0; j < size; ++j) + if (i != j) matrix(i, j) = 0.03 * (1 + i - 2 * j); + } + } + void Mult(const mfem::Vector &x, mfem::Vector &y) const override { + ++applications; + y.SetSize(Height()); + matrix.Mult(x, y); + } + [[nodiscard]] bool IsPrepared() const noexcept { return true; } + [[nodiscard]] const SelfCheckOperator &GetProblem() const noexcept { return *this; } + [[nodiscard]] std::uint64_t GetPreparationGeneration() const noexcept { return generation; } + mfem::DenseMatrix matrix; + std::uint64_t generation{1}; + mutable int applications{0}; + }; + + class SelfCheckInverse final : public mfem::Solver { + public: + static constexpr auto applicationContract = + mean_field::preconditioning::ApplicationContract::stationary_linear; + explicit SelfCheckInverse(int size) : mfem::Solver(size) { iterative_mode = false; } + void SetOperator(const mfem::Operator &) override {} + [[nodiscard]] bool IsCurrent() const noexcept { return true; } + void Mult(const mfem::Vector &x, mfem::Vector &y) const override { + ++applications; + y.SetSize(Height()); + for (int i = 0; i < Height(); ++i) y[i] = x[i] / (2.0 + i); + } + mutable int applications{0}; + }; + } // namespace reduced_surface_detail + + // Tiny deterministic algebra only: this exercises different value/row + // offsets and non-proportional scales without constructing any FEM state. + [[nodiscard]] inline ReducedSurfaceAlgebraSelfCheckReport RunReducedSurfaceAlgebraSelfCheck(MPI_Comm comm) { + ReducedSurfaceAlgebraSelfCheckReport report; + const auto check = [&report](double error) { + ++report.checks; + if (!std::isfinite(error) || error > 2.e-12) + throw std::runtime_error("Reduced surface algebra self-check failed."); + report.maximumError = std::max(report.maximumError, std::abs(error)); + }; + const auto vectorError = [&check](const mfem::Vector &a, const mfem::Vector &b) { + if (a.Size() != b.Size()) throw std::runtime_error("Reduced surface self-check size mismatch."); + double error = 0.0; + for (int i = 0; i < a.Size(); ++i) error = std::max(error, std::abs(a[i] - b[i])); + check(error); + }; + mfem::Vector stateFactors(7), residualFactors(7), full(7), reduced(6); + for (int i = 0; i < 7; ++i) { + stateFactors[i] = 0.7 + 0.3 * i; + residualFactors[i] = 1.9 - 0.13 * i; + full[i] = 0.2 - 0.17 * i; + } + for (int i = 0; i < 6; ++i) reduced[i] = 0.11 * (i + 1); + mfem::DenseMatrix physical(3, 2); + physical(0, 0) = 1.0; physical(1, 0) = 0.5; physical(2, 0) = 0.0; + physical(0, 1) = 0.2; physical(1, 1) = -0.3; physical(2, 1) = 1.0; + ReducedSurfaceCoordinates coordinates(physical, stateFactors, residualFactors, 1, 4, comm); + check(coordinates.Report().stateOrthogonalityError); + check(coordinates.Report().residualOrthogonalityError); + mfem::Vector injected, restricted, roundtrip, other; + coordinates.InjectState(reduced, injected); + coordinates.RestrictState(injected, roundtrip); + vectorError(reduced, roundtrip); + coordinates.RestrictState(full, restricted); + check(std::abs(injected * full - reduced * restricted)); + coordinates.InjectResidual(reduced, injected); + coordinates.RestrictResidual(injected, roundtrip); + vectorError(reduced, roundtrip); + coordinates.RestrictResidual(full, restricted); + check(std::abs(injected * full - reduced * restricted)); + + // Verify the physical dual identity independently of InjectResidual: + // T_r^T rhat_shape = B^T (D_r^{-1} Fhat)_shape. + for (int j = 0; j < 2; ++j) { + double lhs = 0.0, rhs = 0.0; + for (int l = 0; l < 2; ++l) lhs += coordinates.ResidualTransform()(l, j) * restricted[4 + l]; + for (int i = 0; i < 3; ++i) rhs += physical(i, j) * full[4 + i] / residualFactors[4 + i]; + check(std::abs(lhs - rhs)); + } + // The two QR bases really differ in this fixture; transpose reuse + // would therefore fail the physical dual identity above. + double mapDifference = 0.0; + for (int j = 0; j < 2; ++j) + for (int i = 0; i < 3; ++i) + mapDifference += std::abs(coordinates.StateBasis()(i, j) - coordinates.ResidualBasis()(i, j)); + if (mapDifference < 1.e-3) throw std::runtime_error("Reduced surface dual-map fixture is ineffective."); + ++report.checks; + + reduced_surface_detail::SelfCheckOperator fullOperator(7); + ReducedSurfaceOperator operation(fullOperator, coordinates); + operation.Mult(reduced, roundtrip); + coordinates.InjectState(reduced, injected); + fullOperator.Mult(injected, other); + coordinates.RestrictResidual(other, restricted); + vectorError(roundtrip, restricted); + reduced_surface_detail::SelfCheckInverse fullInverse(7); + ReducedSurfacePreconditioner inverse(fullInverse, coordinates); + inverse.SetOperator(operation); + const int oldApplications = fullOperator.applications; + inverse.Mult(reduced, roundtrip); + coordinates.InjectResidual(reduced, injected); + fullInverse.Mult(injected, other); + coordinates.RestrictState(other, restricted); + vectorError(roundtrip, restricted); + check(fullOperator.applications - oldApplications); + + ++fullOperator.generation; + try { operation.Mult(reduced, roundtrip); } + catch (const std::invalid_argument &) { report.rejectedStaleOperator = true; } + if (!report.rejectedStaleOperator || operation.IsCurrent()) + throw std::runtime_error("Reduced surface stale-generation guard failed."); + ++report.checks; + operation.Refresh(); + operation.Mult(reduced, roundtrip); + + mfem::DenseMatrix identity(3); + identity = 0.0; + for (int i = 0; i < 3; ++i) identity(i, i) = 1.0; + ReducedSurfaceCoordinates fullCoordinates(identity, stateFactors, residualFactors, 1, 4, comm); + fullCoordinates.InjectState(full, injected); + vectorError(full, injected); + fullCoordinates.RestrictResidual(full, restricted); + vectorError(full, restricted); + ReducedSurfaceOperator fullRecovery(fullOperator, fullCoordinates); + fullRecovery.Mult(full, roundtrip); + fullOperator.Mult(full, other); + vectorError(roundtrip, other); + ReducedSurfacePreconditioner inverseRecovery(fullInverse, fullCoordinates); + inverseRecovery.Mult(full, roundtrip); + fullInverse.Mult(full, other); + vectorError(roundtrip, other); + + mfem::DenseMatrix invalid(physical); + for (int i = 0; i < 3; ++i) invalid(i, 1) = invalid(i, 0); + try { ReducedSurfaceCoordinates rejected(invalid, stateFactors, residualFactors, 1, 4, comm); } + catch (const std::invalid_argument &) { report.rejectedDependentBasis = true; } + invalid = physical; + invalid(0, 0) = std::numeric_limits::quiet_NaN(); + try { ReducedSurfaceCoordinates rejected(invalid, stateFactors, residualFactors, 1, 4, comm); } + catch (const std::invalid_argument &) { report.rejectedNonfiniteBasis = true; } + if (!report.rejectedDependentBasis || !report.rejectedNonfiniteBasis) + throw std::runtime_error("Reduced surface basis validity guards failed."); + report.checks += 2; + return report; + } +} // namespace experiment::plateau_diagnostics