From c10725beceff7787828102ff5e9d61f03a566820 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Thu, 17 Sep 2026 08:18:09 -0400 Subject: [PATCH] experiment: split frozen operator curvature into shape and field contributions --- experiments/plateau_shape_curvature_README.md | 156 +++++++ .../plateau_shape_curvature_diagnostics.hpp | 420 ++++++++++++++++++ 2 files changed, 576 insertions(+) create mode 100644 experiments/plateau_shape_curvature_README.md create mode 100644 experiments/plateau_shape_curvature_diagnostics.hpp diff --git a/experiments/plateau_shape_curvature_README.md b/experiments/plateau_shape_curvature_README.md new file mode 100644 index 0000000..d1957f9 --- /dev/null +++ b/experiments/plateau_shape_curvature_README.md @@ -0,0 +1,156 @@ +# Shape and field Taylor decomposition + +`plateau_shape_curvature_diagnostics.hpp` is an experiment-only observer, included +after `import mean_field;`. It also includes the existing sampled enthalpy-bound +helper. The root driver supplies an actual frozen normalized base and selected +normalized correction: + +```cpp +experiment::plateau_diagnostics::ProbeShapeCurvature( + state, finiteElements, polytropicIndex, normalizedBase, normalizedDirection, + outputDirectory); +``` + +Use a distinct output directory for each state/direction. The inputs are copied +before state workspaces change, so aliases of those workspaces are safe. The +caller must use this context's normalization and matching EOS index; the helper +does not reconstruct a Newton correction from a damped trial. If the driver does +that reconstruction, its alpha and roundoff uncertainty belong in the enclosing +provenance, particularly when the accepted step is tiny. + +## Experiment and preparation budget + +The actual value manifest must contain exactly one nonempty block whose stable +ID starts with `surface`. Only that block enters p_shape; all remaining value +coordinates, including scalar multipliers and angular velocity, enter p_fields. +No value/residual offsets or field sizes are hard-coded. + +The helper prepares the exact supplied base immediately before each of its +three Jacobian actions. It records repeated base residuals and checks +`p=p_shape+p_fields` and `Jp=Jp_shape+Jp_fields`. These measurements are made +before any step-size gate, even when the direction's enthalpy boundary makes +curvature differencing unusable. + +At that same prepared base, it measures geometry limits for all three directions +using the production preflight rules, and scans the exact production closure and +pressure enthalpy quadrature for full and fields directions separately. Their +enthalpy bounds should agree because the shape direction changes no enthalpy +coefficient. The shape-only bound is therefore infinite for an admissible base. +The common maximum alpha is the minimum of one, all geometry-safe bounds, and +0.8 times the finite full/fields enthalpy crossings. The geometric estimator uses +its existing safety fraction 0.9. + +Six common scales are used: alpha_max times 1, 0.1, 0.01, 0.001, 0.0001, 0.00001. +Each scale prepares full, shape-only, and fields-only trials. Every trial vector +is constructed directly from the copied base; trials never accumulate. The +ordinary budget is 22 preparations (three for J actions, 18 trials, final base +restoration), three Jv, and zero preconditioner applications or linear solves. +An unusable common cap skips all trial preparations, retaining the action/bound +diagnostics. Exceptions attempt one additional base restoration before rethrow. + +For each residual block and direction it records the **vector** Taylor remainder + +``` +R_d(alpha) = F(x + alpha p_d) - F(x) - alpha Jp_d. +``` + +The mixed term is formed by vector subtraction before taking its norm: + +``` +R_mixed = R_full - R_shape - R_fields. +``` + +For each block the mixed file also records the linear prediction +`L=F_base+alpha Jp_full`, its norm, the measured full-trial residual norm, and +long-double-accumulated signed products ``, ``, +``, and ``. These distinguish curvature amplitude from +its effect on the merit. The file checks both scalar identities + +``` +||F_full_trial||² = ||L||² + 2 + ||R_full||² + = + + . +``` + +The errors of these identities are exported, not forced to zero. Individual +signed projections remain available when their own trial is valid, even if +another direction at that alpha is inadmissible; mixed projections require all +three. Summing signed products across residual blocks gives the full normalized +merit decomposition. A large norm alone does not establish an adverse merit +contribution. + +It is not the difference of three norms. Six remainder cosines and the sum of +component norms divided by the full norm show cancellation between terms. All +three trials must be admissible before a mixed result exists; otherwise mixed +values are NaN. Individual valid trials remain available even if another trial +at that alpha fails. + +At fixed shape, the two gravity residual blocks should be affine in the remaining +independent fields. The `fields` rows therefore test an arithmetic-level gravity +remainder. The helper records that measurement without forcing it to zero. +Closure, pressure/mechanics, rotation and scalar constraints need not be affine +in these remaining coordinates; the same expectation does not apply to every +residual block. + +## Arithmetic limits + +The state addition error is measured explicitly as +`(fl(x+alpha p)-x)-alpha p`. A nonzero direction's increment is flagged resolved +only if its expected norm exceeds `64 epsilon (1+||x||)`, its actual increment +is nonzero, and its addition error is below 1% of the expected increment. These +are diagnostic numerical heuristics, not physical tolerances. Tiny components +may be unresolved even when the whole direction is resolved. + +Per-block residual roundoff scales use the larger of eight times the observed +pre-action base-restoration drift and +`128 epsilon (1+||F_base,block||+alpha||Jp_block||)`. Mixed scales sum the three +component scales. The `above_roundoff_scale` flag also requires admissibility +and a resolved state increment. This heuristic is not a certified bound on +matrix-free assembly error: terms may cancel before the residual is formed. +Use the alpha trend, the fields-only gravity measurement, and the final +restoration CSV together. If final restoration drift exceeds the pre-action +drift, increase the interpretation floor accordingly; previously written rows +are not retroactively rescaled. + +If the common alpha cap is zero/nonfinite, its largest full increment is below +the update threshold, or the Jacobian split exceeds its recorded roundoff +heuristic, `common_step_regime_usable=0` and no curvature is inferred. Below-floor +remainders are unresolved, not evidence of zero physical curvature. Inadmissible +trials and zero-denominator cosines/relative errors stay NaN, never fabricated +zeros. Root norms remain in the supplied normalization; manifest scale metadata +is not used as a replacement normalization. + +`alpha_squared_valid` requires a positive finite floating-point alpha². If it +underflows or is nonfinite, scaled curvature is NaN and the resolved-curvature +flag is false. Invalid denominators or overflow in scaled ratios also produce +NaN; a nonpositive/nonfinite alpha itself is skipped and recorded in the summary. +Absolute measured remainders may still be available when squared scaling is not. + +## Outputs and lifecycle + +- `shape_curvature_summary.csv`: base/direction sizes and norms, split checks, + arithmetic thresholds, common cap, preparation count, and accepted-vector + invariance checks. +- `shape_curvature_preparations.csv`: every attempted preparation with exact + available rejection reason/stage, names, meaningful EOS code, and the reported + minimum mapping determinant. The production payload does not identify an + offending element or quadrature point; the helper does not invent one. +- `shape_curvature_restorations.csv`: per-block residual differences on each + repeated base preparation, including final restoration. +- `shape_curvature_directions.csv` and `shape_curvature_actions.csv`: norms + split with the actual value/residual manifests and Jacobian reconstruction. +- `shape_curvature_bounds.csv`: geometry and sampled pressure/closure limits, + including the enthalpy scanner's limiting sample coordinates. +- `shape_curvature_remainders.csv`: absolute R, R/alpha², derivative error + R/alpha, relative derivative error, arithmetic flags, and actual state-addition + errors for every block/direction/alpha. +- `shape_curvature_mixed.csv`: the mixed vector's norm/scaling, cosines, linear + prediction and trial norms, signed merit projections, and identity errors. + +No trial is committed and no accepted vector is changed. The helper leaves the +supplied frozen base prepared and restores the entry normalized correction, +physical correction, and volume-displacement-direction workspaces. Preparation +revisions advance, so a caller that subsequently applies an inverse must refresh +it; the enclosing `WithState` hook normally restores its accepted state and +inverse. The helper performs no preconditioner refresh or policy change itself. +Its rejection/event files are flushed incrementally. No physics, normalization, +production API, driver, or build configuration is changed by this helper. diff --git a/experiments/plateau_shape_curvature_diagnostics.hpp b/experiments/plateau_shape_curvature_diagnostics.hpp new file mode 100644 index 0000000..5549708 --- /dev/null +++ b/experiments/plateau_shape_curvature_diagnostics.hpp @@ -0,0 +1,420 @@ +#pragma once + +// Diagnostic-only. Include after `import mean_field;`. No candidate is committed. +#include "plateau_enthalpy_step_diagnostics.hpp" +#include +#include +#include +#include +#include +#include +#include + +namespace experiment::plateau_diagnostics { +namespace shape_curvature_detail { +inline double Norm(const mfem::Vector &v, int offset, int size) { + long double sum = 0.; + for (int i = offset; i < offset + size; ++i) + sum += static_cast(v(i)) * v(i); + return std::sqrt(sum); +} +inline long double Dot(const mfem::Vector &a, const mfem::Vector &b, int offset, int size) { + long double sum = 0.; + for (int i = offset; i < offset + size; ++i) + sum += static_cast(a(i)) * b(i); + return sum; +} +inline double SafeRatio(double numerator, double denominator) { + if (!std::isfinite(numerator) || !std::isfinite(denominator) || denominator <= 0.) + return std::numeric_limits::quiet_NaN(); + const double ratio = numerator / denominator; + return std::isfinite(ratio) ? ratio : std::numeric_limits::quiet_NaN(); +} +inline double Cosine(const mfem::Vector &a, const mfem::Vector &b, int offset, int size) { + const double na = Norm(a, offset, size), nb = Norm(b, offset, size); + if (!(na > 0.) || !(nb > 0.)) + return std::numeric_limits::quiet_NaN(); + long double dot = 0.; + for (int i = offset; i < offset + size; ++i) + dot += static_cast(a(i)) * b(i); + return std::clamp(static_cast(dot / na / nb), -1., 1.); +} +inline bool Finite(const mfem::Vector &v) { + for (int i = 0; i < v.Size(); ++i) + if (!std::isfinite(v(i))) + return false; + return true; +} +inline std::ofstream Open(const std::filesystem::path &path) { + std::ofstream file(path); + file.exceptions(std::ios::failbit | std::ios::badbit); + file << std::setprecision(17); + return file; +} +inline const char *ReasonName(mean_field::operators::StellarEquilibriumPreparationRejectionReason reason) { + using R = mean_field::operators::StellarEquilibriumPreparationRejectionReason; + switch (reason) { + case R::inverted_geometry: return "inverted_geometry"; + case R::non_finite_geometry: return "non_finite_geometry"; + case R::thermodynamic_domain: return "thermodynamic_domain"; + case R::non_finite_thermodynamics: return "non_finite_thermodynamics"; + case R::inadmissible_physics: return "inadmissible_physics"; + case R::non_finite_physics: return "non_finite_physics"; + } + return "unknown"; +} +inline const char *StageName(mean_field::operators::StellarEquilibriumPreparationStage stage) { + using S = mean_field::operators::StellarEquilibriumPreparationStage; + switch (stage) { + case S::unspecified: return "unspecified"; + case S::generated_geometry: return "generated_geometry"; + case S::gravity: return "gravity"; + case S::barotropic_closure: return "barotropic_closure"; + case S::hydrostatic_equilibrium: return "hydrostatic_equilibrium"; + case S::displacement_residual: return "displacement_residual"; + case S::pressure_force: return "pressure_force"; + case S::gravity_displacement_force: return "gravity_displacement_force"; + case S::rotational_displacement_force: return "rotational_displacement_force"; + case S::displacement_composition: return "displacement_composition"; + case S::mass_normalization: return "mass_normalization"; + case S::model_specification: return "model_specification"; + } + return "unknown"; +} +} // namespace shape_curvature_detail + +// The caller supplies the actual frozen base and selected correction in the +// CURRENT normalization. Fields means every value coordinate except the one +// manifest surface block, including scalar multipliers and angular velocity. +// The usual path costs 22 preparations, three Jv, and no inverse applications. +// On return the frozen base is prepared, caller correction workspaces are +// restored, and the physical inverse may be stale until the caller refreshes it. +template +void ProbeShapeCurvature(State &s, const mean_field::fem::FEM &f, double n, + const mfem::Vector &normalizedBase, const mfem::Vector &normalizedDirection, + const std::filesystem::path &output) { + namespace d = shape_curvature_detail; + using RejectionReason = mean_field::operators::StellarEquilibriumPreparationRejectionReason; + const double nan = std::numeric_limits::quiet_NaN(); + const double eps = std::numeric_limits::epsilon(); + int ranks = 0; + MPI_Comm_size(s.Problem().GetCommunicator(), &ranks); + if (ranks != 1) + throw std::invalid_argument("Shape-curvature diagnostics require one MPI rank."); + const int size = s.Problem().StateSize(); + if (size != s.Problem().EquationSize() || normalizedBase.Size() != size || normalizedDirection.Size() != size || + !d::Finite(normalizedBase) || !d::Finite(normalizedDirection)) + throw std::invalid_argument("Shape-curvature diagnostics require finite compatible root vectors."); + // Copy inputs before touching state: a caller may pass aliases of its workspaces. + const mfem::Vector base(normalizedBase); + std::array directions{mfem::Vector(normalizedDirection), mfem::Vector(size), + mfem::Vector(normalizedDirection)}; + constexpr std::array names{"full", "shape", "fields"}; + directions[1] = 0.; + const auto values = s.Problem().GetManifest().valueBlocks(); + const auto rows = s.Problem().GetManifest().residualBlocks(); + int shapeBlocks = 0, shapeOffset = -1, shapeSize = -1; + for (const auto &block : values) { + if (!std::string_view(block.stableId).starts_with("surface")) + continue; + ++shapeBlocks; + shapeOffset = block.offset; + shapeSize = block.size; + for (int i = block.offset; i < block.offset + block.size; ++i) { + directions[1](i) = directions[0](i); + directions[2](i) = 0.; + } + } + if (shapeBlocks != 1 || shapeOffset < 0 || shapeSize <= 0) + throw std::invalid_argument("Shape-curvature diagnostics require exactly one nonempty surface value block."); + std::filesystem::create_directories(output); + if (std::filesystem::exists(output / "shape_curvature_summary.csv")) + throw std::invalid_argument("Shape-curvature output already exists; use a distinct directory per state."); + auto summary = d::Open(output / "shape_curvature_summary.csv"); + auto events = d::Open(output / "shape_curvature_preparations.csv"); + auto restores = d::Open(output / "shape_curvature_restorations.csv"); + auto directionFile = d::Open(output / "shape_curvature_directions.csv"); + auto actionsFile = d::Open(output / "shape_curvature_actions.csv"); + auto boundsFile = d::Open(output / "shape_curvature_bounds.csv"); + auto remainderFile = d::Open(output / "shape_curvature_remainders.csv"); + auto mixedFile = d::Open(output / "shape_curvature_mixed.csv"); + summary << "metric,value\n"; + events << "preparation,kind,alpha,admissible,rejection_reason,rejection_name,stage,stage_name,eos_code_applicable," + "eos_code,reported_minimum_jacobian,residual_norm\n"; + restores << "preparation,kind,block,base_residual_norm,restoration_difference_norm\n"; + directionFile << "direction,value_block,offset,size,normalized_norm\n"; + actionsFile << "block,base_residual_norm,full_action_norm,shape_action_norm,fields_action_norm," + "split_error_norm,relative_split_error,shape_fields_cosine\n"; + boundsFile << "direction,metric,value\n"; + remainderFile << "alpha,direction,block,admissible,state_increment_resolved,alpha_squared_valid,base_residual_norm,trial_residual_norm," + "action_norm,remainder_norm,remainder_over_alpha_squared,difference_quotient_error," + "relative_difference_quotient_error,heuristic_roundoff_scale,remainder_over_roundoff_scale," + "above_roundoff_scale,expected_step_norm,actual_step_norm,state_addition_error_norm\n"; + mixedFile << "alpha,block,all_trials_admissible,all_state_increments_resolved,alpha_squared_valid,full_remainder_norm," + "shape_remainder_norm,fields_remainder_norm,mixed_remainder_norm,mixed_over_alpha_squared," + "heuristic_roundoff_scale,mixed_over_roundoff_scale,above_roundoff_scale," + "full_shape_cosine,full_fields_cosine,full_mixed_cosine,shape_fields_cosine," + "shape_mixed_cosine,fields_mixed_cosine,component_norm_sum_over_full,linear_prediction_norm," + "full_trial_norm,linear_dot_full,linear_dot_shape,linear_dot_fields,linear_dot_mixed," + "merit_identity_error,projection_decomposition_error\n"; + summary << "machine_epsilon," << eps << "\nshape_offset," << shapeOffset << "\nshape_size," << shapeSize + << "\nbase_state_norm," << base.Norml2() << "\nfull_direction_norm," << directions[0].Norml2() << '\n'; + mfem::Vector directionSplitError(directions[0]); + directionSplitError -= directions[1]; + directionSplitError -= directions[2]; + summary << "direction_split_error_norm," << directionSplitError.Norml2() << '\n'; + for (int k = 0; k < 3; ++k) + for (const auto &block : values) + directionFile << names[k] << ',' << block.stableId << ',' << block.offset << ',' << block.size << ',' + << d::Norm(directions[k], block.offset, block.size) << '\n'; + directionFile.flush(); + + const mfem::Vector savedCorrection(s.normalizedCorrection), savedPhysicalCorrection(s.physicalCorrection), + savedVolumeDirection(s.volumeDisplacementDirection), savedAcceptedState(s.acceptedNormalizedState), + savedAcceptedResidual(s.acceptedNormalizedResidual), savedAcceptedPhysical(s.AcceptedPhysicalState()); + auto restoreWorkspaces = [&]() { + s.normalizedCorrection = savedCorrection; + s.physicalCorrection = savedPhysicalCorrection; + s.volumeDisplacementDirection = savedVolumeDirection; + }; + int preparationCount = 0; + mfem::Vector baseResidual(size), restorationEnvelope(size); + restorationEnvelope = 0.; + bool haveBaseResidual = false; + auto prepare = [&](const mfem::Vector &state, const char *kind, double alpha) { + s.trialNormalizedState = state; + ++preparationCount; + auto result = s.PrepareTrial(); + events << preparationCount << ',' << kind << ',' << alpha << ',' << bool(result) << ','; + if (result) { + events << "-1,none,-1,none,0,-1," << result->minimumJacobianDeterminant << ',' + << s.trialNormalizedResidual.Norml2() << '\n'; + } else { + const auto &error = result.error(); + const bool eos = error.reason == RejectionReason::thermodynamic_domain || + error.reason == RejectionReason::non_finite_thermodynamics; + events << int(error.reason) << ',' << d::ReasonName(error.reason) << ',' << int(error.stage) << ',' + << d::StageName(error.stage) << ',' << eos << ',' << (eos ? int(error.thermodynamicErrorCode) : -1) + << ',' << error.minimumJacobianDeterminant << ',' << nan << '\n'; + } + events.flush(); + return result; + }; + auto restoreBase = [&](const char *kind) { + const auto result = prepare(base, kind, 0.); + if (!result || !d::Finite(s.trialNormalizedResidual)) + throw std::runtime_error("Shape-curvature diagnostic could not restore its finite frozen base."); + if (!haveBaseResidual) { + baseResidual = s.trialNormalizedResidual; + haveBaseResidual = true; + } + mfem::Vector difference(s.trialNormalizedResidual); + difference -= baseResidual; + for (int i = 0; i < size; ++i) + restorationEnvelope(i) = std::max(restorationEnvelope(i), std::abs(difference(i))); + for (const auto &block : rows) + restores << preparationCount << ',' << kind << ',' << block.stableId << ',' + << d::Norm(baseResidual, block.offset, block.size) << ',' + << d::Norm(difference, block.offset, block.size) << '\n'; + restores.flush(); + }; + + try { + std::array actions{mfem::Vector(size), mfem::Vector(size), mfem::Vector(size)}; + for (int k = 0; k < 3; ++k) { + restoreBase(names[k]); // Every J action is preceded by an exact base preparation. + s.normalizedOperator->Mult(directions[k], actions[k]); + if (!d::Finite(actions[k])) + throw std::runtime_error("Nonfinite shape-curvature Jacobian action."); + } + mfem::Vector splitError(actions[0]); + splitError -= actions[1]; + splitError -= actions[2]; + for (const auto &block : rows) { + const int o = block.offset, count = block.size; + const double full = d::Norm(actions[0], o, count), shape = d::Norm(actions[1], o, count), + fields = d::Norm(actions[2], o, count), error = d::Norm(splitError, o, count); + const double scale = full + shape + fields; + actionsFile << block.stableId << ',' << d::Norm(baseResidual, o, count) << ',' << full << ',' << shape + << ',' << fields << ',' << error << ',' << (scale > 0. ? error / scale : nan) << ',' + << d::Cosine(actions[1], actions[2], o, count) << '\n'; + } + actionsFile.flush(); + const double actionLinearityScale = + 128. * eps * (1. + actions[0].Norml2() + actions[1].Norml2() + actions[2].Norml2()); + const bool linearityResolved = splitError.Norml2() <= actionLinearityScale; + summary << "base_residual_norm," << baseResidual.Norml2() << "\naction_split_error_norm," + << splitError.Norml2() << "\naction_linearity_heuristic_scale," << actionLinearityScale + << "\naction_linearity_within_heuristic," << linearityResolved << '\n'; + + mfem::Vector physicalBase; + s.normalizedOperator->DenormalizeState(base, physicalBase); + double commonCap = 1.; + for (int k = 0; k < 3; ++k) { + s.normalizedCorrection = directions[k]; + const auto geometry = s.EstimateLargestSafeStepSize(1., .9); + boundsFile << names[k] << ",geometry_safe_step," << geometry.stepSize << '\n' + << names[k] << ",geometry_boundary_step," << geometry.boundaryStepSize << '\n' + << names[k] << ",geometry_limited," << geometry.limitedByGeometry << '\n' + << names[k] << ",geometry_limiting_element," << geometry.limitingElement << '\n' + << names[k] << ",geometry_limiting_rule," << geometry.limitingRule << '\n' + << names[k] << ",geometry_limiting_point," << geometry.limitingQuadraturePoint << '\n'; + if (!std::isfinite(geometry.stepSize) || geometry.stepSize <= 0.) + commonCap = nan; + else if (std::isfinite(commonCap)) + commonCap = std::min(commonCap, geometry.stepSize); + // Shape-only has exactly zero enthalpy direction. Full and fields + // are scanned separately and should have identical sampled bounds. + if (k != 1) { + const auto enthalpy = MeasureEnthalpyStep(s, f, n, &physicalBase); + for (const auto &[key, value] : enthalpy) + boundsFile << names[k] << ',' << key << ',' << value << '\n'; + const double bound = enthalpy.at("alpha_thermodynamic"); + if (std::isnan(bound) || bound <= 0.) + commonCap = nan; + else if (std::isfinite(bound) && std::isfinite(commonCap)) + commonCap = std::min(commonCap, .8 * bound); + } + } + boundsFile.flush(); + summary << "common_alpha_cap," << commonCap << '\n'; + // This is an explicitly heuristic representation threshold, not a + // certificate of residual assembly accuracy or a physical tolerance. + const double updateFloor = 64. * eps * (1. + base.Norml2()); + summary << "heuristic_state_update_floor," << updateFloor << '\n'; + const bool usableCap = linearityResolved && std::isfinite(commonCap) && commonCap > 0. && + commonCap * directions[0].Norml2() > updateFloor; + summary << "common_step_regime_usable," << usableCap << '\n'; + summary.flush(); + if (usableCap) { + constexpr std::array fractions{1., .1, .01, .001, .0001, .00001}; + for (double fraction : fractions) { + const double alpha = commonCap * fraction; + if (!std::isfinite(alpha) || alpha <= 0.) { + summary << "skipped_nonpositive_or_nonfinite_alpha," << alpha << '\n'; + continue; + } + const double alphaSquared = alpha * alpha; + const bool alphaSquaredValid = std::isfinite(alphaSquared) && alphaSquared > 0.; + mfem::Vector linearPrediction(baseResidual), fullTrialResidual(size); + linearPrediction.Add(alpha, actions[0]); + fullTrialResidual = nan; + std::array remainders{mfem::Vector(size), mfem::Vector(size), mfem::Vector(size)}; + std::array admissible{}, resolved{}; + std::array floors{mfem::Vector(static_cast(rows.size())), + mfem::Vector(static_cast(rows.size())), mfem::Vector(static_cast(rows.size()))}; + for (int k = 0; k < 3; ++k) { + mfem::Vector trial(base); + trial.Add(alpha, directions[k]); + mfem::Vector step(trial); + step -= base; + mfem::Vector additionError(step); + additionError.Add(-alpha, directions[k]); + const double expectedStep = alpha * directions[k].Norml2(); + const double actualStep = step.Norml2(), additionErrorNorm = additionError.Norml2(); + resolved[k] = directions[k].Norml2() == 0. || + (expectedStep > updateFloor && actualStep > 0. && additionErrorNorm < .01 * expectedStep); + const auto result = prepare(trial, names[k], alpha); + admissible[k] = bool(result) && d::Finite(s.trialNormalizedResidual); + if (admissible[k]) { + if (k == 0) + fullTrialResidual = s.trialNormalizedResidual; + remainders[k] = s.trialNormalizedResidual; + remainders[k] -= baseResidual; + remainders[k].Add(-alpha, actions[k]); + } else { + remainders[k] = nan; + } + int rowIndex = 0; + for (const auto &block : rows) { + const int o = block.offset, count = block.size; + const double fn = d::Norm(baseResidual, o, count), jn = d::Norm(actions[k], o, count); + const double floor = std::max(8. * d::Norm(restorationEnvelope, o, count), + 128. * eps * (1. + fn + alpha * jn)); + floors[k](rowIndex++) = floor; + const double rn = admissible[k] ? d::Norm(remainders[k], o, count) : nan; + const double trialNorm = admissible[k] ? d::Norm(s.trialNormalizedResidual, o, count) : nan; + remainderFile << alpha << ',' << names[k] << ',' << block.stableId << ',' << admissible[k] + << ',' << resolved[k] << ',' << alphaSquaredValid << ',' << fn << ',' << trialNorm << ',' << jn << ',' << rn + << ',' << d::SafeRatio(rn, alphaSquared) << ',' << d::SafeRatio(rn, alpha) << ',' + << d::SafeRatio(d::SafeRatio(rn, alpha), jn) << ',' << floor << ',' << d::SafeRatio(rn, floor) + << ',' << (admissible[k] && resolved[k] && alphaSquaredValid && rn > floor) << ',' << expectedStep + << ',' << actualStep << ',' << additionErrorNorm << '\n'; + } + remainderFile.flush(); + } + const bool allAdmissible = admissible[0] && admissible[1] && admissible[2]; + const bool allResolved = resolved[0] && resolved[1] && resolved[2]; + mfem::Vector mixed(size); + mixed = nan; + if (allAdmissible) { + mixed = remainders[0]; + mixed -= remainders[1]; + mixed -= remainders[2]; // Subtract vectors before norms, preserving cancellation. + } + int rowIndex = 0; + for (const auto &block : rows) { + const int o = block.offset, count = block.size; + const double rf = allAdmissible ? d::Norm(remainders[0], o, count) : nan, + rs = allAdmissible ? d::Norm(remainders[1], o, count) : nan, + rp = allAdmissible ? d::Norm(remainders[2], o, count) : nan, + rm = allAdmissible ? d::Norm(mixed, o, count) : nan; + const double floor = floors[0](rowIndex) + floors[1](rowIndex) + floors[2](rowIndex); + ++rowIndex; + mixedFile << alpha << ',' << block.stableId << ',' << allAdmissible << ',' << allResolved << ',' + << alphaSquaredValid << ',' << rf << ',' << rs << ',' << rp << ',' << rm << ',' + << d::SafeRatio(rm, alphaSquared) << ',' << floor << ',' << d::SafeRatio(rm, floor) << ',' + << (allAdmissible && allResolved && alphaSquaredValid && rm > floor); + mixedFile << ',' << (allAdmissible ? d::Cosine(remainders[0], remainders[1], o, count) : nan) + << ',' << (allAdmissible ? d::Cosine(remainders[0], remainders[2], o, count) : nan) + << ',' << (allAdmissible ? d::Cosine(remainders[0], mixed, o, count) : nan) + << ',' << (allAdmissible ? d::Cosine(remainders[1], remainders[2], o, count) : nan) + << ',' << (allAdmissible ? d::Cosine(remainders[1], mixed, o, count) : nan) + << ',' << (allAdmissible ? d::Cosine(remainders[2], mixed, o, count) : nan) + << ',' << d::SafeRatio(rs + rp + rm, rf); + const long double lf = admissible[0] ? d::Dot(linearPrediction, remainders[0], o, count) : nan, + ls = admissible[1] ? d::Dot(linearPrediction, remainders[1], o, count) : nan, + lp = admissible[2] ? d::Dot(linearPrediction, remainders[2], o, count) : nan, + lm = allAdmissible ? d::Dot(linearPrediction, mixed, o, count) : nan; + const long double identityError = admissible[0] ? + d::Dot(fullTrialResidual, fullTrialResidual, o, count) - + d::Dot(linearPrediction, linearPrediction, o, count) - 2. * lf - + d::Dot(remainders[0], remainders[0], o, count) : nan; + mixedFile << ',' << d::Norm(linearPrediction, o, count) << ',' + << (admissible[0] ? d::Norm(fullTrialResidual, o, count) : nan) << ',' + << lf << ',' << ls << ',' << lp << ',' << lm << ',' << identityError << ',' + << (allAdmissible ? lf - ls - lp - lm : nan) << '\n'; + } + mixedFile.flush(); + } + } + restoreBase("final_restore"); + summary << "preparations," << preparationCount << "\nrestoration_envelope_norm," + << restorationEnvelope.Norml2() << '\n'; + for (const auto &entry : std::array, 3>{{ + {"accepted_normalized_state_change", &s.acceptedNormalizedState}, + {"accepted_residual_change", &s.acceptedNormalizedResidual}, + {"accepted_physical_state_change", &s.AcceptedPhysicalState()}}}) { + const mfem::Vector *reference = entry.second == &s.acceptedNormalizedState ? &savedAcceptedState : + (entry.second == &s.acceptedNormalizedResidual ? &savedAcceptedResidual : &savedAcceptedPhysical); + mfem::Vector difference(*entry.second); + difference -= *reference; + summary << entry.first << ',' << difference.Norml2() << '\n'; + if (difference.Norml2() != 0.) + throw std::runtime_error("Shape-curvature diagnostics detected a changed accepted vector."); + } + summary.flush(); + restoreWorkspaces(); + } catch (...) { + const auto original = std::current_exception(); + try { + restoreBase("exception_restore"); + } catch (...) { + restoreWorkspaces(); + throw; + } + restoreWorkspaces(); + std::rethrow_exception(original); + } +} +} // namespace experiment::plateau_diagnostics