From 9bd4af699650e2a6bbe795c618e6ef3fbcb48e12 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Wed, 16 Sep 2026 14:42:41 -0400 Subject: [PATCH] Record reproducible plateau campaigns and selected-direction probes --- experiments/plateau_control_cases.json | 47 +++++++ experiments/plateau_diagnostic.cpp | 40 ++++-- experiments/plateau_diagnostic_README.md | 50 +++++++- experiments/run_plateau_campaign.py | 154 +++++++++++++++++++++++ 4 files changed, 274 insertions(+), 17 deletions(-) create mode 100644 experiments/plateau_control_cases.json create mode 100644 experiments/run_plateau_campaign.py diff --git a/experiments/plateau_control_cases.json b/experiments/plateau_control_cases.json new file mode 100644 index 0000000..d933320 --- /dev/null +++ b/experiments/plateau_control_cases.json @@ -0,0 +1,47 @@ +[ + { + "id": "n1_J0p05_no_central_R40", + "arguments": [ + "--index", + "1", + "--angular-momentum", + ".05", + "--no-central", + "--restart", + "40", + "--max-newton", + "8" + ], + "notes": "Matched central-phase removal control; baseline solver policy." + }, + { + "id": "n1_J0p05_no_central_R160", + "arguments": [ + "--index", + "1", + "--angular-momentum", + ".05", + "--no-central", + "--restart", + "160", + "--max-newton", + "8" + ], + "notes": "Only Krylov restart differs from preceding model." + }, + { + "id": "n1p5_J0_no_central_R160", + "arguments": [ + "--index", + "1.5", + "--angular-momentum", + "0", + "--no-central", + "--restart", + "160", + "--max-newton", + "8" + ], + "notes": "Higher-index live control before initial-guess safeguard." + } +] diff --git a/experiments/plateau_diagnostic.cpp b/experiments/plateau_diagnostic.cpp index e3f1b14..9e3428a 100644 --- a/experiments/plateau_diagnostic.cpp +++ b/experiments/plateau_diagnostic.cpp @@ -23,7 +23,7 @@ namespace { using Clock = std::chrono::steady_clock; struct Options { double n = 1., angularMomentum = 0., atol = 0., rtol = 1e-8, linearTolerance = .03; - int maxNewton = 12, maxLinear = 200, restart = 40; + int maxNewton = 12, maxLinear = 200, restart = 40, probeRestart = 40; bool probes = false, noCentral = false; std::string mesh = "sandbox.smesh"; std::filesystem::path output, probeState; @@ -51,6 +51,8 @@ Options Parse(int argc, char **argv) { o.maxLinear = std::stoi(v()); else if (a == "--restart") o.restart = std::stoi(v()); + else if (a == "--probe-restart") + o.probeRestart = std::stoi(v()); else if (a == "--linear-tol") o.linearTolerance = std::stod(v()); else if (a == "--atol") @@ -68,6 +70,8 @@ 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.probeRestart != 40 && o.probeRestart != 160) + throw std::invalid_argument("--probe-restart must be 40 or 160"); if (std::filesystem::exists(o.output)) throw std::runtime_error("Output directory already exists"); std::filesystem::create_directories(o.output); @@ -130,7 +134,9 @@ void FrozenProbes(State &s, const fem::FEM &f, const Options &o, const mfem::Vec histories << "restart,tolerance,iteration,reported_norm,final\n"; auto blocks = File(o.output / "frozen_blocks.csv"); blocks << "restart,tolerance,block,F_norm,Jp_plus_F_norm,global_correction_norm\n"; - mfem::Vector best, baselineAction(base.Size()); + mfem::Vector selectedCorrection, baselineAction(base.Size()); + // These are two cold direct MFEM solves, not the configured production + // backend or any recycled Krylov history. for (int restart : {40, 160}) { mfem::FGMRESSolver linear(MPI_COMM_WORLD); solver::ResidualHistoryMonitor monitor; @@ -164,17 +170,18 @@ void FrozenProbes(State &s, const fem::FEM &f, const Options &o, const mfem::Vec << ',' << Norm(action, r.offset, r.size) << ',' << p.Norml2() << '\n'; blocks.flush(); std::cout << "Frozen restart=" << restart << " true=" << rel << " seconds=" << seconds << std::endl; - if (restart == 40) - best = p; + if (restart == o.probeRestart) + selectedCorrection = p; } - s.normalizedCorrection = best; - s.normalizedOperator->Mult(best, baselineAction); + std::cout << "Frozen probe direction restart=" << o.probeRestart << std::endl; + s.normalizedCorrection = selectedCorrection; + s.normalizedOperator->Mult(selectedCorrection, baselineAction); auto positivity = File(o.output / "enthalpy_step_bound.csv"); positivity << "metric,value\n"; for (const auto &[key, value] : experiment::plateau_diagnostics::MeasureEnthalpyStep(s, f, o.n, physicalBase)) positivity << key << ',' << value << '\n'; positivity.flush(); - SaveVector(o.output / "probe_normalized_correction.txt", best); + SaveVector(o.output / "probe_normalized_correction.txt", selectedCorrection); auto fd = File(o.output / "finite_differences.csv"); fd << "epsilon,admissible,block,absolute_error,relative_to_Jp,action_norm\n"; const auto safe = s.EstimateLargestSafeStepSize(1., .9); @@ -186,7 +193,7 @@ void FrozenProbes(State &s, const fem::FEM &f, const Options &o, const mfem::Vec for (double a : {1., .5, .1, .01, .001, .0001, .00001}) { const double epsilon = std::min(1., safe.stepSize) * a; s.trialNormalizedState = frozenState; - s.trialNormalizedState.Add(epsilon, best); + s.trialNormalizedState.Add(epsilon, selectedCorrection); const auto prep = s.PrepareTrial(); if (!prep) { trials << epsilon << ",0," << int(prep.error().reason) << ",nan,nan\n"; @@ -217,7 +224,9 @@ template int RunContext(Context &context, const Options &o, cons auto metadata = File(o.output / "metadata.txt"); metadata << "n=" << o.n << "\nK=" << scales.K << "\nrho_c=" << scales.rho << "\nJ=" << o.angularMomentum << "\nM=1\nR_reference=1\nxi1=" << scales.xi << "\nmu1=" << scales.mu << "\ncentral_phase=" << !o.noCentral - << "\nrestart=" << o.restart << "\nlinear_tol=" << o.linearTolerance << "\nmax_linear=" << o.maxLinear + << "\nrestart=" << o.restart << "\nprobe_restart=" << o.probeRestart + << "\nfrozen_probe_backend=direct_mfem_fgmres_cold" + << "\nlinear_tol=" << o.linearTolerance << "\nmax_linear=" << o.maxLinear << "\nmax_newton=" << o.maxNewton << "\natol=" << o.atol << "\nrtol=" << o.rtol << "\ncompiler=" << __VERSION__ << "\ncompiled=" << __DATE__ << ' ' << __TIME__ << '\n'; metadata.flush(); @@ -289,7 +298,8 @@ template int RunContext(Context &context, const Options &o, cons }); auto history = File(o.output / "newton.csv"); history << "iteration,accepted,F,alpha,trials,linear_status,linear_iterations,linear_true,seconds,linear_seconds," - "prepare_seconds,refresh_seconds,geometry_limited,geometry_safe_step,central_border,unbordered_F\n"; + "prepare_seconds,refresh_seconds,geometry_limited,geometry_safe_step,central_border,unbordered_F," + "linear_initial_residual,linear_rhs_norm,linear_initial_over_rhs\n"; auto blocklog = File(o.output / "blocks.csv"); blocklog << "iteration,stage,block,normalized_norm,physical_norm,max_abs_normalized\n"; auto triallog = File(o.output / "trials.csv"); @@ -327,6 +337,13 @@ template int RunContext(Context &context, const Options &o, cons unbordered(i) = e.normalizedResidual[i]; for (int i : centralRows) unbordered(i) -= residualFactors(i) * border; + const double linearInitial = e.linearSolve ? e.linearSolve->initialResidualNorm : NAN; + const double linearRhs = e.linearSolve ? e.linearSolve->rightHandSideNorm : NAN; + double linearInitialRelative = NAN; + if (linearRhs > 0.) + linearInitialRelative = linearInitial / linearRhs; + else if (linearRhs == 0. && std::isfinite(linearInitial)) + linearInitialRelative = linearInitial == 0. ? 0. : std::numeric_limits::infinity(); history << e.iteration << ',' << e.stepAccepted << ',' << e.residualNorm << ',' << e.acceptedStepLength << ',' << e.lineSearchTrials << ',' << (e.linearSolve ? int(e.linearSolve->status) : -1) << ',' << (e.linearSolve ? e.linearSolve->iterations : 0) << ',' @@ -335,7 +352,8 @@ template int RunContext(Context &context, const Options &o, cons << ',' << e.preconditionerRefreshSeconds << ',' << (e.geometryPreflight ? e.geometryPreflight->limitedByGeometry : false) << ',' << (e.geometryPreflight ? e.geometryPreflight->stepSize : NAN) << ',' << border << ',' - << unbordered.Norml2() << '\n'; + << unbordered.Norml2() << ',' << linearInitial << ',' << linearRhs << ',' + << linearInitialRelative << '\n'; history.flush(); std::cout << " alpha=" << e.acceptedStepLength << " F=" << e.residualNorm << " linear=" << (e.linearSolve ? e.linearSolve->relativeTrueResidualNorm : NAN) diff --git a/experiments/plateau_diagnostic_README.md b/experiments/plateau_diagnostic_README.md index 863f67d..ce9a6bc 100644 --- a/experiments/plateau_diagnostic_README.md +++ b/experiments/plateau_diagnostic_README.md @@ -14,7 +14,7 @@ Run a calibrated spherical family with `G=M=R_reference=1`: ```sh ./cmake-build-release-homebrew/plateau_diagnostic \ - --index 1.5 --angular-momentum 0 --max-newton 8 \ + --index 1.5 --angular-momentum 0 --no-central --max-newton 8 \ --output /absolute/path/to/new-results ``` @@ -28,6 +28,9 @@ analytic constants. Other indices use the production radial integrator, so this calibration is seed-consistent, not an independent exact reference. Angular momentum sets a physical target, not a fraction of critical rotation. `--no-central` omits the central phase for controlled comparisons; it changes the model specification. +The central phase remains enabled by default to reproduce earlier runs. Pass +`--no-central` for the n=1 and n=1.5 campaign models; retain the optional phase for +the n=3 homology-control experiment. Default solver controls are restart 40, 200 linear iterations, linear relative tolerance 0.03, Newton relative tolerance 1e-8 and absolute tolerance zero. Override @@ -35,13 +38,20 @@ with `--restart`, `--max-linear`, `--linear-tol`, `--rtol`, `--atol`, and `--max-newton`. Exit 0 means the requested Newton criterion was met, 3 means it was not met (or a seed-only inspection was requested), and 2 means setup/diagnostics failed. A nonconverged result still saves the retained accepted state and its -physical measurements. +physical measurements. `newton.csv` also records `linear_initial_residual`, +`linear_rhs_norm`, and `linear_initial_over_rhs` to measure the incoming correction +warm start. The initial residual is `||b - J p_initial||`; the ratio is zero for an +exactly solved zero RHS and infinity for a nonzero initial residual with zero RHS. +Missing solve reports produce NaN. ## Frozen investigations Append `--probes` to inspect the retained state after a solve. This compares cold FGMRES windows 40 and 160 at the same iteration budget and independently recomputes -both true residuals. It also exports normalization consistency checks, the scalar +both true residuals. These are direct `mfem::FGMRESSolver` calls with empty Krylov +history, not calls through the configured production linear backend; they do not +test a modified backend or its recycling behavior. It also exports normalization +consistency checks, the scalar Schur matrix and its coordinate labels, inverse-diagonal quality, and full-Jacobian actions on preconditioner-completed scalar directions. @@ -49,7 +59,8 @@ To inspect a previously saved physical root vector without repeating Newton: ```sh ./cmake-build-release-homebrew/plateau_diagnostic \ - --index 1.5 --angular-momentum 0 --max-newton 0 \ + --index 1.5 --angular-momentum 0 --no-central --max-newton 0 \ + --probe-restart 160 \ --probe-state /absolute/path/to/previous/accepted_state.txt \ --output /absolute/path/to/new-frozen-results ``` @@ -60,8 +71,12 @@ diagnostic hook; the context's accepted state is unchanged. This is a frozen inspection facility, not a Newton continuation/restart facility. The saved vector is dimension-checked, but does not contain a mesh/model fingerprint. -Finite-difference and nonlinear-trial probes always use the **window-40 direction**, -even if window 160 obtains a better correction. They compute the exact sampled +Set `--probe-restart 40` or `--probe-restart 160` to select the correction used for +positivity, finite-difference, and nonlinear-trial probes. The default is 40 for +reproducibility; both cold solves still run, and the choice does not change the +production Newton `--restart` setting. `metadata.txt` records `probe_restart` and +`frozen_probe_backend=direct_mfem_fgmres_cold`; the selected correction is saved as +`probe_normalized_correction.txt`. The probes compute the exact sampled enthalpy step bounds for both production closure and pressure-force quadrature. Every trial is restored and none is committed. Finite differences are one-sided; interpret their epsilon trend and absolute error, particularly for tiny Jacobian @@ -91,3 +106,26 @@ Quadrature samples are not a global positivity certificate. from the hydrostatic equation. A small bordered residual can otherwise hide a pointwise modification at the central vertex. Virial success and Newton stopping success are deliberately separate measurements. + +## Reproducible serial campaigns + +`run_plateau_campaign.py` requires a committed source revision and an already +configured CMake build directory. It rebuilds the diagnostic target, records the +source commit, build settings, binary and mesh hashes, and copies the executable +into the result directory before running cases serially. Every exact command, +exit code, log, state and compact summary is retained. A model exit code 3 is an +observed nonconvergence; campaign exit 0 means the observations completed, not +that every model converged. Output directories must be new. + +```sh +python3 experiments/run_plateau_campaign.py \ + --source . --build cmake-build-release-homebrew --mesh sandbox.smesh \ + --cases experiments/plateau_control_cases.json \ + --output /absolute/path/to/new-campaign +``` + +Use `--cmake` if CMake is not on PATH. Build products, raw states and full result +folders should remain outside source commits; retain compact result tables and +reproduction instructions with each hypothesis branch. The controls omit the +central phase at n=1 and n=1.5. A separate n=3 control retains it to fix homology; +the production model API does not insert it implicitly. diff --git a/experiments/run_plateau_campaign.py b/experiments/run_plateau_campaign.py new file mode 100644 index 0000000..68e85e0 --- /dev/null +++ b/experiments/run_plateau_campaign.py @@ -0,0 +1,154 @@ +#!/usr/bin/env python3 +"""Build a clean source revision and run a serial, fully recorded numerical campaign.""" +import argparse +import csv +import datetime +import hashlib +import json +import math +import os +from pathlib import Path +import platform +import shutil +import subprocess +import time + + +def digest(path): + result = hashlib.sha256() + with path.open('rb') as stream: + for chunk in iter(lambda: stream.read(1024 * 1024), b''): + result.update(chunk) + return result.hexdigest() + + +def utc(): + return datetime.datetime.now(datetime.timezone.utc).isoformat() + + +def json_safe(value): + if isinstance(value, float) and not math.isfinite(value): + return None + if isinstance(value, dict): + return {key: json_safe(item) for key, item in value.items()} + if isinstance(value, list): + return [json_safe(item) for item in value] + return value + + +def write_json(path, value): + path.write_text(json.dumps(json_safe(value), indent=2, allow_nan=False) + '\n') + + +def read_rows(path): + if not path.exists(): + return [] + with path.open() as stream: + return list(csv.DictReader(stream)) + + +def summarize(folder): + result = {'newton': read_rows(folder / 'newton.csv'), + 'frozen_linear': read_rows(folder / 'frozen_linear.csv'), + 'trials': read_rows(folder / 'trials.csv')} + meta = folder / 'metadata.txt' + if meta.exists(): + result['metadata'] = dict(line.split('=', 1) for line in meta.read_text().splitlines() if '=' in line) + for stage in ['initial', 'final']: + rows = read_rows(folder / f'{stage}_physical.csv') + result[stage + '_physical'] = { + q: {row['metric']: float(row['value']) for row in rows if row['quadrature'] == q} + for q in sorted({row['quadrature'] for row in rows})} + return result + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument('--source', required=True, type=Path) + parser.add_argument('--build', required=True, type=Path) + parser.add_argument('--mesh', required=True, type=Path) + parser.add_argument('--cases', required=True, type=Path) + parser.add_argument('--output', required=True, type=Path) + parser.add_argument('--cmake', default='cmake') + parser.add_argument('--jobs', type=int, default=2) + parser.add_argument('--timeout', type=float, default=1800) + args = parser.parse_args() + source, build, mesh, output = [p.resolve() for p in [args.source, args.build, args.mesh, args.output]] + if output.exists(): + parser.error('Output directory already exists; use a new directory for every run.') + def git(*words): + return subprocess.check_output(['git', *words], cwd=source, text=True).strip() + if git('status', '--porcelain', '--untracked-files=no'): + parser.error('Tracked source must be committed before running a campaign.') + cache = dict(line.split('=', 1) for line in (build / 'CMakeCache.txt').read_text().splitlines() + if '=' in line and not line.startswith(('#', '//'))) + if Path(cache['CMAKE_HOME_DIRECTORY:INTERNAL']).resolve() != source: + parser.error('Build directory belongs to a different source checkout.') + cases = json.loads(args.cases.read_text()) + seen = set() + for case in cases: + name = case['id'] + if not name or any(c not in 'abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ0123456789_-' for c in name): + parser.error('Case IDs may contain only letters, digits, underscore and hyphen.') + if name in seen: + parser.error('Duplicate case ID.') + seen.add(name) + if not isinstance(case['arguments'], list) or not all(isinstance(a, str) for a in case['arguments']): + parser.error('Case arguments must be a list of strings.') + if '--mesh' in case['arguments'] or '--output' in case['arguments']: + parser.error('Mesh and output are controlled by the campaign runner.') + output.mkdir(parents=True) + build_command = [args.cmake, '--build', str(build), '--target', 'plateau_diagnostic', '-j', str(args.jobs)] + provenance = {'started_utc': utc(), 'source': str(source), 'commit': git('rev-parse', 'HEAD'), + 'branch': git('branch', '--show-current'), 'platform': platform.platform(), + 'build_command': build_command, 'build_settings': cache, + 'mesh': str(mesh), 'mesh_sha256': digest(mesh), 'cases': cases, + 'runtime_environment': {'OMPI_MCA_btl': 'self,sm'}} + write_json(output / 'provenance.json', provenance) + shutil.copy2(args.cases, output / 'cases.json') + print('BUILD', provenance['commit'], flush=True) + with (output / 'build.log').open('w') as stream: + subprocess.run(build_command, cwd=source, stdout=stream, stderr=subprocess.STDOUT, check=True) + if git('rev-parse', 'HEAD') != provenance['commit'] or git('status', '--porcelain', '--untracked-files=no'): + raise RuntimeError('Source changed during build.') + # A private immutable executable prevents later branch/build operations from changing queued cases. + executable = output / 'plateau_diagnostic' + shutil.copy2(build / 'plateau_diagnostic', executable) + provenance['binary_sha256'] = digest(executable) + provenance['build_finished_utc'] = utc() + write_json(output / 'provenance.json', provenance) + summaries = [] + for case in cases: + name = case['id'] + folder = output / name + command = [str(executable), *case['arguments'], '--mesh', str(mesh), '--output', str(folder)] + state_files = {} + for i, word in enumerate(case['arguments']): + if word == '--probe-state': + path = Path(case['arguments'][i + 1]).resolve() + state_files[str(path)] = digest(path) + print('START', name, flush=True) + start = time.monotonic() + record = {'id': name, 'command': command, 'started_utc': utc(), 'input_state_hashes': state_files} + write_json(output / (name + '_run.json'), record) + with (output / (name + '.log')).open('w') as stream: + try: + run = subprocess.run(command, cwd=source, stdout=stream, stderr=subprocess.STDOUT, + env={**os.environ, 'OMPI_MCA_btl': 'self,sm'}, timeout=args.timeout) + record['exit_code'] = run.returncode + except subprocess.TimeoutExpired: + record['timed_out'] = True + record['exit_code'] = None + record.update(seconds=time.monotonic() - start, finished_utc=utc()) + write_json(output / (name + '_run.json'), record) + summaries.append({**record, **summarize(folder)}) + write_json(output / 'summary.json', summaries) + print('FINISH', name, record['exit_code'], round(record['seconds'], 2), flush=True) + if record['exit_code'] not in (0, 3): + raise RuntimeError('Case failed to execute; see its retained log and run record.') + # Exit 0 denotes completed observations; exit 3 from a model remains explicitly recorded as nonconvergence. + return 0 + + +if __name__ == '__main__': + raise SystemExit(main())