experiment: split frozen operator curvature into shape and field contributions
This commit is contained in:
156
experiments/plateau_shape_curvature_README.md
Normal file
156
experiments/plateau_shape_curvature_README.md
Normal file
@@ -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 `<L,R_full>`, `<L,R_shape>`,
|
||||
`<L,R_fields>`, and `<L,R_mixed>`. These distinguish curvature amplitude from
|
||||
its effect on the merit. The file checks both scalar identities
|
||||
|
||||
```
|
||||
||F_full_trial||² = ||L||² + 2<L,R_full> + ||R_full||²
|
||||
<L,R_full> = <L,R_shape> + <L,R_fields> + <L,R_mixed>.
|
||||
```
|
||||
|
||||
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.
|
||||
420
experiments/plateau_shape_curvature_diagnostics.hpp
Normal file
420
experiments/plateau_shape_curvature_diagnostics.hpp
Normal file
@@ -0,0 +1,420 @@
|
||||
#pragma once
|
||||
|
||||
// Diagnostic-only. Include after `import mean_field;`. No candidate is committed.
|
||||
#include "plateau_enthalpy_step_diagnostics.hpp"
|
||||
#include <array>
|
||||
#include <exception>
|
||||
#include <filesystem>
|
||||
#include <fstream>
|
||||
#include <iomanip>
|
||||
#include <string_view>
|
||||
#include <utility>
|
||||
|
||||
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<long double>(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<long double>(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<double>::quiet_NaN();
|
||||
const double ratio = numerator / denominator;
|
||||
return std::isfinite(ratio) ? ratio : std::numeric_limits<double>::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<double>::quiet_NaN();
|
||||
long double dot = 0.;
|
||||
for (int i = offset; i < offset + size; ++i)
|
||||
dot += static_cast<long double>(a(i)) * b(i);
|
||||
return std::clamp(static_cast<double>(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 <class State>
|
||||
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<double>::quiet_NaN();
|
||||
const double eps = std::numeric_limits<mfem::real_t>::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<mfem::Vector, 3> directions{mfem::Vector(normalizedDirection), mfem::Vector(size),
|
||||
mfem::Vector(normalizedDirection)};
|
||||
constexpr std::array<const char *, 3> 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<mfem::Vector, 3> 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<double, 6> 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<mfem::Vector, 3> remainders{mfem::Vector(size), mfem::Vector(size), mfem::Vector(size)};
|
||||
std::array<bool, 3> admissible{}, resolved{};
|
||||
std::array<mfem::Vector, 3> floors{mfem::Vector(static_cast<int>(rows.size())),
|
||||
mfem::Vector(static_cast<int>(rows.size())), mfem::Vector(static_cast<int>(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<std::pair<const char *, const mfem::Vector *>, 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
|
||||
Reference in New Issue
Block a user