Add frozen smooth-surface restriction study with physical dual maps
This commit is contained in:
64
experiments/REDUCED_SURFACE_STUDY_2026-09-18.md
Normal file
64
experiments/REDUCED_SURFACE_STUDY_2026-09-18.md
Normal file
@@ -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.
|
||||
517
experiments/audit_reduced_surface.py
Normal file
517
experiments/audit_reduced_surface.py
Normal file
@@ -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())
|
||||
@@ -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;
|
||||
|
||||
339
experiments/reduced_surface_basis.hpp
Normal file
339
experiments/reduced_surface_basis.hpp
Normal file
@@ -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 <algorithm>
|
||||
#include <array>
|
||||
#include <chrono>
|
||||
#include <cmath>
|
||||
#include <cstring>
|
||||
#include <filesystem>
|
||||
#include <fstream>
|
||||
#include <iomanip>
|
||||
#include <limits>
|
||||
#include <numbers>
|
||||
#include <stdexcept>
|
||||
#include <string>
|
||||
#include <string_view>
|
||||
#include <vector>
|
||||
|
||||
#include <mfem.hpp>
|
||||
#include <mpi.h>
|
||||
#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<double>::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<double> Harmonics(int degree, double x, double y, double z) {
|
||||
std::array<std::array<double, 9>, 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<double> 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<long double>(a(row, i)) * b(row, j);
|
||||
return value;
|
||||
}
|
||||
|
||||
struct Face {
|
||||
int boundary{};
|
||||
mfem::Array<int> scalarDofs;
|
||||
mfem::Vector displacement;
|
||||
mfem::DenseMatrix coefficients;
|
||||
};
|
||||
} // namespace reduced_surface_basis_detail
|
||||
|
||||
template <class State>
|
||||
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<d::Face> faces;
|
||||
constexpr int attribute = Schema::template boundary_attribute<mean_field::utils::domain::StellarSurface>();
|
||||
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<int> 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<long double>(shape(i)) * face.coefficients(i, j);
|
||||
samples(sample, j) = static_cast<double>(value);
|
||||
weighted(sample, j) = std::sqrt(weight) * samples(sample, j);
|
||||
}
|
||||
}
|
||||
}
|
||||
result.physicalSurfaceArea = static_cast<double>(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<double>(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<double> pivots(k), relativePivots(k);
|
||||
for (int j = 0; j < k; ++j) {
|
||||
result.coefficientTransform(j, j) = 1.;
|
||||
const double originalNorm = std::sqrt(static_cast<double>(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<double>(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<double>(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<double>(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<double>(std::chrono::steady_clock::now() - start).count() << '\n';
|
||||
return result;
|
||||
}
|
||||
} // namespace experiment::plateau_diagnostics
|
||||
270
experiments/reduced_surface_diagnostics.hpp
Normal file
270
experiments/reduced_surface_diagnostics.hpp
Normal file
@@ -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 <memory>
|
||||
|
||||
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 <class State>
|
||||
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<ReducedSurfaceCoordinates> coordinates;
|
||||
if (!identity) {
|
||||
auto basis = MakeReducedSurfaceBasis(state, fem, maximumDegree, output);
|
||||
coordinates = std::make_unique<ReducedSurfaceCoordinates>(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<std::pair<const char *, const mfem::Vector *>>{
|
||||
{"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<State> 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
|
||||
527
experiments/reduced_surface_operators.hpp
Normal file
527
experiments/reduced_surface_operators.hpp
Normal file
@@ -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 <algorithm>
|
||||
#include <cmath>
|
||||
#include <concepts>
|
||||
#include <cstdint>
|
||||
#include <limits>
|
||||
#include <stdexcept>
|
||||
#include <string>
|
||||
#include <type_traits>
|
||||
#include <vector>
|
||||
|
||||
#include <mfem.hpp>
|
||||
#include <mpi.h>
|
||||
|
||||
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<long double>(matrix(i, column)) * value[i];
|
||||
return static_cast<double>(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<long double>(basis(i, j)) * basis(i, l);
|
||||
report.orthogonalityError = std::max(
|
||||
report.orthogonalityError, std::abs(static_cast<double>(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<long double>(basis(i, l)) * transform(l, j);
|
||||
inputNorm = std::hypot(inputNorm, input(i, j));
|
||||
errorNorm = std::hypot(errorNorm, static_cast<double>(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<long double>(basis(i, j)) * reduced[offset + j];
|
||||
full[offset + i] = static_cast<double>(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<long double>(basis(i, j)) * full[offset + i];
|
||||
reduced[offset + j] = static_cast<double>(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 <std::derived_from<mfem::Operator> 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 <class FullInverse>
|
||||
requires std::derived_from<FullInverse, mfem::Solver> &&
|
||||
mean_field::solver::LinearPreconditionerApplicationContractAvailable<FullInverse>
|
||||
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<FullInverse>;
|
||||
|
||||
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<double>::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
|
||||
Reference in New Issue
Block a user