diff --git a/experiments/SHAPE_FIELD_WEAK_RESPONSE_2026-09-18.md b/experiments/SHAPE_FIELD_WEAK_RESPONSE_2026-09-18.md new file mode 100644 index 0000000..6e9ec1e --- /dev/null +++ b/experiments/SHAPE_FIELD_WEAK_RESPONSE_2026-09-18.md @@ -0,0 +1,65 @@ +# Shape–field weak response investigation — 2026-09-18 + +## Prespecified scope + +Follow up the large, nearly cancelling shape/field Newton direction for n=1.5, +J=0 without a central-density constraint. This study distinguishes within-element +physical field transport from surface variation and tests the specific hypothesis +that the Q3 displacement/Q4 reference geometry mismatch contributes to weak +response. Cancellation by itself is not evidence of a gauge or bad element orders. + +The baseline branch is `research/shape-field-weak-response-20260918`, based on +`db854c0f7cdf5137e377273499f0ac5b6053b752`. The isolated order variant will be +`research/shape-field-displacement-q4-20260918`. Preserve original checkout, +immutable executables, source revisions, mesh/input hashes and serial run records. + +## Observation and controls + +The diagnostic prepares a frozen state with fixed normalization, independently +verifies its full-J linear defect, saves every value-block Jacobian contribution, +and checks restoration of accepted buffers and operator action. Three small raw +Taylor trials test the same original operator. Physical observations use independent +grid functions and the prepared displacement; no accepted state is modified. + +At common quadrature orders 10 and 14, record physical scalar pullback changes, +transport `grad(f).v`, and Eulerian changes `delta f - grad(f).v`. Density and +potential gradients are broken/local: facet jumps and RT Piola variations are not +included, so these observations cannot certify a pure relabelling nullspace. +Surface radial parameters are unweighted nodal samples. An independent affine +identity projection checks ambient displacement representation, separately from +the actual r-squared radial extension. Small-span SVDs describe only the sampled +block directions, never the full Jacobian rank. + +Baseline frozen cases: the recorded n=1.5 seed and late directions at eta=0.03, +and a fresh cold solve at the six-update coupled endpoint. Then a fresh six-update +baseline with gravity-Schur inverse, gravity-curvature globalization, restart 160, +maximum 200 Krylov iterations, eta=0.03, and a fresh final direction. Match this +live run with the Q4 variant from the same analytic seed, including a Q4 seed +probe. Cases, builds and expensive audits run serially. A failed linear gate or +preparation remains a failure; do not silently relax controls. + +## Order intervention (conditional on baseline observation) + +Only change `Displacement::vectorOrder` from `3 + uniform increment` to +`4 + uniform increment` on the variant branch; keep uniform increment zero. +The loaded reference geometry remains Q4, density/potential L2 Q2, gravity RT +index 2, and enthalpy H1 Q3. The shared displacement collection means the scalar +surface trace also becomes Q4. This introduces shape modes and changes derived +quadrature orders and cost; it is not a pure interpolation-error intervention. +Document exact source line, date, rationale, DOFs, and compiled quadrature changes. +Use a separate build directory and new analytic seeds: existing coefficient +snapshots cannot be transferred unchanged between orders. + +## Decision criteria + +Compare verified linear accuracy, accepted step lengths, within-variant growth +of weak-direction indicators, common physical norms, virials, positivity/EOS +metrics and elapsed cost. Across orders, coefficient norms and discrete residuals +are not physical error measures. Q10/Q14 agreement diagnoses integration +sensitivity, not convergence certification. Virial <=1e-5 remains required but +cannot replace internal pointwise checks. Six updates are a bounded mechanistic +comparison, not evidence of high-n/high-rotation convergence. + +## Results + +Pending recorded runs. No production default has been changed on the baseline. diff --git a/experiments/plateau_diagnostic.cpp b/experiments/plateau_diagnostic.cpp index 2f8e3ec..935d124 100644 --- a/experiments/plateau_diagnostic.cpp +++ b/experiments/plateau_diagnostic.cpp @@ -30,6 +30,7 @@ import mean_field; #include "plateau_gravity_application_diagnostics.hpp" #include "plateau_positivity_strategy_diagnostics.hpp" #include "plateau_polynomial_certificate.hpp" +#include "shape_field_weak_diagnostics.hpp" using namespace mean_field; namespace { @@ -40,6 +41,7 @@ struct Options { bool probes = false, noCentral = false; bool strategyProbe = false, positivityProbe = false, positivityCertificate = false; bool preconditionerProbe = false, nonlinearStepProbe = false; + bool weakShapeProbe = false, weakShapeFinal = false; bool coupledTrialProbe = false, coupledSeedControl = false; bool dampedCoupledProbe = false; std::string coupledSearch = "near"; @@ -113,6 +115,10 @@ Options Parse(int argc, char **argv) { o.preconditionerProbe = true; else if (a == "--nonlinear-step-probe") o.nonlinearStepProbe = true; + else if (a == "--weak-shape-probe") + o.weakShapeProbe = true; + else if (a == "--weak-shape-final") + o.weakShapeFinal = true; else if (a == "--damped-coupled-probe") o.dampedCoupledProbe = true; else if (a == "--coupled-search") @@ -171,6 +177,16 @@ 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.weakShapeProbe || o.weakShapeFinal) { + if ((o.weakShapeProbe && o.weakShapeFinal) || o.probes || o.pressureProbe || o.shapeCurvature || + o.surfaceTrace || o.strategyProbe || o.positivityProbe || o.positivityCertificate || + o.preconditionerProbe || o.nonlinearStepProbe || o.coupledTrialProbe || o.dampedCoupledProbe || + !o.replayCandidate.empty() || !o.captureIterations.empty() || + (o.weakShapeProbe && (o.maxNewton != 0 || o.strategy != "baseline")) || + (o.weakShapeFinal && (o.maxNewton == 0 || !o.probeState.empty() || !o.probeDirection.empty())) || + !std::isfinite(o.linearTolerance) || o.linearTolerance <= 0.) + throw std::invalid_argument("Weak-shape diagnostics require a standalone frozen probe or a final live endpoint probe"); + } if (o.probeRestart != 40 && o.probeRestart != 160) throw std::invalid_argument("--probe-restart must be 40 or 160"); if ((o.pressureProbe || o.shapeCurvature || o.surfaceTrace || o.strategyProbe || o.positivityProbe || o.coupledTrialProbe || o.dampedCoupledProbe) && (o.probeState.empty() || o.probeDirection.empty())) @@ -218,8 +234,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.probeDirection.empty() && !o.coupledTrialProbe && !o.dampedCoupledProbe) || + 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.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" && @@ -513,7 +529,8 @@ 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 ? "prepared_fgmres_cold_with_independent_true_gate" : + (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") << "\nlinear_tol=" << o.linearTolerance << "\nmax_linear=" << o.maxLinear @@ -522,6 +539,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 + << "\nweak_shape_probe=" << o.weakShapeProbe << "\nweak_shape_final=" << o.weakShapeFinal << "\ndamped_coupled_probe=" << o.dampedCoupledProbe << "\ncoupled_search=" << o.coupledSearch << "\ncoupled_beta_steps=" << o.coupledBetaSteps @@ -541,7 +559,29 @@ int RunContext(Context &context, const Options &o, const Scales &scales) { << "\ncontext_gravity_schur_relaxation=" << (contextUsesGravitySchur ? o.schurRelaxation : 0.0) << "\nschur_correction_sweeps=" << o.schurCorrectionSweeps << "\nschur_relaxation=" << o.schurRelaxation << "\ncompiler=" << __VERSION__ << "\ncompiled=" << __DATE__ << ' ' << __TIME__ << '\n'; + if (o.weakShapeProbe || o.weakShapeFinal) + metadata << "weak_shape_direction_space=normalized_current_factors\n" + << "weak_shape_operator=full_monolithic\n" + << "weak_shape_trial_policy=raw_prepare_only\n"; metadata.flush(); + 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'; + metadata.flush(); + solver::detail::StellarEquilibriumContextDiagnostics::WithState(context, [&](auto &s, const fem::FEM &f) { + mfem::Vector base(s.acceptedNormalizedState), direction; + if (!o.probeState.empty()) { + const auto physical = LoadVector(o.probeState, s.AcceptedPhysicalState().Size()); + s.normalizedOperator->NormalizeState(physical, base); + } + if (!o.probeDirection.empty()) direction = LoadVector(o.probeDirection, base.Size()); + experiment::plateau_diagnostics::ProbeWeakShapeField(s, f, base, + o.probeDirection.empty() ? nullptr : &direction, o.output, { + .linear = {.relativeTolerance = o.linearTolerance, .absoluteTolerance = 0., + .maximumIterations = o.maxLinear}, .restartLength = o.restart}); + }); + return 0; + } if (o.dampedCoupledProbe) { metadata << "coupled_trial_input=physical_snapshot_and_supplied_normalized_direction\n" << "coupled_trial_admissibility=sampled_enthalpy_without_repair\n" @@ -979,6 +1019,10 @@ int RunContext(Context &context, const Options &o, const Scales &scales) { experiment::plateau_diagnostics::InspectPreconditioner(s, o.output); FrozenProbes(s, f, o); } + if (o.weakShapeFinal) + experiment::plateau_diagnostics::ProbeWeakShapeField(s, f, s.acceptedNormalizedState, nullptr, o.output, { + .linear = {.relativeTolerance = o.linearTolerance, .absoluteTolerance = 0., + .maximumIterations = o.maxLinear}, .restartLength = o.restart}); }); return converged ? 0 : 3; } diff --git a/experiments/shape_field_physical_diagnostics.hpp b/experiments/shape_field_physical_diagnostics.hpp new file mode 100644 index 0000000..4884606 --- /dev/null +++ b/experiments/shape_field_physical_diagnostics.hpp @@ -0,0 +1,456 @@ +#pragma once + +// Include after import mean_field. These observations use the already prepared +// base and private grid functions; they do not prepare a trial or alter a solver +// buffer. Scalar Eulerian variations are local, broken derivatives for DG +// fields. They do not contain jump distributions and cannot establish a gauge. +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +namespace experiment::plateau_diagnostics { +namespace shape_field_physical_detail { +using DomainSchema = mean_field::utils::domain::CoreEnvelopeVacuumDomainSchema; + +// DomainMapper caches reference geometry even through its const interface. +// A private instance prevents the observation from warming the solver's cache. +// Its exterior policy deliberately has no numerical behavior: all observations +// below must be interior under both this map and the prepared production map. +class RejectExteriorMap final : public mean_field::mapping::compactification::ExteriorDomainMap { +public: + mean_field::mapping::MappingStatus Evaluate( + const mean_field::mapping::compactification::ExteriorMapInput &, + mean_field::mapping::compactification::ExteriorMapResult &) const override { + throw std::logic_error("Shape-field physical diagnostic cannot evaluate the exterior domain."); + } + mean_field::mapping::MappingStatus EvaluateVariation( + const mean_field::mapping::compactification::ExteriorMapInput &, + const mean_field::mapping::compactification::ExteriorMapResult &, + const mean_field::mapping::compactification::ExteriorMapDirection &, + mean_field::mapping::compactification::ExteriorMapVariation &) const override { + throw std::logic_error("Shape-field physical diagnostic cannot vary the exterior domain."); + } + std::string_view GetName() const noexcept override { return "diagnostic_interior_only"; } +}; + +inline void Require(bool condition, const std::string &message) { + if (!condition) throw std::runtime_error("Shape-field physical diagnostic: " + message); +} + +inline bool Finite(const mfem::Vector &vector) { + for (int i = 0; i < vector.Size(); ++i) + if (!std::isfinite(vector(i))) return false; + return true; +} + +inline bool Identical(const mfem::Vector &left, const mfem::Vector &right) { + return left.Size() == right.Size() && + (left.Size() == 0 || std::memcmp(left.GetData(), right.GetData(), + sizeof(mfem::real_t) * left.Size()) == 0); +} + +inline std::string Quote(std::string_view text) { + std::string result = "\""; + for (const char ch : text) { + if (ch == '"') result += '"'; + result += ch; + } + return result + '"'; +} + +inline std::ofstream Output(const std::filesystem::path &directory, std::string_view suffix) { + const auto path = directory / ("weak_shape_physical_" + std::string(suffix) + ".csv"); + Require(!std::filesystem::exists(path), "output already exists: " + path.string()); + std::ofstream output(path); + output.exceptions(std::ios::failbit | std::ios::badbit); + output << std::setprecision(std::numeric_limits::max_digits10); + return output; +} + +struct Block { int offset; int size; }; + +template +Block RequireBlock(const Blocks &blocks, std::string_view name, int size) { + for (const auto &block : blocks) { + if (block.stableId != name) continue; + Require(block.offset >= 0 && block.size >= 0 && block.offset <= size - block.size, + "invalid root block " + std::string(name)); + return {block.offset, block.size}; + } + throw std::invalid_argument("Shape-field physical diagnostic requires root block " + std::string(name)); +} + +template +double Scatter(const Blocks &blocks, std::string_view name, const mfem::Vector &values, + mfem::ParGridFunction &field) { + const auto block = RequireBlock(blocks, name, values.Size()); + const auto adapter = mean_field::field::make_field_dof_grid_function_adapter( + *field.ParFESpace()); + Require(block.size == adapter.dof_map().reduced_size(), "field adapter size disagrees with " + std::string(name)); + mfem::Vector reduced(block.size); + for (int i = 0; i < block.size; ++i) reduced(i) = values(block.offset + i); + adapter.scatter(reduced, field); + const auto gathered = adapter.gather(field); + double maximumError = 0.; + for (int i = 0; i < block.size; ++i) + maximumError = std::max(maximumError, std::abs(gathered(i) - reduced(i))); + Require(Finite(field) && std::isfinite(maximumError), "nonfinite field reconstruction"); + return maximumError; +} + +// Every entry is a raw additive physical-volume integral. Keeping the signed +// cross term lets an independent reader verify E^2 = D^2 + T^2 - 2 D.T. +struct ScalarIntegrals { + long double base{}, pullback{}, transport{}, eulerian{}; + long double baseSquared{}, pullbackSquared{}, transportSquared{}, eulerianSquared{}, cross{}; + + void Add(double weight, double f, double delta, double advection) { + const long double w = weight, b = f, d = delta, t = advection, e = d - t; + base += w * b; + pullback += w * d; + transport += w * t; + eulerian += w * e; + baseSquared += w * b * b; + pullbackSquared += w * d * d; + transportSquared += w * t * t; + eulerianSquared += w * e * e; + cross += w * d * t; + } + + ScalarIntegrals &operator+=(const ScalarIntegrals &other) { + base += other.base; pullback += other.pullback; transport += other.transport; eulerian += other.eulerian; + baseSquared += other.baseSquared; pullbackSquared += other.pullbackSquared; + transportSquared += other.transportSquared; eulerianSquared += other.eulerianSquared; cross += other.cross; + return *this; + } + + void Write(std::ostream &output) const { + for (const auto value : {base, pullback, transport, eulerian, baseSquared, pullbackSquared, + transportSquared, eulerianSquared, cross}) + Require(std::isfinite(value), "nonfinite integrated scalar quantity"); + output << ',' << base << ',' << pullback << ',' << transport << ',' << eulerian << ',' << baseSquared + << ',' << pullbackSquared << ',' << transportSquared << ',' << eulerianSquared << ',' << cross; + } +}; + +inline constexpr std::string_view scalarColumns = + ",base_integral,pullback_integral,transport_integral,eulerian_integral,base_squared_integral," + "pullback_squared_integral,transport_squared_integral,eulerian_squared_integral,pullback_transport_cross_integral"; +} // namespace shape_field_physical_detail + +template +void MeasureShapeFieldPhysical(State &state, const mean_field::fem::FEM &fem, + const mfem::Vector &physicalBase, const mfem::Vector &physicalDirection, + const std::filesystem::path &output) { + namespace d = shape_field_physical_detail; + using namespace mean_field; + const auto &problem = state.Problem(); + int ranks = 0; + MPI_Comm_size(problem.GetCommunicator(), &ranks); + d::Require(ranks == 1, "exactly one MPI rank is required"); + d::Require(fem.mesh && fem.displacementFes && fem.surfaceDeformationFes && fem.densityFes && + fem.enthalpyFes && fem.gravityPotentialFes && fem.gravityFluxFes && fem.compactificationCoordinate, + "incomplete finite-element spaces"); + d::Require(fem.mesh->SpaceDimension() == 3 && fem.mesh->Dimension() == 3, + "three-dimensional volume elements are required"); + d::Require(state.normalizedOperator->IsPrepared() && problem.GetPhysicalOperator().IsPrepared(), + "the base operator is not prepared"); + d::Require(physicalBase.Size() == problem.StateSize() && physicalDirection.Size() == physicalBase.Size() && + d::Finite(physicalBase) && d::Finite(physicalDirection) && + d::Identical(physicalBase, state.normalizedOperator->GetPhysicalState()), + "the finite supplied base must exactly equal the prepared physical state"); + + // Capture what a caller could otherwise mistake for a silent preparation. + // The only production action below is the immutable deformation Jacobian. + 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), + 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()); + mfem::Vector volumeDirection(fem.displacementFes->GetTrueVSize()); + problem.BuildVolumeDisplacementDirection(physicalDirection, volumeDirection); + d::Require(volumeDirection.Size() == fem.displacementFes->GetTrueVSize() && d::Finite(volumeDirection), + "invalid generated volume direction"); + + mfem::ParGridFunction rho(fem.densityFes.get()), drho(fem.densityFes.get()), + h(fem.enthalpyFes.get()), dh(fem.enthalpyFes.get()), + phi(fem.gravityPotentialFes.get()), dphi(fem.gravityPotentialFes.get()), + displacement(fem.displacementFes.get()), velocity(fem.displacementFes.get()), + affine(fem.displacementFes.get()); + const auto blocks = problem.GetManifest().valueBlocks(); + const std::array scatterErrors{ + d::Scatter(blocks, "density", physicalBase, rho), + d::Scatter(blocks, "density", physicalDirection, drho), + d::Scatter(blocks, "specific_enthalpy", physicalBase, h), + d::Scatter(blocks, "specific_enthalpy", physicalDirection, dh), + d::Scatter(blocks, "gravity_potential", physicalBase, phi), + d::Scatter(blocks, "gravity_potential", physicalDirection, dphi)}; + displacement.SetFromTrueDofs(baseDisplacement); + velocity.SetFromTrueDofs(volumeDirection); + mfem::VectorFunctionCoefficient identity(3, [](const mfem::Vector &x, mfem::Vector &value) { value = x; }); + affine.ProjectCoefficient(identity); + d::Require(d::Finite(displacement) && d::Finite(velocity) && d::Finite(affine), "nonfinite private displacement field"); + const auto &productionMapper = problem.GetDiscretization().domainMapper(); + const auto productionCacheBefore = productionMapper.GetGeometryCacheStatistics(); + mapping::DomainMapper privateMapper( + utils::DomainMapperOptions{.dimension = 3, + .vacuum_element_attribute = d::DomainSchema::template material_attribute()}, + std::make_unique()); + privateMapper.SetGeometryCacheEnabled(false); + mapping::GridFunctionMappingEvaluator mappingEvaluator(privateMapper, displacement, *fem.compactificationCoordinate); + mappingEvaluator.InvalidateCache(); + + // The default surface descriptor erases its center. Checking both actual + // displacement traces prevents a different radial prescription from being + // silently interpreted in the zero-center basis used for this CSV. + const auto &domain = physical.GetDomainDeformation(); + mfem::Vector center(3); center = 0.; + const deformation::NodalRadialSurface prescription(std::move(center)); + d::Require(domain.descriptor().surfaceDeformation == prescription.descriptor() && + fem.surfaceDeformationFes->GetVDim() == 1 && fem.displacementFes->GetVDim() == 3 && + fem.displacementFes->GetOrdering() == mfem::Ordering::byNODES && + fem.displacementFes->FEColl() == fem.surfaceDeformationFes->FEColl() && + fem.displacementFes->GetTrueVSize() == 3 * fem.surfaceDeformationFes->GetTrueVSize(), + "the default nodal radial surface and its vector companion are required"); + const auto radial = deformation::compileSurfaceDeformationPrescription(prescription, + deformation::SurfaceDeformationCompilationContext{ + *fem.surfaceDeformationFes, field::make_stellar_surface_scalar_dof_map(*fem.surfaceDeformationFes)}); + const auto shape = d::RequireBlock(blocks, "surface_deformation", physicalBase.Size()); + d::Require(radial.parameterCount() == shape.size && shape.size == domain.parameterCount(), + "surface parameter counts disagree"); + const int scalarSize = fem.surfaceDeformationFes->GetTrueVSize(); + double traceError = 0., traceScale = 1.; + for (int i = 0; i < shape.size; ++i) + for (int component = 0; component < 3; ++component) { + const int dof = radial.surfaceDofMap().volume_true_dof(i) + component * scalarSize; + const double direction = radial.radialDirection(i, component); + const double expectedBase = physicalBase(shape.offset + i) * direction; + const double expectedDirection = physicalDirection(shape.offset + i) * direction; + traceError = std::max({traceError, std::abs(baseDisplacement(dof) - expectedBase), + std::abs(volumeDirection(dof) - expectedDirection)}); + traceScale = std::max({traceScale, std::abs(expectedBase), std::abs(expectedDirection)}); + } + d::Require(std::isfinite(traceError) && traceError <= 512. * std::numeric_limits::epsilon() * traceScale, + "actual generated surface trace disagrees with the reconstructed radial basis"); + + std::filesystem::create_directories(output); + auto metadata = d::Output(output, "metadata"); + auto elements = d::Output(output, "elements"); + auto geometry = d::Output(output, "geometry"); + auto affineOutput = d::Output(output, "affine"); + auto surface = d::Output(output, "surface_nodes"); + auto summary = d::Output(output, "summary"); + metadata << "key,value\n"; + const auto meta = [&](std::string_view key, std::string_view value) { + metadata << d::Quote(key) << ',' << d::Quote(value) << '\n'; + }; + meta("schema", "weak_shape_physical_v1"); + meta("scope", "stellar interior only; prepared base; physical direction supplied; no trial preparation"); + meta("scalar_formula", "delta_f_E = delta_f_at_fixed_reference - grad_physical_f dot volume_displacement_direction"); + meta("dg_derivatives", "density and potential gradients are element-local broken derivatives; jumps excluded"); + meta("interpretation", "local shape-field compensation diagnostic; no nullspace or pure-gauge certification"); + meta("quadrature", "fixed integration-rule orders 10 and 14, independent of field order; compare for integration sensitivity"); + meta("radius", "Euclidean distance from zero; per-element min/max are quadrature-sample extrema"); + meta("surface_samples", "nodal values without area weights; not harmonic coefficients or surface integrals"); + meta("affine_test", "projection of reference physical coordinate x into displacement space; independent of the actual domain lift"); + meta("affine_weights", "reference physical mesh volume; gradient is with respect to reference physical coordinates"); + meta("physical_weights", "integration_point.weight * element_transformation.Weight * mapping_determinant"); + meta("mapping", "private uncached DomainMapper; same interior x+u implementation; exterior evaluation forbidden"); + meta("scalar_integrals", "raw additive integrals; transport = grad_physical_f dot v; no distributional interface contribution"); + meta("gravity_vector", "not evaluated; RT Piola variation is outside this scalar diagnostic"); + meta("reference_mesh_geometry_collection", fem.mesh->GetNodes() + ? fem.mesh->GetNodes()->FESpace()->FEColl()->Name() : "linear_mesh_without_nodal_grid_function"); + const std::array, 6> spaces{{ + {"density", fem.densityFes.get()}, {"specific_enthalpy", fem.enthalpyFes.get()}, + {"gravity_potential", fem.gravityPotentialFes.get()}, {"displacement", fem.displacementFes.get()}, + {"surface_deformation", fem.surfaceDeformationFes.get()}, {"gravity_gradient", fem.gravityFluxFes.get()}}}; + for (const auto &[name, space] : spaces) { + meta(std::string(name) + "_collection", space->FEColl()->Name()); + meta(std::string(name) + "_true_size", std::to_string(space->GetTrueVSize())); + } + meta("density_family_order", std::to_string(field::Density::Scalar::familyOrder)); + meta("enthalpy_family_order", std::to_string(field::Enthalpy::Scalar::familyOrder)); + meta("potential_family_order", std::to_string(field::Gravity::Potential::familyOrder)); + meta("gravity_flux_family_order", std::to_string(field::Gravity::Flux::familyOrder)); + meta("displacement_family_order", std::to_string(field::Displacement::Vector::familyOrder)); + + elements << "quadrature,element,attribute,geometry,field,gradient_kind,points,physical_volume" + << d::scalarColumns << '\n'; + geometry << "quadrature,element,attribute,geometry,points,reference_mesh_geometry_order,density_element_order,enthalpy_element_order," + "potential_element_order,displacement_element_order,reference_radius_min,reference_radius_max," + "physical_radius_min,physical_radius_max,physical_volume,base_displacement_squared_integral," + "direction_displacement_squared_integral\n"; + affineOutput << "quadrature,element,points,reference_volume,value_error_squared_integral," + "reference_value_squared_integral,gradient_error_squared_integral,reference_gradient_squared_integral\n"; + surface << "parameter,volume_scalar_true_dof,unit_radial_x,unit_radial_y,unit_radial_z,reference_radius," + "physical_radius,base_parameter,physical_direction_parameter\n"; + for (int i = 0; i < shape.size; ++i) { + const double radius = radial.referenceRadius(i), base = physicalBase(shape.offset + i); + d::Require(std::isfinite(radius) && radius > 0. && radius + base > 0., "invalid radial surface radius"); + surface << i << ',' << radial.surfaceDofMap().volume_true_dof(i); + for (int component = 0; component < 3; ++component) surface << ',' << radial.radialDirection(i, component); + surface << ',' << radius << ',' << radius + base << ',' << base << ',' << physicalDirection(shape.offset + i) << '\n'; + } + summary << "quadrature,field,metric,value\n"; + const auto measure = [&](int q, std::string_view fieldName, std::string_view metric, long double value) { + d::Require(std::isfinite(value), "nonfinite summary " + std::string(metric)); + summary << q << ',' << fieldName << ',' << metric << ',' << value << '\n'; + }; + const std::array names{"density", "specific_enthalpy", "gravity_potential"}; + const std::array fields{&rho, &h, &phi}, directions{&drho, &dh, &dphi}; + mfem::Vector referenceGradient(3), physicalGradient(3), u(3), v(3), x(3), affineValue(3); + mfem::DenseMatrix affineGradient(3); + mapping::MappingPointContext mapped; + for (const int quadrature : {10, 14}) { + std::array totals; + long double totalVolume = 0., totalU2 = 0., totalV2 = 0., affineVolume = 0., affineError2 = 0., + affineReference2 = 0., affineGradientError2 = 0.; + int stellarElements = 0; + for (int element = 0; element < fem.mesh->GetNE(); ++element) { + const int attribute = fem.mesh->GetAttribute(element); + if (!d::DomainSchema::template attribute_belongs_to(attribute)) continue; + ++stellarElements; + auto *transformation = fem.mesh->GetElementTransformation(element); + d::Require(!productionMapper.IsCompactifiedElement(*transformation), + "production mapper classifies a stellar element as compactified"); + const auto &rule = mfem::IntRules.Get(transformation->GetGeometryType(), quadrature); + std::array integrals; + long double volume = 0., u2 = 0., v2 = 0., referenceVolume = 0., affineE2 = 0., affineR2 = 0., affineGE2 = 0.; + double rmin = std::numeric_limits::infinity(), rmax = 0., + physicalRmin = std::numeric_limits::infinity(), physicalRmax = 0.; + for (int point = 0; point < rule.GetNPoints(); ++point) { + const auto &ip = rule.IntPoint(point); + d::Require(mappingEvaluator.EvaluatePoint(*transformation, ip, mapped) == mapping::MappingStatus::valid && + !mapped.compactified && std::isfinite(mapped.mapping_determinant) && mapped.mapping_determinant > 0., + "invalid or compactified stellar mapping"); + transformation->SetIntPoint(&ip); + const double referenceWeight = ip.weight * transformation->Weight(); + const double weight = referenceWeight * mapped.mapping_determinant; + d::Require(std::isfinite(referenceWeight) && referenceWeight > 0. && std::isfinite(weight) && weight > 0., + "invalid quadrature weight"); + displacement.GetVectorValue(element, ip, u); + velocity.GetVectorValue(element, ip, v); + d::Require(d::Finite(u) && d::Finite(v) && d::Finite(mapped.reference_position) && d::Finite(mapped.physical_position), + "nonfinite sampled displacement or position"); + const double r = mapped.reference_position.Norml2(), physicalR = mapped.physical_position.Norml2(); + d::Require(std::isfinite(r) && std::isfinite(physicalR), "overflowing sampled radius"); + rmin = std::min(rmin, r); rmax = std::max(rmax, r); + physicalRmin = std::min(physicalRmin, physicalR); physicalRmax = std::max(physicalRmax, physicalR); + volume += weight; + for (int component = 0; component < 3; ++component) { + u2 += static_cast(weight) * u(component) * u(component); + v2 += static_cast(weight) * v(component) * v(component); + } + for (std::size_t fieldIndex = 0; fieldIndex < fields.size(); ++fieldIndex) { + const double base = fields[fieldIndex]->GetValue(element, ip); + const double delta = directions[fieldIndex]->GetValue(element, ip); + transformation->SetIntPoint(&ip); + fields[fieldIndex]->GetGradient(*transformation, referenceGradient); + mapped.inverse_mapping_jacobian.MultTranspose(referenceGradient, physicalGradient); + const double transport = physicalGradient * v; + d::Require(std::isfinite(base) && std::isfinite(delta) && std::isfinite(transport) && d::Finite(physicalGradient), + "nonfinite scalar or broken gradient"); + integrals[fieldIndex].Add(weight, base, delta, transport); + } + transformation->SetIntPoint(&ip); + transformation->Transform(ip, x); + affine.GetVectorValue(element, ip, affineValue); + affine.GetVectorGradient(*transformation, affineGradient); + d::Require(d::Finite(x) && d::Finite(affineValue), "nonfinite affine projection sample"); + referenceVolume += referenceWeight; + for (int row = 0; row < 3; ++row) { + const long double error = static_cast(affineValue(row)) - x(row); + affineE2 += referenceWeight * error * error; + affineR2 += static_cast(referenceWeight) * x(row) * x(row); + for (int column = 0; column < 3; ++column) { + const long double gradientError = affineGradient(row, column) - (row == column ? 1. : 0.); + d::Require(std::isfinite(gradientError), "nonfinite affine projection gradient"); + affineGE2 += referenceWeight * gradientError * gradientError; + } + } + } + for (std::size_t fieldIndex = 0; fieldIndex < fields.size(); ++fieldIndex) { + elements << quadrature << ',' << element << ',' << attribute << ',' << int(transformation->GetGeometryType()) + << ',' << names[fieldIndex] << ',' << (fieldIndex == 1 ? "continuous_H1" : "broken_DG") + << ',' << rule.GetNPoints() << ',' << volume; + integrals[fieldIndex].Write(elements); + elements << '\n'; + totals[fieldIndex] += integrals[fieldIndex]; + } + geometry << quadrature << ',' << element << ',' << attribute << ',' << int(transformation->GetGeometryType()) + << ',' << rule.GetNPoints() << ',' << (fem.mesh->GetNodes() + ? fem.mesh->GetNodes()->FESpace()->GetElementOrder(element) : 1) + << ',' << fem.densityFes->GetElementOrder(element) + << ',' << fem.enthalpyFes->GetElementOrder(element) << ',' << fem.gravityPotentialFes->GetElementOrder(element) + << ',' << fem.displacementFes->GetElementOrder(element) << ',' << rmin << ',' << rmax << ',' << physicalRmin + << ',' << physicalRmax << ',' << volume << ',' << u2 << ',' << v2 << '\n'; + affineOutput << quadrature << ',' << element << ',' << rule.GetNPoints() << ',' << referenceVolume << ',' + << affineE2 << ',' << affineR2 << ',' << affineGE2 << ',' << 3. * referenceVolume << '\n'; + totalVolume += volume; totalU2 += u2; totalV2 += v2; + affineVolume += referenceVolume; affineError2 += affineE2; affineReference2 += affineR2; affineGradientError2 += affineGE2; + } + d::Require(stellarElements > 0 && totalVolume > 0. && affineVolume > 0. && affineReference2 > 0., "empty stellar integrals"); + for (std::size_t i = 0; i < totals.size(); ++i) { + const auto &integral = totals[i]; + measure(quadrature, names[i], "base_squared_integral", integral.baseSquared); + measure(quadrature, names[i], "pullback_squared_integral", integral.pullbackSquared); + measure(quadrature, names[i], "transport_squared_integral", integral.transportSquared); + measure(quadrature, names[i], "eulerian_squared_integral", integral.eulerianSquared); + measure(quadrature, names[i], "pullback_transport_cross_integral", integral.cross); + measure(quadrature, names[i], "pullback_l2_norm", std::sqrt(integral.pullbackSquared)); + measure(quadrature, names[i], "transport_l2_norm", std::sqrt(integral.transportSquared)); + measure(quadrature, names[i], "eulerian_l2_norm", std::sqrt(integral.eulerianSquared)); + measure(quadrature, names[i], "squared_norm_identity_error", + integral.eulerianSquared - integral.pullbackSquared - integral.transportSquared + 2. * integral.cross); + } + measure(quadrature, "geometry", "stellar_elements", stellarElements); + measure(quadrature, "geometry", "physical_volume", totalVolume); + measure(quadrature, "geometry", "base_displacement_squared_integral", totalU2); + measure(quadrature, "geometry", "direction_displacement_squared_integral", totalV2); + measure(quadrature, "affine", "reference_volume", affineVolume); + measure(quadrature, "affine", "value_error_squared_integral", affineError2); + measure(quadrature, "affine", "reference_value_squared_integral", affineReference2); + measure(quadrature, "affine", "gradient_error_squared_integral", affineGradientError2); + measure(quadrature, "affine", "reference_gradient_squared_integral", 3. * affineVolume); + measure(quadrature, "affine", "relative_value_l2_error", std::sqrt(affineError2 / affineReference2)); + measure(quadrature, "affine", "relative_gradient_l2_error", std::sqrt(affineGradientError2 / (3. * affineVolume))); + } + measure(0, "reconstruction", "maximum_scatter_roundtrip_error", *std::max_element(scatterErrors.begin(), scatterErrors.end())); + measure(0, "reconstruction", "maximum_surface_trace_error", traceError); + const auto productionCacheAfter = productionMapper.GetGeometryCacheStatistics(); + const bool cacheUnchanged = productionCacheBefore.hits == productionCacheAfter.hits && + productionCacheBefore.misses == productionCacheAfter.misses && productionCacheBefore.points == productionCacheAfter.points; + const bool unchanged = cacheUnchanged && + generation == problem.GetPreparationGeneration() && state.normalizedOperator->IsPrepared() && + d::Identical(physicalBase, state.normalizedOperator->GetPhysicalState()) && + d::Identical(preparedResidual, state.normalizedOperator->GetPhysicalResidual()) && + d::Identical(accepted, state.acceptedNormalizedState) && d::Identical(acceptedResidual, state.acceptedNormalizedResidual) && + d::Identical(acceptedPhysical, state.AcceptedPhysicalState()) && d::Identical(trial, state.trialNormalizedState) && + d::Identical(trialResidual, state.trialNormalizedResidual) && d::Identical(candidatePhysical, state.candidatePhysicalState) && + d::Identical(baseDisplacement, physical.GetGeneratedVolumeDisplacement()) && + d::Identical(stateFactors, state.normalizedOperator->GetNormalization().StateFactors()) && + d::Identical(residualFactors, state.normalizedOperator->GetNormalization().ResidualFactors()); + measure(0, "invariant", "accepted_trial_prepared_normalization_unchanged", unchanged); + measure(0, "invariant", "production_mapping_cache_unchanged", cacheUnchanged); + measure(0, "invariant", "trial_preparations", 0); + measure(0, "invariant", "volume_direction_actions", 1); + d::Require(unchanged, "solver state or normalization changed during physical observation"); +} +} // namespace experiment::plateau_diagnostics diff --git a/experiments/shape_field_weak_diagnostics.hpp b/experiments/shape_field_weak_diagnostics.hpp new file mode 100644 index 0000000..b0e5ab9 --- /dev/null +++ b/experiments/shape_field_weak_diagnostics.hpp @@ -0,0 +1,209 @@ +#pragma once + +// Experiment-only: one frozen Jacobian and fixed normalization for every split. +#include "plateau_nonlinear_step_diagnostics.hpp" +#include "shape_field_physical_diagnostics.hpp" +#include + +namespace experiment::plateau_diagnostics { +template +void ProbeWeakShapeField(State &state, const mean_field::fem::FEM &fem, + const mfem::Vector &normalizedBase, const mfem::Vector *suppliedDirection, + const std::filesystem::path &output, NonlinearStepProbeOptions options = {}) { + using namespace mean_field; + namespace d = nonlinear_step_detail; + 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 || + (suppliedDirection && (suppliedDirection->Size() != size || !d::Finite(*suppliedDirection)))) + throw std::invalid_argument("Weak shape probes require one rank, finite compatible vectors and valid controls."); + const auto begin = d::Clock::now(); + 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); + std::filesystem::create_directories(output); + const auto file = [&](const std::string &name) { return output / ("weak_shape_" + 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 manifest = d::Open(file("layout.csv")); + manifest << "kind,block,offset,size\n"; + const auto values = state.Problem().GetManifest().valueBlocks(); + const auto rows = state.Problem().GetManifest().residualBlocks(); + for (const auto &b : values) manifest << "value," << b.stableId << ',' << b.offset << ',' << b.size << '\n'; + for (const auto &b : rows) manifest << "residual," << b.stableId << ',' << b.offset << ',' << b.size << '\n'; + manifest.flush(); + mfem::Vector residual(size), direction(size), action(size), physicalBase, physicalDirection; + bool haveBase = false, finalized = false; + int preparations = 0, applications = 0; + double maximumRestoreError = 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("Weak shape probe cannot restore its base."); + if (!haveBase) { residual = state.trialNormalizedResidual; haveBase = true; } + if (!d::Finite(state.trialNormalizedResidual)) throw std::runtime_error("Nonfinite frozen residual."); + maximumRestoreError = std::max(maximumRestoreError, d::Difference(residual, state.trialNormalizedResidual)); + }; + const auto apply = [&](const mfem::Vector &v, mfem::Vector &a) { + state.normalizedOperator->Mult(v, a); + ++applications; + if (!d::Finite(a)) throw std::runtime_error("Nonfinite weak shape Jacobian action."); + }; + const auto finish = [&]() { + // Finalization is attempted once, including on exceptional paths; an + // invariant failure must not be obscured by overwriting its artifacts. + 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()); + 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("maximum_restore_error", maximumRestoreError); + d::Save(file("restored_residual.txt"), state.trialNormalizedResidual); + metric("preparations", preparations); + metric("direct_jacobian_applications", applications); + metric("total_probe_seconds", d::Seconds(begin)); + summary.flush(); + if (!acceptedSame || !factorsSame || maximumRestoreError != 0.) + throw std::runtime_error("Weak shape probe changed frozen invariants."); + }; + try { + restore(); + state.physicalInverse->Refresh(); + state.normalizedOperator->DenormalizeState(base, physicalBase); + metric("schema_version", 1); + metric("state_size", size); + metric("supplied_direction", suppliedDirection != nullptr); + metric("accepted_commits", 0); + metric("outer_tolerance", options.linear.relativeTolerance); + metric("outer_absolute_tolerance", options.linear.absoluteTolerance); + metric("maximum_iterations", options.linear.maximumIterations); + metric("restart", options.restartLength); + if (suppliedDirection) direction = *suppliedDirection; + else { + mfem::Vector rhs(residual); + rhs *= -1.; + direction = 0.; + auto backend = solver::linear::prepareLinearBackend( + solver::linear::FGMRES({.restartLength = options.restartLength, .printLevel = -1}), + *state.normalizedOperator, *state.normalizedInverse, state.Problem().GetCommunicator()); + const auto solve = backend.Solve(rhs, direction, options.linear); + metric("linear_status", int(solve.status)); + metric("linear_iterations", solve.iterations); + metric("linear_seconds", solve.solveSeconds); + metric("linear_reported_true_norm", solve.trueResidualNorm); + restore(); + } + if (!d::Finite(direction) || d::Norm(direction) == 0.) + throw std::runtime_error("Weak shape probe requires a nonzero finite direction."); + apply(direction, action); + mfem::Vector defect(action); + defect += residual; + const double threshold = options.linear.ConvergenceThreshold(d::Norm(residual)); + for (const double norm : {d::Norm(residual), d::Norm(direction), d::Norm(action), d::Norm(defect), threshold}) + if (!std::isfinite(norm)) throw std::runtime_error("Overflowing weak shape norm or linear threshold."); + metric("base_norm", d::Norm(residual)); + metric("direction_norm", d::Norm(direction)); + metric("action_norm", d::Norm(action)); + metric("verified_defect_norm", d::Norm(defect)); + metric("verified_relative_defect", d::Ratio(d::Norm(defect), d::Norm(residual))); + metric("verified_threshold", threshold); + metric("verified_gate", d::Norm(defect) <= threshold); + // A finite failed solve is still a directional observation, never an + // accepted Newton step or an equal-accuracy order comparison. + d::Save(file("base.txt"), base); + d::Save(file("physical_base.txt"), physicalBase); + d::Save(file("residual.txt"), residual); + d::Save(file("direction.txt"), direction); + d::Save(file("action.txt"), action); + d::Save(file("defect.txt"), defect); + d::Save(file("state_factors.txt"), stateFactors); + d::Save(file("residual_factors.txt"), residualFactors); + state.normalizedOperator->DenormalizeState(direction, physicalDirection); + d::Save(file("physical_direction.txt"), physicalDirection); + + auto split = d::Open(file("blocks.csv")); + split << "column,value_block,value_offset,value_size,residual_block,residual_offset,residual_size,direction_norm,action_norm\n"; + mfem::Vector sum(size), blockDirection(size), blockAction(size); + sum = 0.; + int column = 0; + for (const auto &b : values) { + blockDirection = 0.; + for (int i = b.offset; i < b.offset + b.size; ++i) blockDirection(i) = direction(i); + apply(blockDirection, blockAction); + sum += blockAction; + d::Save(file("action_b" + std::to_string(column) + ".txt"), blockAction); + for (const auto &r : rows) + split << column << ',' << b.stableId << ',' << b.offset << ',' << b.size << ',' << r.stableId << ',' + << r.offset << ',' << r.size << ',' << d::Norm(blockDirection) << ',' + << d::Norm(blockAction, r.offset, r.size) << '\n'; + ++column; + } + metric("split_columns", column); + metric("split_action_sum_error", d::Difference(sum, action)); + split.flush(); + MeasureShapeFieldPhysical(state, fem, physicalBase, physicalDirection, output); + + auto trials = d::Open(file("trials.csv")); + trials << "trial,alpha,prepared,reason,stage,residual_norm,remainder_norm,remainder_over_alpha,linear_prediction_norm\n"; + int trial = 0; + for (const double alpha : {1./128., 1./512., 1./2048.}) { + mfem::Vector candidate(base), prediction(residual); + candidate.Add(alpha, direction); + prediction.Add(alpha, action); + const auto prepared = prepare(candidate); + d::Save(file("candidate_a" + std::to_string(trial) + ".txt"), candidate); + const double nan = std::numeric_limits::quiet_NaN(); + double residualNorm = nan, remainderNorm = nan; + if (prepared) { + if (!d::Finite(state.trialNormalizedResidual)) throw std::runtime_error("Nonfinite trial residual."); + mfem::Vector remainder(state.trialNormalizedResidual); + remainder -= prediction; + residualNorm = d::Norm(state.trialNormalizedResidual); + remainderNorm = d::Norm(remainder); + d::Save(file("residual_a" + std::to_string(trial) + ".txt"), state.trialNormalizedResidual); + } + trials << trial++ << ',' << alpha << ',' << bool(prepared) << ',' + << (prepared ? -1 : int(prepared.error().reason)) << ',' + << (prepared ? -1 : int(prepared.error().stage)) << ',' << residualNorm << ',' << remainderNorm << ',' + << remainderNorm / alpha << ',' << d::Norm(prediction) << '\n'; + trials.flush(); + } + restore(); + mfem::Vector restoredAction(size); + apply(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("Frozen Jacobian action changed."); + finish(); + metric("probe_exception", 0); + summary.flush(); + } catch (...) { + const auto error = std::current_exception(); + metric("probe_exception", 1); + if (!finalized) finish(); + summary.flush(); + std::rethrow_exception(error); + } +} +} // namespace experiment::plateau_diagnostics