Add frozen physical shape-field compensation diagnostics and controlled order study plan
This commit is contained in:
65
experiments/SHAPE_FIELD_WEAK_RESPONSE_2026-09-18.md
Normal file
65
experiments/SHAPE_FIELD_WEAK_RESPONSE_2026-09-18.md
Normal file
@@ -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.
|
||||
@@ -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;
|
||||
}
|
||||
|
||||
456
experiments/shape_field_physical_diagnostics.hpp
Normal file
456
experiments/shape_field_physical_diagnostics.hpp
Normal file
@@ -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 <algorithm>
|
||||
#include <array>
|
||||
#include <cmath>
|
||||
#include <cstring>
|
||||
#include <filesystem>
|
||||
#include <fstream>
|
||||
#include <iomanip>
|
||||
#include <limits>
|
||||
#include <memory>
|
||||
#include <stdexcept>
|
||||
#include <string>
|
||||
#include <string_view>
|
||||
#include <utility>
|
||||
|
||||
#include <mfem.hpp>
|
||||
#include <mpi.h>
|
||||
|
||||
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<long double>::max_digits10);
|
||||
return output;
|
||||
}
|
||||
|
||||
struct Block { int offset; int size; };
|
||||
|
||||
template <class Blocks>
|
||||
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 <class Field, class Blocks>
|
||||
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, DomainSchema>(
|
||||
*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 <class State>
|
||||
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<double, 6> scatterErrors{
|
||||
d::Scatter<field::Density>(blocks, "density", physicalBase, rho),
|
||||
d::Scatter<field::Density>(blocks, "density", physicalDirection, drho),
|
||||
d::Scatter<field::Enthalpy>(blocks, "specific_enthalpy", physicalBase, h),
|
||||
d::Scatter<field::Enthalpy>(blocks, "specific_enthalpy", physicalDirection, dh),
|
||||
d::Scatter<field::Gravity>(blocks, "gravity_potential", physicalBase, phi),
|
||||
d::Scatter<field::Gravity>(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<utils::domain::Vacuum>()},
|
||||
std::make_unique<d::RejectExteriorMap>());
|
||||
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<double>::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<std::pair<std::string_view, const mfem::ParFiniteElementSpace *>, 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<std::string_view, 3> names{"density", "specific_enthalpy", "gravity_potential"};
|
||||
const std::array<mfem::ParGridFunction *, 3> 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<d::ScalarIntegrals, 3> 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<utils::domain::Stellar>(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<d::ScalarIntegrals, 3> integrals;
|
||||
long double volume = 0., u2 = 0., v2 = 0., referenceVolume = 0., affineE2 = 0., affineR2 = 0., affineGE2 = 0.;
|
||||
double rmin = std::numeric_limits<double>::infinity(), rmax = 0.,
|
||||
physicalRmin = std::numeric_limits<double>::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<long double>(weight) * u(component) * u(component);
|
||||
v2 += static_cast<long double>(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<long double>(affineValue(row)) - x(row);
|
||||
affineE2 += referenceWeight * error * error;
|
||||
affineR2 += static_cast<long double>(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
|
||||
209
experiments/shape_field_weak_diagnostics.hpp
Normal file
209
experiments/shape_field_weak_diagnostics.hpp
Normal file
@@ -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 <vector>
|
||||
|
||||
namespace experiment::plateau_diagnostics {
|
||||
template <class State>
|
||||
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<double>::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
|
||||
Reference in New Issue
Block a user