tboudreaux/MeanField
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Record reproducible plateau campaigns and selected-direction probes

Browse Source
This commit is contained in:
Emily Boudreaux
2026-09-16 14:42:41 -04:00
parent 80cc3d83d7
commit 9bd4af6996
4 changed files with 274 additions and 17 deletions
Download Patch File Download Diff File Expand all files Collapse all files

47
experiments/plateau_control_cases.json Normal file
View File

@@ -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."
}
]

40
experiments/plateau_diagnostic.cpp
View File

@@ -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 <class Context> 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 <class Context> 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 <class Context> 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<double>::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 <class Context> 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)

50
experiments/plateau_diagnostic_README.md
View File

@@ -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.

154
experiments/run_plateau_campaign.py Normal file
View File

@@ -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())