tboudreaux/MeanField
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
7a07caeee33e2addb897b89682678b4cc41b1f3e
MeanField/experiments/system_isolation_profile_harness.cpp

1683 lines
109 KiB
C++
Raw Blame History

#include <algorithm>
#include <chrono>
#include <cmath>
#include <filesystem>
#include <fstream>
#include <functional>
#include <iomanip>
#include <iostream>
#include <limits>
#include <memory>
#include <stdexcept>
#include <string>
#include <cstdint>
#include <utility>
#include <vector>
#include <mfem.hpp>
import mean_field;
namespace {
using namespace mean_field;
using Clock = std::chrono::steady_clock;
using Apply = std::function<void(int, mfem::Vector &)>;
struct Options {
std::string mesh{"sandbox.smesh"};
std::filesystem::path output{"system_isolation_results"};
std::string system{"hdiv-displacement"};
std::string implementation{"both"};
std::string check{"full"};
int refinements{0}, samples{5}, iterations{3}, warmup{1}, directions{3};
int prepareSamples{3}, prepareIterations{1};
double deformation{0.7}, finiteDifferenceStep{1.0e-5}, polytropicIndex{1.0};
bool help{false};
};
int Integer(const std::string &value) {
std::size_t consumed = 0;
const int result = std::stoi(value, &consumed);
if (consumed != value.size()) throw std::invalid_argument("Invalid integer: " + value);
return result;
}
double Real(const std::string &value) {
std::size_t consumed = 0;
const double result = std::stod(value, &consumed);
if (consumed != value.size() || !std::isfinite(result))
throw std::invalid_argument("Invalid finite number: " + value);
return result;
}
bool IsLinearDerivativeSystem(const std::string &name);
Options Parse(int argc, char **argv) {
Options options;
for (int index = 1; index < argc; ++index) {
const std::string argument = argv[index];
auto value = [&]() -> std::string {
if (++index == argc) throw std::invalid_argument("Missing value after " + argument);
return argv[index];
};
if (argument == "--mesh") options.mesh = value();
else if (argument == "--output") options.output = value();
else if (argument == "--system") options.system = value();
else if (argument == "--implementation") options.implementation = value();
else if (argument == "--check") options.check = value();
else if (argument == "--refinements") options.refinements = Integer(value());
else if (argument == "--samples") options.samples = Integer(value());
else if (argument == "--iterations") options.iterations = Integer(value());
else if (argument == "--warmup") options.warmup = Integer(value());
else if (argument == "--directions") options.directions = Integer(value());
else if (argument == "--prepare-samples") options.prepareSamples = Integer(value());
else if (argument == "--prepare-iterations") options.prepareIterations = Integer(value());
else if (argument == "--deformation") options.deformation = Real(value());
else if (argument == "--fd-step") options.finiteDifferenceStep = Real(value());
else if (argument == "--polytropic-index") options.polytropicIndex = Real(value());
else if (argument == "--help") options.help = true;
else throw std::invalid_argument("Unknown argument: " + argument);
}
if (options.system != "hdiv-displacement" && options.system != "hdiv-mass" &&
options.system != "mass-normalization" && options.system != "angular-momentum" &&
options.system != "hydrostatic-displacement" && options.system != "stellar-derivatives" &&
options.system != "all" && options.system != "runtime-optimizations" &&
options.system != "geometry-preparation" && options.system != "normalization-diagonals" &&
options.system != "geometry-inspection" && options.system != "gravity-source-preparation" &&
options.system != "final-optimizations" && options.system != "gravity-source-displacement" &&
options.system != "hydrostatic-preparation" && options.system != "gravity-schur-preparation" &&
!IsLinearDerivativeSystem(options.system))
throw std::invalid_argument("Unknown system; use --help for the available systems.");
if (options.implementation != "both" && options.implementation != "optimized" && options.implementation != "reference")
throw std::invalid_argument("Implementation must be both, optimized, or reference.");
if (options.check != "full" && options.check != "reference" && options.check != "none")
throw std::invalid_argument("Check must be full, reference, or none.");
if (options.refinements < 0 || options.samples < 1 || options.iterations < 1 || options.warmup < 0 ||
options.directions < 1 || options.prepareSamples < 0 || options.prepareIterations < 1 ||
options.finiteDifferenceStep <= 0.0 || options.polytropicIndex <= 0.0)
throw std::invalid_argument("Invalid sample count, refinement, or finite-difference step.");
return options;
}
std::ofstream File(const std::filesystem::path &path) {
std::ofstream stream(path);
stream.exceptions(std::ios::badbit | std::ios::failbit);
stream << std::setprecision(17);
return stream;
}
mfem::Vector Displacement(const fem::FEM &f, double scale, double phase = 0.0) {
mfem::ParGridFunction field(f.displacementFes.get());
mfem::VectorFunctionCoefficient coefficient(3, [=](const mfem::Vector &x, mfem::Vector &u) {
u.SetSize(3);
u(0) = scale * (0.04 * x(0) + 0.01 * x(1) * x(2) + phase * (0.012 + 0.007 * x(2)));
u(1) = scale * (-0.03 * x(1) + 0.008 * x(0) * x(2) + phase * (0.009 * x(0) - 0.003));
u(2) = scale * (0.02 * x(2) - 0.006 * x(0) * x(1) - phase * 0.011 * x(1));
});
field.ProjectCoefficient(coefficient);
mfem::Vector result;
field.GetTrueDofs(result);
return result;
}
mfem::Vector Flux(int size, double phase) {
mfem::Vector result(size);
for (int i = 0; i < size; ++i)
result(i) = std::sin(0.37 * (i + 1) + phase) + 0.31 * std::cos(0.19 * (i + 1) - 0.5 * phase);
return result;
}
double Norm(const mfem::Vector &vector) {
for (int i = 0; i < vector.Size(); ++i)
if (!std::isfinite(vector(i))) throw std::runtime_error("An operator returned non-finite output.");
const double result = vector.Norml2();
if (!std::isfinite(result)) throw std::runtime_error("An output norm is non-finite.");
return result;
}
double RelativeError(const mfem::Vector &actual, const mfem::Vector &expected) {
if (actual.Size() != expected.Size()) throw std::runtime_error("Comparison vector sizes differ.");
mfem::Vector difference(actual);
difference -= expected;
return Norm(difference) / std::max({Norm(actual), Norm(expected), 1.0e-30});
}
double Median(std::vector<double> values) {
std::sort(values.begin(), values.end());
const std::size_t middle = values.size() / 2;
return values.size() % 2 ? values[middle] : 0.5 * (values[middle - 1] + values[middle]);
}
struct CheckReport {
std::ofstream stream;
bool passed{true};
void Add(const std::string &name, double error, double tolerance) {
const bool okay = std::isfinite(error) && error <= tolerance;
stream << name << ',' << error << ',' << tolerance << ',' << okay << '\n';
stream.flush();
passed = passed && okay;
std::cout << "Check " << name << ": " << error << (okay ? " PASS" : " FAIL") << std::endl;
}
};
// All systems use this timing driver. Construction, checks, warmup and output
// reductions are excluded. Alternate order to reduce systematic A/B bias.
void Benchmark(const Options &options, const std::string &system,
const std::string &referenceName, const Apply &optimized, const Apply &reference,
std::ofstream &timings, std::ofstream &summary, double extraPrepareSeconds = 0.0) {
const bool runOptimized = options.implementation != "reference";
const bool runReference = options.implementation != "optimized";
mfem::Vector output;
for (int iteration = 0; iteration < options.warmup; ++iteration) {
if (runOptimized) { optimized(iteration % options.directions, output); Norm(output); }
if (runReference) { reference(iteration % options.directions, output); Norm(output); }
}
std::vector<double> optimizedTimes, referenceTimes;
auto sample = [&](const std::string &name, const Apply &apply, std::vector<double> &times, int index) {
const auto start = Clock::now();
for (int iteration = 0; iteration < options.iterations; ++iteration)
apply(iteration % options.directions, output);
const double seconds = std::chrono::duration<double>(Clock::now() - start).count();
const double perCall = seconds / options.iterations;
times.push_back(perCall);
timings << system << ',' << name << ',' << index << ',' << options.iterations << ','
<< seconds << ',' << perCall << ',' << Norm(output) << '\n';
timings.flush();
std::cout << system << '/' << name << " sample " << index + 1 << ": " << perCall << " s/apply" << std::endl;
};
for (int index = 0; index < options.samples; ++index) {
if (index % 2 == 0 && runReference) sample(referenceName, reference, referenceTimes, index);
if (runOptimized) sample("optimized", optimized, optimizedTimes, index);
if (index % 2 != 0 && runReference) sample(referenceName, reference, referenceTimes, index);
}
if (runOptimized) summary << system << ".optimized_median_seconds=" << Median(optimizedTimes) << '\n';
if (runReference) summary << system << '.' << referenceName << "_median_seconds=" << Median(referenceTimes) << '\n';
if (runOptimized && runReference) {
const double speedup = Median(referenceTimes) / Median(optimizedTimes);
summary << system << ".reference_over_optimized=" << speedup << '\n';
if (extraPrepareSeconds > 0.0) {
const double saving = Median(referenceTimes) - Median(optimizedTimes);
summary << system << ".row_prepare_seconds=" << extraPrepareSeconds << '\n';
summary << system << ".row_prepare_break_even_applications="
<< (saving > 0.0 ? std::ceil(extraPrepareSeconds / saving) : -1.0) << '\n';
for (const int applications : {1, 10, 25})
summary << system << ".reference_over_optimized_including_row_prepare_" << applications << '='
<< applications * Median(referenceTimes) /
(extraPrepareSeconds + applications * Median(optimizedTimes)) << '\n';
}
std::cout << system << " reference/optimized median ratio: " << speedup << std::endl;
}
summary.flush();
}
// State/context updates happen before either timer. Each timed operation
// must rebuild once; changing revisions prevents accidental no-op timings.
void BenchmarkPreparation(
const Options &options, const std::string &system, const std::function<void()> &advance,
const std::function<void()> &optimized, const std::function<void()> &reference,
const std::function<void(mfem::Vector &)> &optimizedResidual,
const std::function<void(mfem::Vector &)> &referenceResidual,
std::ofstream &timings, std::ofstream &summary
) {
if (options.prepareSamples == 0) return;
const bool runOptimized = options.implementation != "reference";
const bool runReference = options.implementation != "optimized";
std::vector<double> optimizedTimes, referenceTimes;
mfem::Vector residual;
for (int sample = 0; sample < options.prepareSamples; ++sample) {
double fastSeconds = 0.0, referenceSeconds = 0.0;
auto measure = [&](const std::function<void()> &operation, double &seconds) {
const auto start = Clock::now();
operation();
seconds += std::chrono::duration<double>(Clock::now() - start).count();
};
for (int iteration = 0; iteration < options.prepareIterations; ++iteration) {
advance();
if (sample % 2 == 0 && runReference) measure(reference, referenceSeconds);
if (runOptimized) measure(optimized, fastSeconds);
if (sample % 2 != 0 && runReference) measure(reference, referenceSeconds);
}
auto record = [&](const std::string &name, double seconds, auto &values, const auto &getResidual) {
getResidual(residual); // Residual copy/norm is outside preparation timing.
const double perCall = seconds / options.prepareIterations;
values.push_back(perCall);
timings << system << ',' << name << ',' << sample << ',' << options.prepareIterations << ','
<< seconds << ',' << perCall << ',' << Norm(residual) << '\n';
std::cout << system << '/' << name << " sample " << sample + 1 << ": "
<< perCall << " s/prepare" << std::endl;
};
if (runOptimized) record("optimized", fastSeconds, optimizedTimes, optimizedResidual);
if (runReference) record("reference", referenceSeconds, referenceTimes, referenceResidual);
timings.flush();
}
if (runOptimized) summary << system << ".optimized_median_seconds=" << Median(optimizedTimes) << '\n';
if (runReference) summary << system << ".reference_median_seconds=" << Median(referenceTimes) << '\n';
if (runOptimized && runReference)
summary << system << ".reference_over_optimized=" << Median(referenceTimes) / Median(optimizedTimes) << '\n';
summary.flush();
}
mfem::Vector ScalarField(mfem::ParFiniteElementSpace &space, double offset, double phase) {
mfem::ParGridFunction field(&space);
mfem::FunctionCoefficient coefficient([=](const mfem::Vector &x) {
return offset + 0.09 * std::sin(0.71 * x(0) + phase) +
0.05 * std::cos(0.63 * x(1) - phase) + 0.02 * x(2);
});
field.ProjectCoefficient(coefficient);
mfem::Vector result;
field.GetTrueDofs(result);
return result;
}
// Mixed relative/absolute error also handles zero scalar columns without
// amplifying harmless cancellation into an artificial relative failure.
double ScaledError(const mfem::Vector &actual, const mfem::Vector &expected) {
if (actual.Size() != expected.Size()) throw std::runtime_error("Comparison vector sizes differ.");
mfem::Vector difference(actual);
difference -= expected;
return Norm(difference) / std::max({1.0, Norm(actual), Norm(expected)});
}
int RunStellarDerivatives(const Options &options, fem::FEM &f, std::ofstream &metadata,
std::ofstream &summary, std::ofstream &timings, CheckReport &checks) {
const bool group = options.system == "all" || options.system == "stellar-derivatives";
const bool runMass = group || options.system == "mass-normalization";
const bool runAngular = group || options.system == "angular-momentum";
const bool runHydro = group || options.system == "hydrostatic-displacement";
using GravityContext = operators::context::gravity_field::GravityFieldLinearizationContext;
std::unique_ptr<GravityContext> gravity;
std::unique_ptr<operators::PreparedMassNormalizationOperator> mass;
std::unique_ptr<operators::PreparedAngularMomentumOperator> angular;
std::unique_ptr<operators::PreparedHydrostaticEquilibriumOperator> hydro;
const auto constructStart = Clock::now();
if (runMass || runAngular) gravity = std::make_unique<GravityContext>(f, *f.domainMapperStateless);
if (runMass) mass = std::make_unique<operators::PreparedMassNormalizationOperator>(f, *f.domainMapperStateless, *gravity);
if (runAngular) angular = std::make_unique<operators::PreparedAngularMomentumOperator>(
f, *f.domainMapperStateless, *gravity,
models::compileConstraint(integral::FixedAngularMomentum(
{.Jtotal = dimensions::AngularMomentumValue{0.41}, .axis = {1.0, -2.0, 3.0}, .center = {0.12, -0.07, 0.03}})));
if (runHydro) hydro = std::make_unique<operators::PreparedHydrostaticEquilibriumOperator>(f, *f.domainMapperStateless);
summary << "stellar.operator_construction_seconds=" << std::chrono::duration<double>(Clock::now() - constructStart).count() << '\n';
const auto baseDensity = ScalarField(*f.densityFes, 0.94, 0.31);
const auto baseDisplacement = Displacement(f, options.deformation);
const auto enthalpy = ScalarField(*f.enthalpyFes, 0.73, 0.27);
const auto potential = ScalarField(*f.gravityPotentialFes, -0.38, 0.61);
mfem::Vector flux(f.gravityFluxFes->GetTrueVSize());
flux = 0.0;
operators::MassNormalizationDependencies md{
.discretization = {1, 1}, .density = {2, 1}, .displacement = {3, 1}, .targetMass = {4, 1}};
operators::AngularMomentumDependencies ad{
.discretization = {1, 1}, .density = {2, 1}, .displacement = {3, 1}, .rotation = {5, 1}};
operators::context::hydrostatic::HydrostaticEquilibriumDependencies hd{
.discretization = {1, 1}, .enthalpy = {6, 1}, .gravityPotential = {7, 1},
.displacement = {3, 1}, .rotation = {5, 1}, .bernoulliConstant = {8, 1}};
double targetMass = 1.23;
auto prepare = [&](const mfem::Vector &density, const mfem::Vector &displacement, double omega,
bool geometryChanged, bool densityChanged, bool rotationChanged,
const std::string &timingLabel) {
if (geometryChanged) ++md.displacement.revision;
if (densityChanged) ++md.density.revision;
if (rotationChanged) ++ad.rotation.revision;
ad.density.revision = md.density.revision;
ad.displacement.revision = md.displacement.revision;
hd.displacement.revision = md.displacement.revision;
hd.rotation.revision = ad.rotation.revision;
auto time = [&](const std::string &name, auto &&operation) {
const auto start = Clock::now();
operation();
if (!timingLabel.empty()) summary << name << '.' << timingLabel << "_seconds="
<< std::chrono::duration<double>(Clock::now() - start).count() << '\n';
};
if (gravity) time("stellar.shared_gravity_prepare", [&] {
gravity->Prepare(
{.density = gravity->GetDensityMap().gather(density),
.displacement = gravity->GetDisplacementMap().gather(displacement),
.gravity_gradient = gravity->GetGravityGradientMap().gather(flux),
.gravity_potential = gravity->GetGravityPotentialMap().gather(potential)},
{.discretization = {1}, .displacement = {md.displacement.revision},
.density = {md.density.revision}, .gravity_gradient = {1}, .gravity_potential = {1}});
});
if (mass) time("mass-normalization.prepare", [&] { mass->Prepare({.targetMass = targetMass}, md); });
if (angular) time("angular-momentum.prepare", [&] { angular->Prepare(omega, ad); });
if (hydro) time("hydrostatic-displacement.prepare", [&] {
mfem::Vector velocity(3), center(3);
velocity(0) = omega / std::sqrt(14.0); velocity(1) = -2.0 * velocity(0); velocity(2) = 3.0 * velocity(0);
center(0) = 0.12; center(1) = -0.07; center(2) = 0.03;
hydro->Prepare({.enthalpy = hydro->GetEnthalpyMap().gather(enthalpy),
.gravityPotential = hydro->GetGravityPotentialMap().gather(potential),
.displacement = hydro->GetDisplacementMap().gather(displacement),
.bernoulliConstant = 0.39}, hd, physics::RigidRotation(velocity, center));
});
};
prepare(baseDensity, baseDisplacement, 0.63, true, true, true, "initial");
std::vector<mfem::Vector> densityDirections, fullDensityDirections, displacementDirections,
hydroDirections, fullDisplacementDirections;
std::vector<double> omegaDirections;
for (int i = 0; i < options.directions; ++i) {
fullDensityDirections.push_back(ScalarField(*f.densityFes, 0.17, 0.4 + i));
fullDisplacementDirections.push_back(Displacement(f, 0.43, 0.2 + 0.7 * i));
// Include a direction with variation at the element scale as well
// as smooth shear/translation directions. It is not a base state.
if (i == 1) fullDisplacementDirections.back().Add(0.002, Flux(f.displacementFes->GetTrueVSize(), 0.37));
if (gravity) {
densityDirections.push_back(gravity->GetDensityMap().gather(fullDensityDirections.back()));
displacementDirections.push_back(gravity->GetDisplacementMap().gather(fullDisplacementDirections.back()));
}
if (hydro) hydroDirections.push_back(hydro->GetDisplacementMap().gather(fullDisplacementDirections.back()));
omegaDirections.push_back(-0.13 + 0.07 * i);
}
const Apply massFast = [&](int i, mfem::Vector &out) { mass->ApplyCompleteJacobianAction(densityDirections[i], displacementDirections[i], out); };
const Apply massReference = [&](int i, mfem::Vector &out) { mass->ApplyCompleteJacobianActionReference(densityDirections[i], displacementDirections[i], out); };
const Apply angularFast = [&](int i, mfem::Vector &out) { angular->ApplyCompleteJacobianAction(densityDirections[i], displacementDirections[i], omegaDirections[i], out); };
const Apply angularReference = [&](int i, mfem::Vector &out) { angular->ApplyCompleteJacobianActionReference(densityDirections[i], displacementDirections[i], omegaDirections[i], out); };
const Apply hydroFast = [&](int i, mfem::Vector &out) { hydro->ApplyDisplacementJacobianAction(hydroDirections[i], out); };
const Apply hydroReference = [&](int i, mfem::Vector &out) { hydro->ApplyDisplacementJacobianActionCachedReference(hydroDirections[i], out); };
const Apply hydroGeneralReference = [&](int i, mfem::Vector &out) { hydro->ApplyDisplacementJacobianActionReference(hydroDirections[i], out); };
auto compare = [&](const std::string &name, const Apply &fast, const Apply &reference) {
for (int i = 0; i < options.directions; ++i) {
mfem::Vector a, b;
reference(i, b);
fast(i, a);
checks.Add(name + "_direction" + std::to_string(i), ScaledError(a, b), 2.0e-11);
}
};
if (options.check != "none") {
const auto validationStart = Clock::now();
if (mass) checks.Add("mass_rows_initially_lazy", mass->GetJacobianRowPreparationCount(), 0.0);
if (angular) checks.Add("angular_rows_initially_lazy", angular->GetJacobianRowPreparationCount(), 0.0);
const int states = options.check == "full" ? 4 : 1;
for (int state = 0; state < states; ++state) {
const auto density = state >= 2 ? ScalarField(*f.densityFes, 1.13, 0.82) : baseDensity;
const auto displacement = state == 0 && states > 1 ? Displacement(f, 0.0) : baseDisplacement;
prepare(density, displacement, state == 3 ? -0.27 : 0.63, state < 2, state != 3, true, "");
const std::string label = "state" + std::to_string(state);
const auto massCount = mass ? mass->GetJacobianRowPreparationCount() : 0;
const auto angularCount = angular ? angular->GetJacobianRowPreparationCount() : 0;
if (mass) {
compare(label + "_mass", massFast, massReference);
checks.Add(label + "_mass_row_refresh", mass->GetJacobianRowPreparationCount() == massCount + (state != 3) ? 0.0 : 1.0, 0.0);
}
if (angular) {
compare(label + "_angular", angularFast, angularReference);
checks.Add(label + "_angular_row_refresh", angular->GetJacobianRowPreparationCount() == angularCount + (state != 3) ? 0.0 : 1.0, 0.0);
}
if (hydro) {
compare(label + "_hydro_cached", hydroFast, hydroReference);
compare(label + "_hydro_general", hydroFast, hydroGeneralReference);
mfem::Vector zero(hydroDirections[0].Size()), action;
zero = 0.0;
hydro->ApplyDisplacementJacobianAction(zero, action);
checks.Add(label + "_hydro_zero", Norm(action), 0.0);
}
}
prepare(baseDensity, baseDisplacement, 0.63, true, true, true, "");
mfem::Vector massAnalytic, angularAnalytic, hydroAnalytic;
if (mass) massFast(0, massAnalytic);
if (angular) angularFast(0, angularAnalytic);
if (hydro) hydroFast(0, hydroAnalytic);
if (options.check == "full") {
// Check individual scalar columns, including zero input, so
// cancellation in a combined direction cannot hide an error.
if (gravity) {
mfem::Vector zeroDensity(densityDirections[0].Size()), zeroDisplacement(displacementDirections[0].Size());
zeroDensity = 0.0; zeroDisplacement = 0.0;
for (int column = 0; column < 4; ++column) {
const auto &rho = column == 0 ? densityDirections[0] : zeroDensity;
const auto &u = column == 1 ? displacementDirections[0] : zeroDisplacement;
const double omega = column == 2 ? omegaDirections[0] : 0.0;
mfem::Vector a, b;
if (mass) {
mass->ApplyCompleteJacobianAction(rho, u, a);
mass->ApplyCompleteJacobianActionReference(rho, u, b);
checks.Add("mass_column" + std::to_string(column), ScaledError(a, b), 2.0e-11);
}
if (angular) {
angular->ApplyCompleteJacobianAction(rho, u, omega, a);
angular->ApplyCompleteJacobianActionReference(rho, u, omega, b);
checks.Add("angular_column" + std::to_string(column), ScaledError(a, b), 2.0e-11);
}
}
}
if (mass) {
mfem::Vector rhoDual, displacementDual;
mass->ApplyCompleteJacobianTransposeAction(0.73, rhoDual, displacementDual);
mfem::Vector expected(1), actual(1);
expected(0) = 0.73 * massAnalytic(0);
actual(0) = rhoDual * densityDirections[0] + displacementDual * displacementDirections[0];
checks.Add("mass_transpose_identity", ScaledError(actual, expected), 2.0e-11);
const auto count = mass->GetJacobianRowPreparationCount();
targetMass = 1.37; ++md.targetMass.revision;
mass->Prepare({.targetMass = targetMass}, md);
mfem::Vector afterTarget;
massFast(0, afterTarget);
checks.Add("mass_target_only_action", ScaledError(afterTarget, massAnalytic), 2.0e-11);
checks.Add("mass_target_only_reuses_rows", mass->GetJacobianRowPreparationCount() == count ? 0.0 : 1.0, 0.0);
}
for (const double step : {options.finiteDifferenceStep, 0.5 * options.finiteDifferenceStep}) {
mfem::Vector massPlus, angularPlus, hydroPlus, massMinus, angularMinus, hydroMinus;
for (const int sign : {1, -1}) {
mfem::Vector rho(baseDensity), displacement(baseDisplacement);
rho.Add(sign * step, fullDensityDirections[0]);
displacement.Add(sign * step, fullDisplacementDirections[0]);
// Scalar FD varies rho/u/Omega together. Hydro FD keeps
// all fields and rotation fixed except displacement.
prepare(rho, displacement, 0.63, true, true, true, "");
if (mass) mass->BuildResidual(sign == 1 ? massPlus : massMinus);
if (angular) {
++ad.rotation.revision;
angular->Prepare(0.63 + sign * step * omegaDirections[0], ad);
angular->BuildResidual(sign == 1 ? angularPlus : angularMinus);
}
if (hydro) hydro->BuildResidual(sign == 1 ? hydroPlus : hydroMinus);
}
auto difference = [&](const std::string &name, const mfem::Vector &analytic, mfem::Vector &plus, const mfem::Vector &minus) {
plus -= minus; plus /= 2.0 * step;
checks.Add(name + "_finite_difference_" + std::to_string(step), ScaledError(analytic, plus), 2.0e-6);
};
if (mass) difference("mass", massAnalytic, massPlus, massMinus);
if (angular) difference("angular", angularAnalytic, angularPlus, angularMinus);
if (hydro) difference("hydro", hydroAnalytic, hydroPlus, hydroMinus);
}
}
summary << "stellar.validation_seconds=" << std::chrono::duration<double>(Clock::now() - validationStart).count() << '\n';
if (!checks.passed) { summary << "stellar.checks_passed=0\nstellar.benchmark_skipped=1\n"; return 3; }
}
prepare(baseDensity, baseDisplacement, 0.63, true, true, true, "repeat");
auto rows = [&](const std::string &name, auto &operation) {
if (!operation) return 0.0;
const auto start = Clock::now();
const bool supported = operation->PrepareJacobianRows();
const double seconds = std::chrono::duration<double>(Clock::now() - start).count();
metadata << name << ".prepared_rows_supported=" << supported << '\n';
summary << name << ".explicit_row_prepare_seconds=" << seconds << '\n';
return seconds;
};
const double massRowSeconds = rows("mass-normalization", mass);
const double angularRowSeconds = rows("angular-momentum", angular);
if (hydro) {
metadata << "hydrostatic-displacement.fast_elements=" << hydro->GetFastDisplacementElementCount() << '\n'
<< "hydrostatic-displacement.tensor_elements=" << hydro->GetTensorDisplacementElementCount() << '\n'
<< "hydrostatic-displacement.reference=previous_cached_pointwise\n";
std::size_t nq = 0;
hydro->VisitMappedGeometryRules([&](int, const mfem::IntegrationRule &rule) { nq += rule.GetNPoints(); });
metadata << "hydrostatic-displacement.quadrature_points=" << nq << '\n';
}
metadata << "stellar.check_error_metric=norm_difference/max(1,norm_actual,norm_reference)\n";
metadata.flush();
summary << "stellar.checks=" << (options.check == "none" ? "not_run" : "passed") << '\n';
if (mass) Benchmark(options, "mass-normalization", "prepared_reference", massFast, massReference, timings, summary, massRowSeconds);
if (angular) Benchmark(options, "angular-momentum", "prepared_reference", angularFast, angularReference, timings, summary, angularRowSeconds);
if (hydro) Benchmark(options, "hydrostatic-displacement", "cached_pointwise_reference", hydroFast, hydroReference, timings, summary);
return 0;
}
bool IsLinearDerivativeSystem(const std::string &name) {
return name == "application-optimizations" || name == "linear-derivatives" || name == "barotropic-closure" ||
name == "barotropic-enthalpy" || name == "barotropic-displacement" ||
name == "pressure-displacement" || name == "pressure-enthalpy" ||
name == "gravity-force" || name == "rotation-force" || name == "preparation-optimizations";
}
int RunLinearDerivatives(const Options &options, fem::FEM &f, std::ofstream &metadata,
std::ofstream &summary, std::ofstream &timings, CheckReport &checks) {
const bool preparationGroup = options.system == "preparation-optimizations";
const bool group = preparationGroup || options.system == "all" || options.system == "linear-derivatives";
const bool applicationGroup = options.system == "application-optimizations";
const bool runClosure = applicationGroup || (group && !preparationGroup) || options.system.starts_with("barotropic-");
const bool runPressure = group || options.system.starts_with("pressure-");
const bool runGravity = group || options.system == "gravity-force";
const bool runRotation = applicationGroup || group || options.system == "rotation-force";
using Schema = utils::domain::CoreEnvelopeVacuumDomainSchema;
const auto densityMap = field::make_field_dof_map<field::Density, Schema>(*f.densityFes);
const auto enthalpyMap = field::make_field_dof_map<field::Enthalpy, Schema>(*f.enthalpyFes);
const auto displacementMap = field::make_field_dof_map<field::Displacement, Schema>(*f.displacementFes);
const eos::Polytrope equationOfState(options.polytropicIndex, 1.0);
using GravityContext = operators::context::gravity_field::GravityFieldLinearizationContext;
const auto constructStart = Clock::now();
std::unique_ptr<GravityContext> gravityContext;
std::unique_ptr<operators::PreparedBarotropicClosureOperator> closure;
std::unique_ptr<operators::PreparedPressureForceOperator> pressure, pressureReference;
std::unique_ptr<operators::PreparedGravityDisplacementForceOperator> gravity, gravityReference;
std::unique_ptr<operators::PreparedRotationalDisplacementForceOperator> rotation, rotationReference;
if (runClosure) closure = std::make_unique<operators::PreparedBarotropicClosureOperator>(
f, *f.domainMapperStateless, equationOfState);
if (runPressure) {
pressure = std::make_unique<operators::PreparedPressureForceOperator>(f, *f.domainMapperStateless, equationOfState);
pressureReference = std::make_unique<operators::PreparedPressureForceOperator>(f, *f.domainMapperStateless, equationOfState);
}
if (runGravity) {
gravityContext = std::make_unique<GravityContext>(f, *f.domainMapperStateless);
gravity = std::make_unique<operators::PreparedGravityDisplacementForceOperator>(
f, *f.domainMapperStateless, *gravityContext);
gravityReference = std::make_unique<operators::PreparedGravityDisplacementForceOperator>(
f, *f.domainMapperStateless, *gravityContext);
}
if (runRotation) {
rotation = std::make_unique<operators::PreparedRotationalDisplacementForceOperator>(f, *f.domainMapperStateless);
rotationReference = std::make_unique<operators::PreparedRotationalDisplacementForceOperator>(f, *f.domainMapperStateless);
}
summary << "linear.operator_construction_seconds=" <<
std::chrono::duration<double>(Clock::now() - constructStart).count() << '\n';
const auto baseDensity = densityMap.gather(ScalarField(*f.densityFes, 0.94, 0.31));
const auto baseEnthalpy = enthalpyMap.gather(ScalarField(*f.enthalpyFes, 0.83, 0.27));
const auto baseDisplacement = displacementMap.gather(Displacement(f, options.deformation));
mfem::Vector baseFlux;
mfem::Vector potential;
if (gravityContext) {
baseFlux = gravityContext->GetGravityGradientMap().gather(Flux(f.gravityFluxFes->GetTrueVSize(), 0.29));
baseFlux *= 0.04;
potential = gravityContext->GetGravityPotentialMap().gather(
ScalarField(*f.gravityPotentialFes, -0.38, 0.61));
}
std::vector<mfem::Vector> rhoDirections, hDirections, uDirections, gDirections;
for (int i = 0; i < options.directions; ++i) {
rhoDirections.push_back(densityMap.gather(ScalarField(*f.densityFes, 0.17, 0.4 + i)));
hDirections.push_back(enthalpyMap.gather(ScalarField(*f.enthalpyFes, 0.11, 0.6 + i)));
auto u = Displacement(f, 0.43, 0.2 + 0.7 * i);
if (i == 1) u.Add(0.002, Flux(f.displacementFes->GetTrueVSize(), 0.37));
uDirections.push_back(displacementMap.gather(u));
if (gravityContext) {
gDirections.push_back(gravityContext->GetGravityGradientMap().gather(
Flux(f.gravityFluxFes->GetTrueVSize(), 0.81 + i)));
gDirections.back() *= 0.03;
}
}
mfem::Vector zeroRho(baseDensity.Size()), zeroH(baseEnthalpy.Size()),
zeroU(baseDisplacement.Size()), zeroG(baseFlux.Size());
zeroRho = 0.0; zeroH = 0.0; zeroU = 0.0; zeroG = 0.0;
// The last direction is reserved for checks and is never timed.
rhoDirections.push_back(zeroRho); hDirections.push_back(zeroH); uDirections.push_back(zeroU);
if (gravityContext) gDirections.push_back(zeroG);
std::uint64_t densityRevision = 1, enthalpyRevision = 1, displacementRevision = 1,
fluxRevision = 1, rotationRevision = 1;
const bool needPreparationReferences = options.check != "none" || options.implementation != "optimized";
bool prepareReferences = needPreparationReferences;
auto prepare = [&](const mfem::Vector &rho, const mfem::Vector &h, const mfem::Vector &u,
const mfem::Vector &g, double omega, unsigned changed, const std::string &label) {
if (changed & 1) ++densityRevision;
if (changed & 2) ++enthalpyRevision;
if (changed & 4) ++displacementRevision;
if (changed & 8) ++fluxRevision;
if (changed & 16) ++rotationRevision;
auto time = [&](const std::string &name, auto &&operation) {
const auto start = Clock::now(); operation();
if (!label.empty()) summary << name << ".prepare." << label << "_seconds="
<< std::chrono::duration<double>(Clock::now() - start).count() << '\n';
};
if (gravityContext) time("linear.shared-gravity", [&] {
gravityContext->Prepare(
{.density = rho, .displacement = u, .gravity_gradient = g, .gravity_potential = potential},
{.discretization = {1}, .displacement = {displacementRevision}, .density = {densityRevision},
.gravity_gradient = {fluxRevision}, .gravity_potential = {1}});
});
if (closure) time("barotropic-closure", [&] {
closure->Prepare({.density = rho, .enthalpy = h, .displacement = u},
{.discretization = {1, 1}, .density = {2, densityRevision}, .enthalpy = {3, enthalpyRevision},
.displacement = {4, displacementRevision}});
});
if (pressure) time("pressure-displacement", [&] {
pressure->Prepare({.enthalpy = h, .displacement = u},
{.discretization = {1, 1}, .enthalpy = {3, enthalpyRevision}, .displacement = {4, displacementRevision}});
});
if (pressureReference && prepareReferences) time("pressure-displacement.reference", [&] {
pressureReference->PrepareCachedReference({.enthalpy = h, .displacement = u},
{.discretization = {1, 1}, .enthalpy = {3, enthalpyRevision}, .displacement = {4, displacementRevision}});
});
if (gravity) time("gravity-force", [&] { gravity->Prepare(); });
if (gravityReference && prepareReferences) time("gravity-force.reference", [&] { gravityReference->PrepareReference(); });
if (rotation) time("rotation-force", [&] {
mfem::Vector velocity(3), center(3);
velocity(0) = omega / std::sqrt(14.0); velocity(1) = -2.0 * velocity(0); velocity(2) = 3.0 * velocity(0);
center(0) = 0.12; center(1) = -0.07; center(2) = 0.03;
rotation->Prepare({.density = rho, .displacement = u},
{.discretization = {1, 1}, .density = {2, densityRevision},
.displacement = {4, displacementRevision}, .rotation = {5, rotationRevision}},
physics::RigidRotation(velocity, center));
});
if (rotationReference && prepareReferences) time("rotation-force.reference", [&] {
mfem::Vector velocity(3), center(3);
velocity(0) = omega / std::sqrt(14.0); velocity(1) = -2.0 * velocity(0); velocity(2) = 3.0 * velocity(0);
center(0) = 0.12; center(1) = -0.07; center(2) = 0.03;
rotationReference->PrepareReference({.density = rho, .displacement = u},
{.discretization = {1, 1}, .density = {2, densityRevision},
.displacement = {4, displacementRevision}, .rotation = {5, rotationRevision}},
physics::RigidRotation(velocity, center));
});
};
prepare(baseDensity, baseEnthalpy, baseDisplacement, baseFlux, 0.63, 31, "initial");
if (gravity) metadata << "gravity_tensor_elements=" << gravity->GetTensorElementCount() << '\n';
struct System {
std::string name;
Apply fast, reference;
std::function<void(mfem::Vector &)> residual;
unsigned columns; // density=1, enthalpy=2, displacement=4, gravity=8
};
std::vector<System> systems;
auto add = [&](System system) {
if (applicationGroup || group || options.system == system.name) systems.push_back(std::move(system));
};
if (closure) {
auto residual = [&](mfem::Vector &out) { closure->BuildResidual(out); };
add({"barotropic-closure",
[&](int i, mfem::Vector &out) { closure->Mult(rhoDirections[i], hDirections[i], uDirections[i], out); },
[&](int i, mfem::Vector &out) { closure->MultCachedReference(rhoDirections[i], hDirections[i], uDirections[i], out); },
residual, 7});
add({"barotropic-enthalpy",
[&](int i, mfem::Vector &out) { closure->ApplyEnthalpyJacobianAction(hDirections[i], out); },
[&](int i, mfem::Vector &out) { closure->ApplyEnthalpyJacobianAction(hDirections[i], out, false); },
residual, 2});
add({"barotropic-displacement",
[&](int i, mfem::Vector &out) { closure->ApplyDisplacementJacobianAction(uDirections[i], out); },
[&](int i, mfem::Vector &out) { closure->ApplyDisplacementJacobianAction(uDirections[i], out, false); },
residual, 4});
}
if (pressure) add({"pressure-displacement",
[&](int i, mfem::Vector &out) { pressure->ApplyDisplacementJacobianAction(uDirections[i], out); },
[&](int i, mfem::Vector &out) { pressure->ApplyDisplacementJacobianActionCachedReference(uDirections[i], out); },
[&](mfem::Vector &out) { pressure->BuildResidual(out); }, 4});
if (pressure) add({"pressure-enthalpy",
[&](int i, mfem::Vector &out) { pressure->ApplyEnthalpyJacobianAction(hDirections[i], out); },
[&](int i, mfem::Vector &out) { pressureReference->ApplyEnthalpyJacobianAction(hDirections[i], out); },
[&](mfem::Vector &out) { pressure->BuildResidual(out); }, 2});
if (gravity) add({"gravity-force",
[&](int i, mfem::Vector &out) { gravity->ApplyCompleteJacobianAction(rhoDirections[i], uDirections[i], gDirections[i], out); },
[&](int i, mfem::Vector &out) { gravity->ApplyCompleteJacobianActionCachedReference(rhoDirections[i], uDirections[i], gDirections[i], out); },
[&](mfem::Vector &out) { gravity->BuildResidual(out); }, 13});
if (rotation) add({"rotation-force",
[&](int i, mfem::Vector &out) { rotation->ApplyCompleteJacobianAction(rhoDirections[i], uDirections[i], out); },
[&](int i, mfem::Vector &out) { rotation->ApplyCompleteJacobianActionCachedReference(rhoDirections[i], uDirections[i], out); },
[&](mfem::Vector &out) { rotation->BuildResidual(out); }, 5});
auto compare = [&](const std::string &label, const mfem::Vector &a, const mfem::Vector &b) {
checks.Add(label, ScaledError(a, b), 2.0e-11);
if (std::max(Norm(a), Norm(b)) > 1.0e-14)
checks.Add(label + "_relative", RelativeError(a, b), 2.0e-10);
};
if (options.check != "none") {
const auto start = Clock::now();
const int states = options.check == "full" ? 5 : 1;
for (int state = 0; state < states; ++state) {
const auto rho = state >= 2 ? densityMap.gather(ScalarField(*f.densityFes, 1.13, 0.82)) : baseDensity;
const auto h = state >= 2 ? enthalpyMap.gather(ScalarField(*f.enthalpyFes, 0.91, 0.42)) : baseEnthalpy;
const auto &u = state == 0 && states > 1 ? zeroU : baseDisplacement;
auto g = baseFlux;
if (state >= 3) g *= -0.71;
const unsigned revisions = state == 0 ? 31u : state == 1 ? 4u : state == 2 ? 3u : state == 3 ? 24u : 16u;
const double omega = state == 4 ? 0.0 : state == 3 ? -0.27 : 0.63;
prepare(rho, h, u, g, omega, revisions, "");
auto compareResidual = [&](const std::string &name, auto &fast, auto &reference) {
if (!fast) return;
mfem::Vector a, b;
fast->BuildResidual(a); reference->BuildResidual(b);
compare("state" + std::to_string(state) + "_" + name + "_prepared_residual", a, b);
};
compareResidual("pressure", pressure, pressureReference);
compareResidual("gravity", gravity, gravityReference);
compareResidual("rotation", rotation, rotationReference);
if (pressure) {
for (int i = 0; i <= options.directions; ++i) {
mfem::Vector a, b;
pressure->ApplyEnthalpyJacobianAction(hDirections[i], a);
pressureReference->ApplyEnthalpyJacobianAction(hDirections[i], b);
compare("state" + std::to_string(state) + "_pressure_assembled_enthalpy_" + std::to_string(i), a, b);
}
}
for (const auto &system : systems) {
for (int i = 0; i <= options.directions; ++i) {
mfem::Vector a, b;
system.fast(i, a); system.reference(i, b);
compare("state" + std::to_string(state) + "_" + system.name + "_direction" + std::to_string(i), a, b);
if (system.name == "pressure-displacement") {
pressure->ApplyDisplacementJacobianActionReference(uDirections[i], b);
compare("state" + std::to_string(state) + "_pressure_general_" + std::to_string(i), a, b);
}
if (system.name == "gravity-force") {
gravity->ApplyCompleteJacobianActionReference(rhoDirections[i], uDirections[i], gDirections[i], b);
compare("state" + std::to_string(state) + "_gravity_general_" + std::to_string(i), a, b);
}
if (system.name == "rotation-force") {
rotation->ApplyCompleteJacobianActionReference(rhoDirections[i], uDirections[i], b);
compare("state" + std::to_string(state) + "_rotation_general_" + std::to_string(i), a, b);
}
if (system.name.starts_with("barotropic-")) {
closure->MultReference(system.columns & 1 ? rhoDirections[i] : zeroRho,
system.columns & 2 ? hDirections[i] : zeroH,
system.columns & 4 ? uDirections[i] : zeroU, b);
compare("state" + std::to_string(state) + "_" + system.name + "_general_" + std::to_string(i), a, b);
}
if (i == options.directions) checks.Add(system.name + "_zero", Norm(a), 0.0);
}
}
}
prepare(baseDensity, baseEnthalpy, baseDisplacement, baseFlux, 0.63, 31, "");
if (options.check == "full") {
// Reference preparation was checked at every state above.
// Finite differences below need only the production residual.
prepareReferences = false;
// Compare separate columns against the old complete action.
// Individual gravity/rotation column methods remain independent
// stateless kernels and serve as an additional cross-check.
for (unsigned column : {1u, 2u, 4u, 8u}) {
const auto &rho = column == 1 ? rhoDirections[0] : zeroRho;
const auto &h = column == 2 ? hDirections[0] : zeroH;
const auto &u = column == 4 ? uDirections[0] : zeroU;
const auto &g = column == 8 && gravity ? gDirections[0] : zeroG;
mfem::Vector a, b, independent;
if (closure) {
closure->Mult(rho, h, u, a); closure->MultReference(rho, h, u, b);
compare("closure_column" + std::to_string(column), a, b);
}
if (gravity && column != 2) {
gravity->ApplyCompleteJacobianAction(rho, u, g, a);
gravity->ApplyCompleteJacobianActionReference(rho, u, g, b);
compare("gravity_column" + std::to_string(column), a, b);
if (column == 1) gravity->ApplyDensityJacobianAction(rho, independent);
if (column == 4) gravity->ApplyDisplacementJacobianAction(u, independent);
if (column == 8) gravity->ApplyGravityGradientJacobianAction(g, independent);
compare("gravity_stateless_column" + std::to_string(column), a, independent);
}
if (rotation && (column == 1 || column == 4)) {
rotation->ApplyCompleteJacobianAction(rho, u, a);
rotation->ApplyCompleteJacobianActionReference(rho, u, b);
compare("rotation_column" + std::to_string(column), a, b);
if (column == 1) rotation->ApplyDensityJacobianAction(rho, independent);
else rotation->ApplyDisplacementJacobianAction(u, independent);
compare("rotation_stateless_column" + std::to_string(column), a, independent);
}
}
for (const auto &system : systems) {
// The complete closure FD plus independent column checks
// covers its two couplings in a grouped run.
if (group && (system.name == "barotropic-enthalpy" || system.name == "barotropic-displacement")) continue;
prepare(baseDensity, baseEnthalpy, baseDisplacement, baseFlux, 0.63, 31, "");
mfem::Vector analytic;
system.fast(0, analytic);
for (double step : {options.finiteDifferenceStep, 0.5 * options.finiteDifferenceStep}) {
mfem::Vector plus, minus;
for (int sign : {1, -1}) {
auto rho = baseDensity, h = baseEnthalpy, u = baseDisplacement, g = baseFlux;
if (system.columns & 1) rho.Add(sign * step, rhoDirections[0]);
if (system.columns & 2) h.Add(sign * step, hDirections[0]);
if (system.columns & 4) u.Add(sign * step, uDirections[0]);
if (system.columns & 8) g.Add(sign * step, gDirections[0]);
prepare(rho, h, u, g, 0.63, 31, "");
system.residual(sign == 1 ? plus : minus);
}
plus -= minus; plus /= 2.0 * step;
checks.Add(system.name + "_finite_difference_" + std::to_string(step),
ScaledError(analytic, plus), 2.0e-6);
if (Norm(analytic) > 1.0e-12)
checks.Add(system.name + "_finite_difference_relative_" + std::to_string(step),
RelativeError(analytic, plus), 2.0e-6);
}
}
}
prepareReferences = needPreparationReferences;
summary << "linear.validation_seconds=" << std::chrono::duration<double>(Clock::now() - start).count() << '\n';
if (!checks.passed) { summary << "linear.checks=failed\n"; return 3; }
}
prepare(baseDensity, baseEnthalpy, baseDisplacement, baseFlux, 0.63, 31, "repeat");
// Unchanged dependencies must preserve the prepared state and results.
const auto count = [&] {
return (closure ? closure->GetPreparationCount() : 0) +
(pressure ? pressure->GetResidualPreparationCount() : 0) +
(gravity ? gravity->GetResidualPreparationCount() : 0) +
(rotation ? rotation->GetResidualPreparationCount() : 0);
};
const auto beforeNoop = count();
prepare(baseDensity, baseEnthalpy, baseDisplacement, baseFlux, 0.63, 0, "unchanged");
if (options.check != "none") {
checks.Add("linear_unchanged_dependencies_reuse", count() == beforeNoop ? 0.0 : 1.0, 0.0);
if (!checks.passed) { summary << "linear.checks=failed\n"; return 3; }
}
if (gravity) metadata << "gravity-force.fast_elements=" << gravity->GetFastElementCount() << '\n';
if (rotation) metadata << "rotation-force.fast_elements=" << rotation->GetFastElementCount()
<< "\nrotation_tensor_elements=" << rotation->GetTensorElementCount() << '\n';
if (pressure) metadata << "pressure_tensor_elements=" << pressure->GetTensorElementCount() << '\n';
if (closure) metadata << "closure_tensor_elements=" << closure->GetTensorElementCount() << '\n';
metadata << "linear.polytropic_index=" << options.polytropicIndex << "\nlinear.polytropic_constant=1\n"
<< "linear.check_error_metric=norm_difference/max(1,norm_actual,norm_reference)\n"
<< "linear.relative_check_metric=norm_difference/max(norm_actual,norm_reference,1e-30)\n";
metadata.flush();
summary << "linear.checks=" << (options.check == "none" ? "not_run" : "passed") << '\n';
for (const auto &system : systems)
Benchmark(options, system.name, (system.name == "gravity-force" || system.name == "rotation-force" || system.name == "pressure-displacement" || system.name.starts_with("barotropic-")) ?
"cached_pointwise_reference" : "prepared_reference",
system.fast, system.reference, timings, summary);
// Use the same physical state for both paths, with fresh dependency
// revisions on each pair. Shared gravity context preparation is excluded.
auto requireRebuild = [](auto &operation, auto &&prepareOperation) {
const auto before = operation.GetResidualPreparationCount();
prepareOperation();
if (operation.GetResidualPreparationCount() != before + 1)
throw std::runtime_error("Preparation benchmark timed an unchanged-state reuse.");
};
if (pressure) {
auto dependencies = operators::context::pressure_force::PressureForceDependencies{
.discretization = {1, 1}, .enthalpy = {3, enthalpyRevision}, .displacement = {4, displacementRevision}};
const auto state = operators::context::pressure_force::PressureForceStateView{
.enthalpy = baseEnthalpy, .displacement = baseDisplacement};
BenchmarkPreparation(options, "pressure-preparation",
[&] { dependencies.enthalpy.revision = ++enthalpyRevision;
dependencies.displacement.revision = ++displacementRevision; },
[&] { requireRebuild(*pressure, [&] { pressure->Prepare(state, dependencies); }); },
[&] { requireRebuild(*pressureReference, [&] { pressureReference->PrepareCachedReference(state, dependencies); }); },
[&](mfem::Vector &out) { pressure->BuildResidual(out); },
[&](mfem::Vector &out) { pressureReference->BuildResidual(out); }, timings, summary);
}
if (gravity) {
BenchmarkPreparation(options, "gravity-force-preparation",
[&] {
gravityContext->Prepare(
{.density = baseDensity, .displacement = baseDisplacement,
.gravity_gradient = baseFlux, .gravity_potential = potential},
{.discretization = {1}, .displacement = {++displacementRevision},
.density = {++densityRevision}, .gravity_gradient = {++fluxRevision}, .gravity_potential = {1}});
},
[&] { requireRebuild(*gravity, [&] { gravity->Prepare(); }); },
[&] { requireRebuild(*gravityReference, [&] { gravityReference->PrepareReference(); }); },
[&](mfem::Vector &out) { gravity->BuildResidual(out); },
[&](mfem::Vector &out) { gravityReference->BuildResidual(out); }, timings, summary);
}
if (rotation) {
auto dependencies = operators::context::rotational_displacement_force::RotationalDisplacementForceDependencies{
.discretization = {1, 1}, .density = {2, densityRevision},
.displacement = {4, displacementRevision}, .rotation = {5, rotationRevision}};
const auto state = operators::context::rotational_displacement_force::RotationalDisplacementForceStateView{
.density = baseDensity, .displacement = baseDisplacement};
mfem::Vector velocity(3), center(3);
velocity(0) = 0.63 / std::sqrt(14.0); velocity(1) = -2.0 * velocity(0); velocity(2) = 3.0 * velocity(0);
center(0) = 0.12; center(1) = -0.07; center(2) = 0.03;
const physics::RigidRotation rigidRotation(velocity, center);
BenchmarkPreparation(options, "rotation-force-preparation",
[&] { dependencies.density.revision = ++densityRevision;
dependencies.displacement.revision = ++displacementRevision;
dependencies.rotation.revision = ++rotationRevision; },
[&] { requireRebuild(*rotation, [&] { rotation->Prepare(state, dependencies, rigidRotation); }); },
[&] { requireRebuild(*rotationReference, [&] { rotationReference->PrepareReference(state, dependencies, rigidRotation); }); },
[&](mfem::Vector &out) { rotation->BuildResidual(out); },
[&](mfem::Vector &out) { rotationReference->BuildResidual(out); }, timings, summary);
}
return 0;
}
int RunNormalizationDiagonals(
const Options &options, fem::FEM &f, std::ofstream &summary, std::ofstream &timings, CheckReport &checks
) {
mfem::Array<int> stellarMarker(f.mesh->attributes.Max());
stellarMarker = 0;
for (int e = 0; e < f.mesh->GetNE(); ++e)
if (!f.domainMapperStateless->IsCompactifiedElement(*f.mesh->GetElementTransformation(e)))
stellarMarker[f.mesh->GetAttribute(e)-1] = 1;
const auto surfaceMap = field::make_stellar_surface_scalar_dof_map(*f.surfaceDeformationFes);
struct Diagonal { std::string name; std::function<mfem::Vector(bool)> assemble; };
const std::vector<Diagonal> systems{
{"normalization-hdiv", [&](bool ref) { return normalization::AssembleReferenceHDivMassDiagonal(*f.gravityFluxFes, ref); }},
{"normalization-potential", [&](bool ref) { return normalization::AssembleReferenceScalarMassDiagonal(*f.gravityPotentialFes, nullptr, ref); }},
{"normalization-density", [&](bool ref) { return normalization::AssembleReferenceScalarMassDiagonal(*f.densityFes, &stellarMarker, ref); }},
{"normalization-enthalpy", [&](bool ref) { return normalization::AssembleReferenceScalarMassDiagonal(*f.enthalpyFes, &stellarMarker, ref); }},
{"normalization-surface", [&](bool ref) { return normalization::AssembleReferenceSurfaceMassDiagonal(f, surfaceMap, ref); }}
};
for (const auto &system : systems) {
mfem::Vector fast, reference;
if (options.check != "none") {
const auto start = Clock::now();
fast = system.assemble(false);
summary << system.name << ".first_optimized_seconds=" <<
std::chrono::duration<double>(Clock::now()-start).count() << '\n';
reference = system.assemble(true);
summary << system.name << ".optimized_diagonal_norm=" << Norm(fast) << '\n'
<< system.name << ".reference_diagonal_norm=" << Norm(reference) << '\n';
checks.Add(system.name + "_exact_true_diagonal", RelativeError(fast, reference), 2.0e-11);
if (system.name == "normalization-surface") {
int nonpositiveEntries = 0;
for (int i = 0; i < fast.Size(); ++i)
if (fast(i) <= 0.0) ++nonpositiveEntries;
checks.Add(system.name + "_positive_surface_diagonal", nonpositiveEntries, 0.0);
}
}
if (!checks.passed) { summary << "checks_passed=0\nbenchmark_skipped=1\n"; return 3; }
BenchmarkPreparation(options, system.name, [] {},
[&] { fast = system.assemble(false); }, [&] { reference = system.assemble(true); },
[&](mfem::Vector &out) { out = fast; }, [&](mfem::Vector &out) { out = reference; }, timings, summary);
}
return 0;
}
int RunGeometryPreparation(
const Options &options, fem::FEM &f, std::ofstream &metadata, std::ofstream &summary,
std::ofstream &timings, CheckReport &checks
) {
const auto &mapper = *f.domainMapperStateless;
struct RestoreCacheMode {
const mapping::DomainMapper &mapper;
bool enabled;
~RestoreCacheMode() { mapper.SetGeometryCacheEnabled(enabled); }
} restore{mapper, mapper.IsGeometryCacheEnabled()};
operators::PreparedMappedHDivMassOperator fast(f, mapper), reference(f, mapper);
const auto flux = Flux(f.gravityFluxFes->GetTrueVSize(), 0.37);
const auto reducedFlux = fast.GetFluxMap().gather(flux);
const auto direction = Displacement(f, 0.21, 0.32);
const std::vector<mfem::Vector> states{
fast.GetDisplacementMap().gather(Displacement(f, 0.0)),
fast.GetDisplacementMap().gather(Displacement(f, options.deformation)),
fast.GetDisplacementMap().gather(Displacement(f, 0.5*options.deformation, 0.4))
};
int state = 1;
auto prepare = [&](auto &op, bool cache) {
mapper.SetGeometryCacheEnabled(cache);
const auto count = op.GetPreparationCount();
op.Prepare(states[state]);
if (op.GetPreparationCount() != count + 1)
throw std::runtime_error("Geometry preparation benchmark did not rebuild its operator.");
};
mapper.ClearGeometryCache();
const auto start = Clock::now();
prepare(fast, true);
summary << "geometry-preparation.cold_cache_prepare_seconds=" <<
std::chrono::duration<double>(Clock::now()-start).count() << '\n';
if (options.check != "none") {
for (state = 0; state < static_cast<int>(states.size()); ++state) {
prepare(fast, true); prepare(reference, false);
mfem::Vector a, b;
fast.Mult(reducedFlux, a); reference.Mult(reducedFlux, b);
checks.Add("geometry_state" + std::to_string(state) + "_mass", RelativeError(a,b), 2.0e-11);
fast.MultDisplacementVariationTrue(flux, direction, a);
reference.MultDisplacementVariationTrue(flux, direction, b);
checks.Add("geometry_state" + std::to_string(state) + "_derivative", RelativeError(a,b), 2.0e-11);
}
}
if (!checks.passed) { summary << "checks_passed=0\nbenchmark_skipped=1\n"; return 3; }
state = 0;
const auto before = mapper.GetGeometryCacheStatistics();
BenchmarkPreparation(options, "geometry-preparation",
[&] { state = (state + 1) % static_cast<int>(states.size()); },
[&] { prepare(fast, true); }, [&] { prepare(reference, false); },
[&](mfem::Vector &out) { fast.Mult(reducedFlux, out); },
[&](mfem::Vector &out) { reference.Mult(reducedFlux, out); }, timings, summary);
const auto after = mapper.GetGeometryCacheStatistics();
metadata << "geometry_cache_points=" << after.points << "\ngeometry_cache_hits=" << after.hits
<< "\ngeometry_cache_misses=" << after.misses << '\n';
if (options.check != "none" && options.prepareSamples > 0 && options.implementation != "reference") {
checks.Add("geometry_warm_preparation_new_fixed_points", double(after.misses-before.misses), 0.0);
checks.Add("geometry_warm_preparation_reuses_fixed_points", after.hits > before.hits ? 0.0 : 1.0, 0.0);
}
return checks.passed ? 0 : 3;
}
int RunGeometryInspection(
const Options &options, fem::FEM &f, std::ofstream &metadata,
std::ofstream &summary, std::ofstream &timings, CheckReport &checks
) {
fem::PreparedGeometryInspector inspector;
std::vector<mfem::Vector> states{
Displacement(f, 0.0), Displacement(f, options.deformation),
Displacement(f, 0.5 * options.deformation, 0.4)
};
const auto start = Clock::now();
double fast = inspector.Inspect(*f.displacementFes, states[1]), reference = 0.0;
summary << "geometry-inspection.cold_cache_seconds="
<< std::chrono::duration<double>(Clock::now()-start).count() << '\n';
if (options.check != "none") {
for (int i = 0; i < static_cast<int>(states.size()); ++i) {
fast = inspector.Inspect(*f.displacementFes, states[i]);
reference = inspector.Inspect(*f.displacementFes, states[i], true);
checks.Add("geometry_inspection_state" + std::to_string(i),
std::abs(fast-reference) / std::max(1.0, std::abs(reference)), 2.0e-11);
checks.Add("geometry_inspection_state" + std::to_string(i) + "_positive",
fast > 0.0 ? 0.0 : 1.0, 0.0);
}
if (options.check == "full") {
mfem::ParGridFunction folded(f.displacementFes.get());
mfem::VectorFunctionCoefficient fold(3, [](const mfem::Vector &x, mfem::Vector &u) {
u.SetSize(3); u = 0.0; u(0) = -2.0*x(0);
});
folded.ProjectCoefficient(fold);
mfem::Vector foldedTrue; folded.GetTrueDofs(foldedTrue);
fast = inspector.Inspect(*f.displacementFes, foldedTrue);
reference = inspector.Inspect(*f.displacementFes, foldedTrue, true);
checks.Add("geometry_inspection_fold_matches_reference", std::abs(fast-reference), 2.0e-11);
checks.Add("geometry_inspection_fold_rejected", fast < 0.0 ? 0.0 : 1.0, 0.0);
auto nonfinite = states[1];
nonfinite(0) = std::numeric_limits<double>::quiet_NaN();
fast = inspector.Inspect(*f.displacementFes, nonfinite);
reference = inspector.Inspect(*f.displacementFes, nonfinite, true);
checks.Add("geometry_inspection_nonfinite_rejected",
!std::isfinite(fast) && !std::isfinite(reference) ? 0.0 : 1.0, 0.0);
}
}
if (!checks.passed) { summary << "checks_passed=0\nbenchmark_skipped=1\n"; return 3; }
int state = 0;
const auto builds = inspector.GetGeometryBuildCount();
BenchmarkPreparation(options, "geometry-inspection",
[&] { state = (state + 1) % static_cast<int>(states.size()); },
[&] { fast = inspector.Inspect(*f.displacementFes, states[state]); },
[&] { reference = inspector.Inspect(*f.displacementFes, states[state], true); },
[&](mfem::Vector &out) { out.SetSize(1); out(0) = fast; },
[&](mfem::Vector &out) { out.SetSize(1); out(0) = reference; }, timings, summary);
metadata << "geometry_inspection_tensor_elements=" << inspector.GetTensorElementCount()
<< "\ngeometry_inspection_geometry_builds=" << inspector.GetGeometryBuildCount() << '\n';
if (options.check != "none")
checks.Add("geometry_inspection_warm_cache_reuse", double(inspector.GetGeometryBuildCount()-builds), 0.0);
summary << "geometry-inspection.checks=" << (options.check == "none" ? "not_run" : checks.passed ? "passed" : "failed") << '\n';
return checks.passed ? 0 : 3;
}
int RunGravitySourcePreparation(
const Options &options, fem::FEM &f, std::ofstream &summary,
std::ofstream &timings, CheckReport &checks
) {
operators::PreparedMappedGravitySourceOperator fast(f, *f.domainMapperStateless),
reference(f, *f.domainMapperStateless);
const auto rhoTrue = ScalarField(*f.densityFes, 0.94, 0.31);
const auto rho = fast.GetDensityMap().gather(rhoTrue);
const auto potential = fast.GetPotentialMap().gather(ScalarField(*f.gravityPotentialFes, -0.38, 0.61));
const auto direction = Displacement(f, 0.21, 0.32);
const std::vector<mfem::Vector> states{
fast.GetDisplacementMap().gather(Displacement(f, 0.0)),
fast.GetDisplacementMap().gather(Displacement(f, options.deformation)),
fast.GetDisplacementMap().gather(Displacement(f, 0.5 * options.deformation, 0.4))
};
int state = 1;
auto prepare = [&](auto &op, bool oldWeights) {
const auto count = op.GetPreparationCount();
if (oldWeights) op.PrepareReference(states[state]); else op.Prepare(states[state]);
if (op.GetPreparationCount() != count + 1)
throw std::runtime_error("Gravity-source preparation benchmark did not rebuild.");
};
const auto start = Clock::now();
prepare(fast, false);
summary << "gravity-source-preparation.first_optimized_seconds="
<< std::chrono::duration<double>(Clock::now()-start).count() << '\n';
if (options.check != "none") {
for (state = 0; state < static_cast<int>(states.size()); ++state) {
prepare(fast, false); prepare(reference, true);
mfem::Vector a, b;
fast.Mult(rho, a); reference.Mult(rho, b);
checks.Add("source_state" + std::to_string(state) + "_mass", RelativeError(a,b), 2.0e-11);
fast.MultTranspose(potential, a); reference.MultTranspose(potential, b);
checks.Add("source_state" + std::to_string(state) + "_transpose", RelativeError(a,b), 2.0e-11);
fast.MultDisplacementVariationTrue(rhoTrue, direction, a);
reference.MultDisplacementVariationTrue(rhoTrue, direction, b);
checks.Add("source_state" + std::to_string(state) + "_derivative", RelativeError(a,b), 2.0e-11);
}
}
if (!checks.passed) { summary << "checks_passed=0\nbenchmark_skipped=1\n"; return 3; }
state = 0;
BenchmarkPreparation(options, "gravity-source-preparation",
[&] { state = (state + 1) % static_cast<int>(states.size()); },
[&] { prepare(fast, false); }, [&] { prepare(reference, true); },
[&](mfem::Vector &out) { fast.Mult(rho, out); },
[&](mfem::Vector &out) { reference.Mult(rho, out); }, timings, summary);
summary << "gravity-source-preparation.checks=" << (options.check == "none" ? "not_run" : "passed") << '\n';
return 0;
}
int RunSourceDisplacement(const Options &options, fem::FEM &f, std::ofstream &metadata,
std::ofstream &summary, std::ofstream &timings, CheckReport &checks) {
operators::PreparedMappedGravitySourceOperator source(f, *f.domainMapperStateless);
const auto rho = ScalarField(*f.densityFes, 0.94, 0.31);
const auto reducedRho = source.GetDensityMap().gather(rho);
const auto base = Displacement(f, options.deformation);
std::vector<mfem::Vector> directions;
for (int i = 0; i < options.directions; ++i) {
directions.push_back(Displacement(f, 0.43, 0.2+0.7*i));
if (i == 1) directions.back().Add(0.002, Flux(f.displacementFes->GetTrueVSize(), 0.37));
}
auto prepare = [&](const mfem::Vector &u) { source.Prepare(source.GetDisplacementMap().gather(u)); };
const Apply fast = [&](int i, mfem::Vector &out) { source.MultDisplacementVariationTrue(rho, directions[i], out); };
const Apply reference = [&](int i, mfem::Vector &out) {
source.MultDisplacementVariationTrueReference(rho, directions[i], out);
};
if (options.check != "none") {
for (int state = 0; state < (options.check == "full" ? 3 : 1); ++state) {
const auto u = state == 0 ? base : Displacement(f, state == 1 ? 0.0 : 0.5*options.deformation, 0.3);
prepare(u);
for (int i = 0; i < options.directions; ++i) {
mfem::Vector a,b;
fast(i,a); reference(i,b);
const auto name = "source_displacement_state"+std::to_string(state)+"_direction"+std::to_string(i);
checks.Add(name, RelativeError(a,b), 2.0e-11);
if (options.check == "full" && i == 0) {
operators::kernels::apply_mapped_source_variation(f, *f.domainMapperStateless, rho, u, directions[i], b);
checks.Add(name+"_stateless", RelativeError(a,b), 2.0e-11);
}
}
mfem::Vector zero(directions[0]); zero=0.0;
mfem::Vector a;
source.MultDisplacementVariationTrue(rho,zero,a);
checks.Add("source_displacement_zero_state"+std::to_string(state),Norm(a),0.0);
}
prepare(base);
if (options.check == "full") {
mfem::Vector analytic; fast(0,analytic);
for (double step : {options.finiteDifferenceStep,0.5*options.finiteDifferenceStep}) {
mfem::Vector up(base),um(base),a,b;
up.Add(step,directions[0]); um.Add(-step,directions[0]);
source.PreparePrimal(source.GetDisplacementMap().gather(up)); source.Mult(reducedRho,a);
source.PreparePrimal(source.GetDisplacementMap().gather(um)); source.Mult(reducedRho,b);
a-=b; a/=2.0*step;
auto expected=source.GetPotentialMap().gather(analytic);
checks.Add("source_displacement_fd_"+std::to_string(step),RelativeError(a,expected),2.0e-6);
}
}
}
if (!checks.passed) return 3;
prepare(base);
metadata << "source_tensor_variation_elements=" << source.GetTensorVariationElementCount() << '\n';
Benchmark(options,"gravity-source-displacement","cached_pointwise_reference",fast,reference,timings,summary);
return 0;
}
int RunHydrostaticPreparation(const Options &options, fem::FEM &f, std::ofstream &metadata,
std::ofstream &summary, std::ofstream &timings, CheckReport &checks) {
operators::PreparedHydrostaticEquilibriumOperator fast(f,*f.domainMapperStateless),reference(f,*f.domainMapperStateless);
auto h = fast.GetEnthalpyMap().gather(ScalarField(*f.enthalpyFes,0.73,0.27));
auto phi = fast.GetGravityPotentialMap().gather(ScalarField(*f.gravityPotentialFes,-0.38,0.61));
const auto dh = fast.GetEnthalpyMap().gather(ScalarField(*f.enthalpyFes,0.13,0.41));
const auto dp = fast.GetGravityPotentialMap().gather(ScalarField(*f.gravityPotentialFes,0.09,0.71));
const auto du = fast.GetDisplacementMap().gather(Displacement(f,0.31,0.4));
const std::vector<mfem::Vector> states{
fast.GetDisplacementMap().gather(Displacement(f,options.deformation)),
fast.GetDisplacementMap().gather(Displacement(f,0.0)),
fast.GetDisplacementMap().gather(Displacement(f,0.5*options.deformation,0.3))};
mfem::Vector velocity(3),center(3); center=0.0;
velocity(0)=0.17; velocity(1)=-0.21; velocity(2)=0.39;
const physics::RigidRotation rotation(velocity,center);
operators::context::hydrostatic::HydrostaticEquilibriumDependencies dependencies{
.discretization={1,1},.enthalpy={2,1},.gravityPotential={3,1},
.displacement={4,1},.rotation={5,1},.bernoulliConstant={6,1}};
int state=0;
auto prepare = [&](auto &op,bool old) {
const auto view=operators::context::hydrostatic::HydrostaticEquilibriumStateView{
.enthalpy=h,.gravityPotential=phi,.displacement=states[state],.bernoulliConstant=0.39};
const auto before=op.GetResidualPreparationCount();
if (old) op.PrepareReference(view,dependencies,rotation); else op.Prepare(view,dependencies,rotation);
if (op.GetResidualPreparationCount()!=before+1)
throw std::runtime_error("Hydrostatic preparation benchmark did not rebuild.");
};
const auto start=Clock::now(); prepare(fast,false);
summary << "hydrostatic-preparation.first_optimized_seconds=" << std::chrono::duration<double>(Clock::now()-start).count() << '\n';
if (options.check!="none") {
for (state=0;state<(options.check=="full"?3:1);++state) {
++dependencies.displacement.revision;
++dependencies.enthalpy.revision; ++dependencies.gravityPotential.revision;
h*=1.01; phi*=0.99;
prepare(fast,false); prepare(reference,true);
auto compare=[&](const std::string &name,auto &&action) {
mfem::Vector a,b; action(fast,a); action(reference,b);
checks.Add("hydro_prepare_state"+std::to_string(state)+"_"+name,RelativeError(a,b),2.0e-11);
};
compare("residual",[](auto &op,auto &out){op.BuildResidual(out);});
compare("enthalpy",[&](auto &op,auto &out){op.ApplyEnthalpyJacobianAction(dh,out);});
compare("potential",[&](auto &op,auto &out){op.ApplyGravityPotentialJacobianAction(dp,out);});
compare("bernoulli",[&](auto &op,auto &out){op.ApplyBernoulliConstantJacobianAction(0.13,out);});
compare("displacement",[&](auto &op,auto &out){op.ApplyDisplacementJacobianAction(du,out);});
compare("complete",[&](auto &op,auto &out){op.ApplyCompleteJacobianAction(dh,dp,0.13,du,out);});
compare("diagonal",[](auto &op,auto &out){op.AssembleEnthalpyJacobianDiagonal(out);});
const auto view=operators::context::hydrostatic::HydrostaticEquilibriumStateView{
.enthalpy=h,.gravityPotential=phi,.displacement=states[state],.bernoulliConstant=0.39};
checks.Add("hydro_prepare_unchanged_state"+std::to_string(state),
fast.Prepare(view,dependencies,rotation).DidAnyWork()?1.0:0.0,0.0);
}
}
if (options.check == "full") {
state = 0;
++dependencies.displacement.revision;
prepare(fast, false);
mfem::Vector analytic;
fast.ApplyCompleteJacobianAction(dh, dp, 0.13, du, analytic);
for (double step : {options.finiteDifferenceStep, 0.5*options.finiteDifferenceStep}) {
mfem::Vector hp(h), hm(h), pp(phi), pm(phi), up(states[0]), um(states[0]), a, b;
hp.Add(step, dh); hm.Add(-step, dh);
pp.Add(step, dp); pm.Add(-step, dp);
up.Add(step, du); um.Add(-step, du);
++dependencies.enthalpy.revision;
++dependencies.gravityPotential.revision;
++dependencies.displacement.revision;
++dependencies.bernoulliConstant.revision;
fast.Prepare({.enthalpy=hp, .gravityPotential=pp, .displacement=up,
.bernoulliConstant=0.39+step*0.13}, dependencies, rotation);
fast.BuildResidual(a);
++dependencies.enthalpy.revision;
++dependencies.gravityPotential.revision;
++dependencies.displacement.revision;
++dependencies.bernoulliConstant.revision;
fast.Prepare({.enthalpy=hm, .gravityPotential=pm, .displacement=um,
.bernoulliConstant=0.39-step*0.13}, dependencies, rotation);
fast.BuildResidual(b);
a -= b; a /= 2.0*step;
checks.Add("hydro_prepare_complete_fd_"+std::to_string(step), RelativeError(a, analytic), 2.0e-6);
}
mfem::ParGridFunction folded(f.displacementFes.get());
mfem::VectorFunctionCoefficient fold(3, [](const mfem::Vector &x, mfem::Vector &u) {
u = x; u *= -2.0;
});
folded.ProjectCoefficient(fold);
mfem::Vector foldedTrue; folded.GetTrueDofs(foldedTrue);
const auto rejectedDisplacement = fast.GetDisplacementMap().gather(foldedTrue);
++dependencies.displacement.revision;
const auto rejectedView = operators::context::hydrostatic::HydrostaticEquilibriumStateView{
.enthalpy=h, .gravityPotential=phi, .displacement=rejectedDisplacement, .bernoulliConstant=0.39};
const auto rejected = fast.TryPrepare(rejectedView, dependencies, rotation);
const auto referenceRejected = reference.TryPrepareReference(rejectedView, dependencies, rotation);
checks.Add("hydro_prepare_fold_rejected",
!rejected && !referenceRejected &&
rejected.error().mappingStatus == mapping::MappingStatus::non_positive_determinant &&
rejected.error().mappingStatus == referenceRejected.error().mappingStatus ? 0.0 : 1.0, 0.0);
checks.Add("hydro_prepare_rejected_state_unavailable",
!fast.IsPrepared() && !reference.IsPrepared() ? 0.0 : 1.0, 0.0);
}
if(!checks.passed)return 3;
state=0;
++dependencies.displacement.revision;
++dependencies.enthalpy.revision; ++dependencies.gravityPotential.revision;
++dependencies.bernoulliConstant.revision;
prepare(fast,false); prepare(reference,true);
metadata << "hydro_preparation_tensor_elements=" << fast.GetTensorDisplacementElementCount() << '\n';
BenchmarkPreparation(options,"hydrostatic-preparation",
[&]{state=(state+1)%static_cast<int>(states.size());++dependencies.displacement.revision;},
[&]{prepare(fast,false);},[&]{prepare(reference,true);},
[&](auto &out){fast.BuildResidual(out);},[&](auto &out){reference.BuildResidual(out);},timings,summary);
return 0;
}
int RunGravitySchurPreparation(const Options &options,fem::FEM &f,std::ofstream &metadata,
std::ofstream &summary,std::ofstream &timings,CheckReport &checks) {
preconditioning::PreparedGravityDivergenceCache cache;
operators::PreparedMappedHDivMassOperator mass(f,*f.domainMapperStateless);
std::vector<mfem::Vector> diagonals(3);
for(int i=0;i<3;++i) {
mass.PreparePrimal(mass.GetDisplacementMap().gather(Displacement(f,0.5*i*options.deformation,0.2)));
mass.AssembleTrueDiagonal(diagonals[i]);
}
const auto x=Flux(f.gravityPotentialFes->GetTrueVSize(),0.37);
std::unique_ptr<mfem::HypreParMatrix> fast,reference;
int state=0;
auto prepareFast=[&]{fast=preconditioning::assembleGravityPotentialSchurSurrogate(f,diagonals[state],&cache);};
auto prepareReference=[&]{reference=preconditioning::assembleGravityPotentialSchurSurrogate(f,diagonals[state]);};
const auto start=Clock::now(); prepareFast();
summary << "gravity-schur-preparation.cold_cache_seconds=" << std::chrono::duration<double>(Clock::now()-start).count() << '\n';
auto compare=[&](const std::string &label){
mfem::Vector a(fast->Height()),b(reference->Height());
fast->Mult(x,a);reference->Mult(x,b);checks.Add(label,RelativeError(a,b),2.0e-11);
};
if(options.check!="none") {
for(state=0;state<3;++state){prepareFast();prepareReference();compare("schur_state"+std::to_string(state));}
// An exact scaling check also works with --deformation 0, and catches
// both a stale Schur matrix and accidental cumulative scaling of B.
const auto originalDiagonal = diagonals[0];
state = 0; prepareFast();
mfem::Vector originalAction(fast->Height()), changedAction(fast->Height());
fast->Mult(x, originalAction);
diagonals[0] *= 2.0; prepareFast(); prepareReference(); compare("schur_changed_mass_diagonal");
fast->Mult(x, changedAction);
originalAction *= 0.5;
checks.Add("schur_inverse_mass_scaling", RelativeError(changedAction, originalAction), 2.0e-11);
diagonals[0] = originalDiagonal; prepareFast(); prepareReference(); compare("schur_restored_mass_diagonal");
checks.Add("divergence_shared_across_geometry",cache.GetBuildCount()==1?0.0:1.0,0.0);
state=0;
if(options.check=="full") {
cache.Clear();prepareFast();prepareReference();compare("schur_after_clear");
checks.Add("divergence_clear_rebuild",cache.GetBuildCount()==2?0.0:1.0,0.0);
f.mesh->NodesUpdated();prepareFast();prepareReference();compare("schur_after_nodes_update");
checks.Add("divergence_nodes_update_rebuild",cache.GetBuildCount()==3?0.0:1.0,0.0);
}
}
if(!checks.passed)return 3;
state=0;const auto builds=cache.GetBuildCount();
BenchmarkPreparation(options,"gravity-schur-preparation",
[&]{state=(state+1)%3;},prepareFast,prepareReference,
[&](auto &out){out.SetSize(fast->Height());fast->Mult(x,out);},
[&](auto &out){out.SetSize(reference->Height());reference->Mult(x,out);},timings,summary);
if(options.check!="none")checks.Add("divergence_warm_rebuilds",double(cache.GetBuildCount()-builds),0.0);
metadata << "fixed_divergence_assemblies=" << cache.GetBuildCount() << '\n';
return checks.passed?0:3;
}
int Run(const Options &options) {
if (options.help) {
std::cout << "system_isolation_profile_harness\n"
" --system hdiv-displacement|hdiv-mass|mass-normalization|angular-momentum|\n"
" hydrostatic-displacement|stellar-derivatives|barotropic-closure|barotropic-enthalpy|\n"
" barotropic-displacement|pressure-displacement|pressure-enthalpy|gravity-force|rotation-force|\n"
" linear-derivatives|preparation-optimizations|runtime-optimizations|\n"
" geometry-preparation|normalization-diagonals|application-optimizations|\n"
" final-optimizations|gravity-source-displacement|hydrostatic-preparation|gravity-schur-preparation|\n"
" geometry-inspection|gravity-source-preparation|all (default hdiv-displacement)\n"
" --implementation both|optimized|reference (default both)\n"
" --check full|reference|none (default full)\n"
" --mesh sandbox.smesh --output NEW_DIRECTORY --refinements 0\n"
" --samples 5 --iterations 3 --warmup 1 --directions 3\n"
" --prepare-samples 3 --prepare-iterations 1 (prepare-samples 0 skips preparation timing)\n"
" --deformation 0.7 --fd-step 1e-5 --polytropic-index 1\n"
"Single MPI rank, CPU. No Newton solve. Exit 0: requested checks pass;\n"
"2: execution error; 3: check failure. Existing output directories are refused.\n"
"hdiv-displacement compares the optimized and retained prepared derivatives.\n"
"hdiv-mass compares the new CPU application with retained MFEM partial assembly.\n"
"Full checks also compare the independent stateless mass kernel.\n";
return 0;
}
if (!std::filesystem::create_directory(options.output))
throw std::invalid_argument("Output directory already exists: " + options.output.string());
auto metadata = File(options.output / "metadata.txt");
auto summary = File(options.output / "summary.txt");
auto timings = File(options.output / "timings.csv");
CheckReport checks{File(options.output / "checks.csv")};
checks.stream << "check,normalized_error,tolerance,passed\n";
timings << "system,implementation,sample,applications,total_seconds,seconds_per_application,last_output_norm\n";
metadata << "mesh=" << std::filesystem::absolute(options.mesh).string()
<< "\nsystem=" << options.system << "\nimplementation=" << options.implementation
<< "\ncheck=" << options.check << "\nrefinements=" << options.refinements
<< "\nsamples=" << options.samples << "\niterations=" << options.iterations
<< "\nwarmup=" << options.warmup << "\ndirections=" << options.directions
<< "\nprepare_samples=" << options.prepareSamples << "\nprepare_iterations=" << options.prepareIterations
<< "\npreparation_benchmark=forced_revisions_shared_gravity_excluded"
<< "\npreparation_reference=original_matrix_assembly_and_separate_force_residual_traversals"
<< "\nhydrostatic_reference=previous_cached_pointwise"
<< "\nhdiv_displacement_reference=previous_optimized_pointwise"
<< "\ngravity_force_reference=previous_optimized_pointwise"
<< "\ngeometry_preparation_reference=same_operator_with_fixed_mesh_cache_disabled"
<< "\nhdiv_mass_reference=previous_mfem_partial_assembly"
<< "\nclosure_reference=previous_optimized_pointwise"
<< "\nrotation_reference=previous_optimized_pointwise"
<< "\ngeometry_inspection_reference=original_gridfunction_gradient_same_points"
<< "\ngravity_source_preparation_reference=previous_mesh_weight_recomputation"
<< "\npressure_displacement_reference=previous_optimized_pointwise"
<< "\npressure_preparation_reference=previous_blocked_matrix_assembly_and_pointwise_fields"
<< "\ngravity_source_displacement_reference=previous_pointwise"
<< "\nhydrostatic_preparation_reference=previous_pointwise_fields_and_dense_matrix_assembly"
<< "\ngravity_schur_preparation_reference=fixed_divergence_reassembled_every_call"
<< "\ndeformation=" << options.deformation << "\nfd_step=" << options.finiteDifferenceStep
<< "\nmpi_ranks=1\ndevice=cpu\ncompiler=" << __VERSION__
<< "\nuniform_order_increment=" << MEAN_FIELD_UNIFORM_POLYNOMIAL_ORDER_INCREMENT << '\n';
metadata.flush();
utils::Args arguments;
arguments.mesh_file = options.mesh;
const auto setupStart = Clock::now();
auto f = fem::setup_fem(arguments.mesh_file, arguments, options.refinements);
if (!f.okay() || f.mesh->Dimension() != 3) throw std::runtime_error("A valid 3D discretization is required.");
summary << "fem_setup_seconds=" << std::chrono::duration<double>(Clock::now() - setupStart).count() << '\n';
int vacuumElements = 0;
for (int element = 0; element < f.mesh->GetNE(); ++element)
vacuumElements += f.domainMapperStateless->IsCompactifiedElement(*f.mesh->GetElementTransformation(element));
metadata << "elements=" << f.mesh->GetNE() << "\nstellar_elements=" << f.mesh->GetNE() - vacuumElements
<< "\nvacuum_elements=" << vacuumElements << "\nflux_true_dofs=" << f.gravityFluxFes->GetTrueVSize()
<< "\ndisplacement_true_dofs=" << f.displacementFes->GetTrueVSize()
<< "\ndensity_true_dofs=" << f.densityFes->GetTrueVSize()
<< "\nenthalpy_true_dofs=" << f.enthalpyFes->GetTrueVSize() << '\n';
const bool finalGroup = options.system == "final-optimizations";
if (finalGroup || options.system == "gravity-source-displacement" ||
options.system == "hydrostatic-preparation" || options.system == "gravity-schur-preparation") {
int result=0;
if(finalGroup) {
auto pressureOptions=options;pressureOptions.system="pressure-displacement";
result=RunLinearDerivatives(pressureOptions,f,metadata,summary,timings,checks);
if(result)return result;
}
if(finalGroup || options.system=="gravity-source-displacement") {
result=RunSourceDisplacement(options,f,metadata,summary,timings,checks);
if(result)return result;
}
if(finalGroup || options.system=="hydrostatic-preparation") {
result=RunHydrostaticPreparation(options,f,metadata,summary,timings,checks);
if(result)return result;
}
if(finalGroup || options.system=="gravity-schur-preparation") {
result=RunGravitySchurPreparation(options,f,metadata,summary,timings,checks);
if(result)return result;
}
summary << "final.checks=" << (options.check=="none"?"not_run":"passed") << '\n';
std::cout << "Results: " << std::filesystem::absolute(options.output) << std::endl;
return 0;
}
if (options.system == "all") {
int result = RunSourceDisplacement(options, f, metadata, summary, timings, checks);
if (result != 0) return result;
result = RunHydrostaticPreparation(options, f, metadata, summary, timings, checks);
if (result != 0) return result;
result = RunGravitySchurPreparation(options, f, metadata, summary, timings, checks);
if (result != 0) return result;
}
const bool runtimeGroup = options.system == "runtime-optimizations";
const bool applicationGroup = options.system == "application-optimizations";
if (applicationGroup || options.system == "all" || options.system == "geometry-inspection") {
const int result = RunGeometryInspection(options, f, metadata, summary, timings, checks);
if (result != 0) return result;
if (options.system == "geometry-inspection") {
std::cout << "Results: " << std::filesystem::absolute(options.output) << std::endl;
return 0;
}
}
if (applicationGroup || options.system == "all" || options.system == "gravity-source-preparation") {
const int result = RunGravitySourcePreparation(options, f, summary, timings, checks);
if (result != 0) return result;
if (options.system == "gravity-source-preparation") {
std::cout << "Results: " << std::filesystem::absolute(options.output) << std::endl;
return 0;
}
}
if (runtimeGroup || options.system == "all" || options.system == "normalization-diagonals") {
const int result = RunNormalizationDiagonals(options, f, summary, timings, checks);
if (result != 0) return result;
if (options.system == "normalization-diagonals") {
summary << "checks=" << (options.check == "none" ? "not_run" : "passed") << '\n';
std::cout << "Results: " << std::filesystem::absolute(options.output) << std::endl;
return 0;
}
}
if (runtimeGroup || options.system == "all" || options.system == "geometry-preparation") {
const int result = RunGeometryPreparation(options, f, metadata, summary, timings, checks);
if (result != 0) return result;
if (options.system == "geometry-preparation") {
summary << "checks=" << (options.check == "none" ? "not_run" : "passed") << '\n';
std::cout << "Results: " << std::filesystem::absolute(options.output) << std::endl;
return 0;
}
}
if (IsLinearDerivativeSystem(options.system) && !applicationGroup) {
int result = RunLinearDerivatives(options, f, metadata, summary, timings, checks);
if (result == 0 && options.system == "preparation-optimizations") {
auto hydroOptions = options;
hydroOptions.system = "hydrostatic-displacement";
result = RunStellarDerivatives(hydroOptions, f, metadata, summary, timings, checks);
}
std::cout << "Results: " << std::filesystem::absolute(options.output) << std::endl;
return result;
}
if (!runtimeGroup && !applicationGroup && options.system != "hdiv-displacement" && options.system != "hdiv-mass" && options.system != "all") {
const int result = RunStellarDerivatives(options, f, metadata, summary, timings, checks);
std::cout << "Results: " << std::filesystem::absolute(options.output) << std::endl;
return result;
}
const auto constructStart = Clock::now();
operators::PreparedMappedHDivMassOperator mass(f, *f.domainMapperStateless);
summary << "operator_construction_seconds=" << std::chrono::duration<double>(Clock::now() - constructStart).count() << '\n';
const mfem::Vector base = Displacement(f, options.deformation);
std::vector<mfem::Vector> fluxes, directions, reducedFluxes;
for (int index = 0; index < options.directions; ++index) {
fluxes.push_back(Flux(f.gravityFluxFes->GetTrueVSize(), 0.31 + 0.47 * index));
directions.push_back(Displacement(f, 0.43, 0.2 + 0.7 * index));
if (index == 1) directions.back().Add(0.002, Flux(f.displacementFes->GetTrueVSize(), 0.37));
reducedFluxes.push_back(mass.GetFluxMap().gather(fluxes.back()));
}
const bool derivative = options.system != "hdiv-mass" && !applicationGroup;
auto prepare = [&](const mfem::Vector &state) {
const auto reduced = mass.GetDisplacementMap().gather(state);
if (derivative) mass.Prepare(reduced); else mass.PreparePrimal(reduced);
};
auto start = Clock::now();
prepare(base);
summary << "initial_prepare_seconds=" << std::chrono::duration<double>(Clock::now() - start).count() << '\n';
metadata << "fast_mass_elements=" << mass.GetFastMassElementCount() << '\n';
metadata << "fast_variation_elements=" << mass.GetFastVariationElementCount()
<< "\ntensor_variation_elements=" << mass.GetTensorVariationElementCount() << '\n';
std::size_t quadraturePoints = 0;
mass.VisitMappedGeometryRules([&](int, const mfem::IntegrationRule &rule) { quadraturePoints += rule.GetNPoints(); });
metadata << "variation_quadrature_points=" << quadraturePoints << '\n';
metadata.flush();
const Apply fastDerivative = [&](int index, mfem::Vector &out) {
mass.MultDisplacementVariationTrue(fluxes[index], directions[index], out);
};
const Apply referenceDerivative = [&](int index, mfem::Vector &out) {
mass.MultDisplacementVariationTrueCachedReference(fluxes[index], directions[index], out);
};
const Apply fastMass = [&](int index, mfem::Vector &out) { mass.Mult(reducedFluxes[index], out); };
const Apply referenceMass = [&](int index, mfem::Vector &out) { mass.MultReference(reducedFluxes[index], out); };
if (options.check != "none") {
const auto validationStart = Clock::now();
// Reprepare distinct geometries to expose stale geometry caches.
const std::vector<double> scales = options.check == "full" ? std::vector<double>{0.0, options.deformation}
: std::vector<double>{options.deformation};
for (std::size_t state = 0; state < scales.size(); ++state) {
const auto displacement = Displacement(f, scales[state]);
prepare(displacement);
for (int index = 0; index < options.directions; ++index) {
mfem::Vector actual, expected;
const std::string label = "state" + std::to_string(state) + "_direction" + std::to_string(index);
if (derivative) {
fastDerivative(index, actual);
referenceDerivative(index, expected);
checks.Add(label + "_cached_derivative", RelativeError(actual, expected), 2.0e-11);
mass.MultDisplacementVariationTrueReference(fluxes[index], directions[index], expected);
checks.Add(label + "_general_prepared_derivative", RelativeError(actual, expected), 2.0e-11);
if (options.check == "full" && index == 0) {
operators::kernels::apply_mapped_hdiv_mass_variation(
f, *f.domainMapperStateless, fluxes[index], displacement, directions[index], expected);
checks.Add(label + "_stateless_derivative", RelativeError(actual, expected), 2.0e-11);
}
}
if (options.system != "hdiv-displacement" && !runtimeGroup) {
fastMass(index, actual);
referenceMass(index, expected);
checks.Add(label + "_mfem_pa_mass", RelativeError(actual, expected), 2.0e-11);
if (options.check == "full" && index == 0) {
mfem::Vector full;
operators::kernels::apply_mapped_hdiv_mass(f, *f.domainMapperStateless, fluxes[index], displacement, full);
expected = mass.GetFluxMap().gather(full);
checks.Add(label + "_stateless_mass", RelativeError(actual, expected), 2.0e-11);
}
}
}
if (options.system != "hdiv-displacement" && !runtimeGroup) {
const std::string label = "state" + std::to_string(state) + "_mass";
mfem::Vector other = mass.GetFluxMap().gather(
Flux(f.gravityFluxFes->GetTrueVSize(), 1.37));
mfem::Vector firstAction, otherAction, combinedAction, expected;
mass.Mult(reducedFluxes[0], firstAction);
mass.Mult(other, otherAction);
const double xy = reducedFluxes[0] * otherAction;
const double yx = other * firstAction;
checks.Add(label + "_symmetry", std::abs(xy-yx) / std::max({1.0, std::abs(xy), std::abs(yx)}), 2.0e-11);
checks.Add(label + "_positive_energy", reducedFluxes[0] * firstAction > 0.0 ? 0.0 : 1.0, 0.0);
mfem::Vector combined(reducedFluxes[0]);
combined *= 0.7;
combined.Add(-0.4, other);
mass.Mult(combined, combinedAction);
expected = firstAction;
expected *= 0.7;
expected.Add(-0.4, otherAction);
checks.Add(label + "_linearity", RelativeError(combinedAction, expected), 2.0e-11);
combined = 0.0;
mass.Mult(combined, combinedAction);
checks.Add(label + "_zero", combinedAction.Norml2(), 0.0);
mass.MultReference(reducedFluxes[0], expected);
mass.Mult(reducedFluxes[0], combinedAction);
checks.Add(label + "_reference_restores_fast_path", RelativeError(combinedAction, firstAction), 0.0);
}
}
if (derivative && options.check == "full") {
mfem::Vector analyticTrue;
fastDerivative(0, analyticTrue);
const auto analytic = mass.GetFluxMap().gather(analyticTrue);
operators::PreparedMappedHDivMassOperator differenceMass(f, *f.domainMapperStateless);
for (const double step : {options.finiteDifferenceStep, 0.5 * options.finiteDifferenceStep}) {
mfem::Vector plus(base), minus(base), plusAction, minusAction;
plus.Add(step, directions[0]);
minus.Add(-step, directions[0]);
differenceMass.PreparePrimal(differenceMass.GetDisplacementMap().gather(plus));
differenceMass.Mult(reducedFluxes[0], plusAction);
differenceMass.PreparePrimal(differenceMass.GetDisplacementMap().gather(minus));
differenceMass.Mult(reducedFluxes[0], minusAction);
plusAction -= minusAction;
plusAction /= 2.0 * step;
checks.Add("finite_difference_" + std::to_string(step), RelativeError(analytic, plusAction), 2.0e-6);
}
}
summary << "validation_seconds=" << std::chrono::duration<double>(Clock::now() - validationStart).count() << '\n';
if (!checks.passed) {
summary << "checks_passed=0\nbenchmark_skipped=1\n";
return 3;
}
}
start = Clock::now();
prepare(base);
summary << "repeat_prepare_seconds=" << std::chrono::duration<double>(Clock::now() - start).count() << '\n';
summary << "checks=" << (options.check == "none" ? "not_run" : "passed") << '\n';
summary.flush();
if (derivative) Benchmark(options, "hdiv-displacement", "cached_pointwise_reference", fastDerivative, referenceDerivative, timings, summary);
if (options.system != "hdiv-displacement" && !runtimeGroup) Benchmark(options, "hdiv-mass", "mfem_pa_reference", fastMass, referenceMass, timings, summary);
if (applicationGroup) {
auto linearOptions = options;
linearOptions.prepareSamples = 0; // preparation measurements above isolate this batch
const int result = RunLinearDerivatives(linearOptions, f, metadata, summary, timings, checks);
if (result != 0) return result;
}
if (runtimeGroup) {
auto gravityOptions = options;
gravityOptions.system = "gravity-force";
gravityOptions.prepareSamples = 0; // this group measures new application and geometry changes
const int result = RunLinearDerivatives(gravityOptions, f, metadata, summary, timings, checks);
if (result != 0) return result;
}
if (options.system == "all") {
const int result = RunStellarDerivatives(options, f, metadata, summary, timings, checks);
if (result != 0) return result;
const int linearResult = RunLinearDerivatives(options, f, metadata, summary, timings, checks);
if (linearResult != 0) return linearResult;
}
std::cout << "Results: " << std::filesystem::absolute(options.output) << std::endl;
return 0;
}
}
int main(int argc, char **argv) {
// Parse before MPI initialization so --help does not require a runtime.
Options options;
try {
options = Parse(argc, argv);
if (options.help) return Run(options);
} catch (const std::exception &error) {
std::cerr << error.what() << '\n';
return 2;
}
mfem::Mpi::Init(argc, argv);
int result = 0;
try {
int ranks = 0;
MPI_Comm_size(MPI_COMM_WORLD, &ranks);
if (ranks != 1) throw std::invalid_argument("This harness currently requires exactly one MPI rank.");
mfem::Device device("cpu");
std::cout << std::setprecision(12);
result = Run(options);
} catch (const std::exception &error) {
std::cerr << "Isolation harness error: " << error.what() << std::endl;
result = 2;
}
mfem::Mpi::Finalize();
return result;
}
Reference in New Issue View Git Blame Copy Permalink