782 lines
36 KiB
C++
782 lines
36 KiB
C++
#include <catch2/catch_test_macros.hpp>
|
|
|
|
#include <algorithm>
|
|
#include <array>
|
|
#include <cmath>
|
|
#include <limits>
|
|
#include <map>
|
|
#include <string>
|
|
#include <string_view>
|
|
#include <utility>
|
|
#include <vector>
|
|
|
|
#include <mfem.hpp>
|
|
#include <mpi.h>
|
|
|
|
import experiment;
|
|
import experiment.stellar_null_space;
|
|
import mean_field;
|
|
import test_helpers;
|
|
|
|
namespace {
|
|
namespace null_space = experiment::null_space;
|
|
|
|
struct GaugeMode final {
|
|
std::string name;
|
|
std::string family;
|
|
int axis{-1};
|
|
bool requiresGravityCompletion{true};
|
|
mfem::Vector direction;
|
|
};
|
|
|
|
class GravityUnknownJacobian final : public mfem::Operator {
|
|
public:
|
|
explicit GravityUnknownJacobian(
|
|
const mean_field::operators::PreparedStellarEquilibriumOperator &stellarOperator
|
|
)
|
|
: mfem::Operator(
|
|
stellarOperator.GetLayout().size(null_space::gravityGradientValue) +
|
|
stellarOperator.GetLayout().size(null_space::gravityPotentialValue)
|
|
),
|
|
m_stellarOperator(stellarOperator),
|
|
m_gravityGradientSize(stellarOperator.GetLayout().size(null_space::gravityGradientValue)) {
|
|
MFEM_VERIFY(Width() == Height(), "The restricted gravity Jacobian must be square.");
|
|
}
|
|
|
|
void Mult(
|
|
const mfem::Vector &gravityDirection,
|
|
mfem::Vector &gravityAction
|
|
) const override {
|
|
MFEM_VERIFY(gravityDirection.Size() == Width(), "The restricted gravity direction has the wrong size.");
|
|
|
|
const mfem::Vector gravityGradientDirection(
|
|
const_cast<mfem::real_t *>(gravityDirection.GetData()), m_gravityGradientSize
|
|
);
|
|
const mfem::Vector gravityPotentialDirection(
|
|
const_cast<mfem::real_t *>(gravityDirection.GetData()) + m_gravityGradientSize,
|
|
Width() - m_gravityGradientSize
|
|
);
|
|
|
|
m_stellarOperator.GetGravityOperator().ApplyGravityUnknowns(
|
|
gravityGradientDirection, gravityPotentialDirection,
|
|
m_stellarOperator.GetGravityContext().GetGeometryContext(), gravityAction
|
|
);
|
|
}
|
|
|
|
[[nodiscard]] int gravity_gradient_size() const noexcept {
|
|
return m_gravityGradientSize;
|
|
}
|
|
|
|
private:
|
|
const mean_field::operators::PreparedStellarEquilibriumOperator &m_stellarOperator;
|
|
int m_gravityGradientSize;
|
|
};
|
|
|
|
struct GravityCompletionResult final {
|
|
mfem::Vector direction;
|
|
double rightHandSideNorm{0.0};
|
|
double residualNorm{0.0};
|
|
double relativeResidual{0.0};
|
|
double finalNorm{0.0};
|
|
int iterations{0};
|
|
bool solvePerformed{false};
|
|
};
|
|
|
|
void add_block_metrics(
|
|
std::map<
|
|
std::string,
|
|
double> &metrics,
|
|
const std::string &prefix,
|
|
const std::array<
|
|
double,
|
|
6> &norms
|
|
) {
|
|
for (std::size_t block = 0; block < norms.size(); ++block) {
|
|
metrics.emplace(prefix + null_space::residualBlockNames[block] + "_norm", norms[block]);
|
|
}
|
|
}
|
|
|
|
[[nodiscard]] mfem::Vector gravity_residual_blocks(
|
|
const mfem::Vector &completeAction,
|
|
const mean_field::operators::StellarEquilibriumLayout &layout
|
|
) {
|
|
const mfem::Vector gradient =
|
|
null_space::const_residual_view(completeAction, layout, null_space::gravityGradientResidual);
|
|
const mfem::Vector potential =
|
|
null_space::const_residual_view(completeAction, layout, null_space::gravityPotentialResidual);
|
|
|
|
mfem::Vector result(gradient.Size() + potential.Size());
|
|
mfem::Vector(result.GetData(), gradient.Size()) = gradient;
|
|
mfem::Vector(result.GetData() + gradient.Size(), potential.Size()) = potential;
|
|
return result;
|
|
}
|
|
|
|
void assign_gravity_completion(
|
|
mfem::Vector &completeDirection,
|
|
const mean_field::operators::StellarEquilibriumLayout &layout,
|
|
const mfem::Vector &gravityCompletion,
|
|
const int gravityGradientSize
|
|
) {
|
|
const mfem::Vector gravityGradient(
|
|
const_cast<mfem::real_t *>(gravityCompletion.GetData()), gravityGradientSize
|
|
);
|
|
const mfem::Vector gravityPotential(
|
|
const_cast<mfem::real_t *>(gravityCompletion.GetData()) + gravityGradientSize,
|
|
gravityCompletion.Size() - gravityGradientSize
|
|
);
|
|
null_space::assign_value_block(completeDirection, layout, null_space::gravityGradientValue, gravityGradient);
|
|
null_space::assign_value_block(completeDirection, layout, null_space::gravityPotentialValue, gravityPotential);
|
|
}
|
|
|
|
[[nodiscard]] GravityCompletionResult solve_gravity_completion(
|
|
const mfem::Vector &prescribedAction,
|
|
const mean_field::operators::StellarEquilibriumLayout &layout,
|
|
const MPI_Comm communicator,
|
|
GravityUnknownJacobian &gravityJacobian,
|
|
mfem::MINRESSolver &gravitySolver
|
|
) {
|
|
mfem::Vector rightHandSide = gravity_residual_blocks(prescribedAction, layout);
|
|
rightHandSide *= -1.0;
|
|
|
|
GravityCompletionResult result;
|
|
result.direction.SetSize(gravityJacobian.Width());
|
|
result.direction = 0.0;
|
|
result.rightHandSideNorm = null_space::global_norm(rightHandSide, communicator);
|
|
|
|
const double skipThreshold = 100.0 * std::numeric_limits<double>::epsilon();
|
|
if (result.rightHandSideNorm <= skipThreshold) {
|
|
return result;
|
|
}
|
|
|
|
gravitySolver.Mult(rightHandSide, result.direction);
|
|
REQUIRE(gravitySolver.GetConverged());
|
|
|
|
mfem::Vector action;
|
|
gravityJacobian.Mult(result.direction, action);
|
|
action -= rightHandSide;
|
|
|
|
result.residualNorm = null_space::global_norm(action, communicator);
|
|
result.relativeResidual = result.residualNorm / result.rightHandSideNorm;
|
|
result.finalNorm = gravitySolver.GetFinalNorm();
|
|
result.iterations = gravitySolver.GetNumIterations();
|
|
result.solvePerformed = true;
|
|
|
|
REQUIRE(std::isfinite(result.relativeResidual));
|
|
return result;
|
|
}
|
|
|
|
class ExtensionAwareHomologyScalarCoefficient final : public mfem::Coefficient {
|
|
public:
|
|
ExtensionAwareHomologyScalarCoefficient(
|
|
const mfem::ParGridFunction &baseField,
|
|
const mfem::ParGridFunction &coordinateVelocity,
|
|
const mfem::Vector &referenceCenter,
|
|
const double physicalScalingExponent
|
|
)
|
|
: m_baseField(&baseField),
|
|
m_coordinateVelocity(&coordinateVelocity),
|
|
m_referenceCenter(&referenceCenter),
|
|
m_physicalScalingExponent(physicalScalingExponent) {
|
|
}
|
|
|
|
double Eval(
|
|
mfem::ElementTransformation &transformation,
|
|
const mfem::IntegrationPoint &integrationPoint
|
|
) override {
|
|
transformation.SetIntPoint(&integrationPoint);
|
|
|
|
mfem::Vector referencePosition;
|
|
mfem::Vector coordinateVelocity;
|
|
mfem::Vector baseGradient;
|
|
transformation.Transform(integrationPoint, referencePosition);
|
|
m_coordinateVelocity->GetVectorValue(transformation, integrationPoint, coordinateVelocity);
|
|
m_baseField->GetGradient(transformation, baseGradient);
|
|
|
|
coordinateVelocity -= referencePosition;
|
|
coordinateVelocity += *m_referenceCenter;
|
|
|
|
return -m_physicalScalingExponent * m_baseField->GetValue(transformation, integrationPoint) +
|
|
baseGradient * coordinateVelocity;
|
|
}
|
|
|
|
private:
|
|
const mfem::ParGridFunction *m_baseField;
|
|
const mfem::ParGridFunction *m_coordinateVelocity;
|
|
const mfem::Vector *m_referenceCenter;
|
|
double m_physicalScalingExponent;
|
|
};
|
|
|
|
[[nodiscard]] mfem::Vector project_extension_aware_homology_scalar(
|
|
mfem::ParFiniteElementSpace &finiteElementSpace,
|
|
const mean_field::field::FieldDofMap &fieldMap,
|
|
const mfem::Vector &baseReducedField,
|
|
const mfem::ParGridFunction &coordinateVelocity,
|
|
const mfem::Vector &referenceCenter,
|
|
const double physicalScalingExponent
|
|
) {
|
|
mfem::ParGridFunction baseField(&finiteElementSpace);
|
|
baseField.SetFromTrueDofs(fieldMap.scatter(baseReducedField));
|
|
|
|
ExtensionAwareHomologyScalarCoefficient coefficient(
|
|
baseField, coordinateVelocity, referenceCenter, physicalScalingExponent
|
|
);
|
|
mfem::ParGridFunction directionField(&finiteElementSpace);
|
|
directionField.ProjectCoefficient(coefficient);
|
|
|
|
mfem::Vector directionTrue;
|
|
directionField.GetTrueDofs(directionTrue);
|
|
return fieldMap.gather(directionTrue);
|
|
}
|
|
|
|
struct HomologyMassCancellation final {
|
|
double currentMass{0.0};
|
|
double targetMass{0.0};
|
|
double densityContribution{0.0};
|
|
double geometryContribution{0.0};
|
|
double completeDerivative{0.0};
|
|
};
|
|
|
|
[[nodiscard]] HomologyMassCancellation measure_homology_mass_cancellation(
|
|
const mean_field::operators::PreparedMassNormalizationOperator &massOperator,
|
|
const mfem::Vector &densityDirection,
|
|
const mfem::Vector &volumeDirection
|
|
) {
|
|
mfem::Vector densityAction;
|
|
mfem::Vector geometryAction;
|
|
mfem::Vector completeAction;
|
|
massOperator.ApplyDensityJacobianAction(densityDirection, densityAction);
|
|
massOperator.ApplyDisplacementJacobianAction(volumeDirection, geometryAction);
|
|
massOperator.ApplyCompleteJacobianAction(densityDirection, volumeDirection, completeAction);
|
|
|
|
REQUIRE(densityAction.Size() == 1);
|
|
REQUIRE(geometryAction.Size() == 1);
|
|
REQUIRE(completeAction.Size() == 1);
|
|
|
|
const double recomposedDerivative = densityAction(0) + geometryAction(0);
|
|
const double comparisonScale = std::max({1.0, std::abs(recomposedDerivative), std::abs(completeAction(0))});
|
|
CHECK(
|
|
std::abs(completeAction(0) - recomposedDerivative) <=
|
|
64.0 * std::numeric_limits<double>::epsilon() * comparisonScale
|
|
);
|
|
|
|
return {
|
|
.currentMass = massOperator.GetCurrentMass(),
|
|
.targetMass = massOperator.GetTargetMass(),
|
|
.densityContribution = densityAction(0),
|
|
.geometryContribution = geometryAction(0),
|
|
.completeDerivative = completeAction(0)
|
|
};
|
|
}
|
|
|
|
void add_homology_mass_metrics(
|
|
std::map<
|
|
std::string,
|
|
double> &metrics,
|
|
const HomologyMassCancellation &cancellation
|
|
) {
|
|
const double uncancelledMagnitude =
|
|
std::abs(cancellation.densityContribution) + std::abs(cancellation.geometryContribution);
|
|
const double targetScale = std::max(std::abs(cancellation.targetMass), std::numeric_limits<double>::epsilon());
|
|
|
|
metrics.emplace("current_mass", cancellation.currentMass);
|
|
metrics.emplace("target_mass", cancellation.targetMass);
|
|
metrics.emplace("base_mass_residual", cancellation.currentMass - cancellation.targetMass);
|
|
metrics.emplace(
|
|
"relative_base_mass_residual", (cancellation.currentMass - cancellation.targetMass) / targetScale
|
|
);
|
|
metrics.emplace("density_mass_derivative", cancellation.densityContribution);
|
|
metrics.emplace("geometry_mass_derivative", cancellation.geometryContribution);
|
|
metrics.emplace("complete_mass_derivative", cancellation.completeDerivative);
|
|
metrics.emplace("mass_derivative_uncancelled_magnitude", uncancelledMagnitude);
|
|
metrics.emplace(
|
|
"mass_derivative_relative_cancellation_error",
|
|
std::abs(cancellation.completeDerivative) /
|
|
std::max(uncancelledMagnitude, std::numeric_limits<double>::epsilon())
|
|
);
|
|
metrics.emplace("complete_mass_derivative_per_target_mass", cancellation.completeDerivative / targetScale);
|
|
}
|
|
|
|
[[nodiscard]] GaugeMode make_homology_mode(
|
|
null_space::N3Equilibrium &fixture,
|
|
const null_space::SurfaceMode &uniformRadialMode
|
|
) {
|
|
const auto &layout = fixture.stellar_operator().GetLayout();
|
|
const auto &state = fixture.state();
|
|
|
|
mfem::Vector direction(layout.value_offsets().Last());
|
|
direction = 0.0;
|
|
|
|
const mfem::Vector volumeDirection = fixture.lifted_surface_direction(uniformRadialMode.direction);
|
|
mfem::ParGridFunction coordinateVelocity(fixture.fem().displacementFes.get());
|
|
coordinateVelocity.SetFromTrueDofs(volumeDirection);
|
|
|
|
const mean_field::field::FieldDofMap densityMap =
|
|
mean_field::field::make_field_dof_map<mean_field::field::Density, null_space::DomainSchema>(
|
|
*fixture.fem().densityFes
|
|
);
|
|
const mean_field::field::FieldDofMap enthalpyMap =
|
|
mean_field::field::make_field_dof_map<mean_field::field::Enthalpy, null_space::DomainSchema>(
|
|
*fixture.fem().enthalpyFes
|
|
);
|
|
const mfem::Vector &referenceCenter = fixture.model().surfaceDeformationPrescription().referenceCenter();
|
|
|
|
const mfem::Vector densityDirection = project_extension_aware_homology_scalar(
|
|
*fixture.fem().densityFes, densityMap,
|
|
null_space::const_value_view(state, layout, null_space::densityValue), coordinateVelocity, referenceCenter,
|
|
3.0
|
|
);
|
|
null_space::assign_value_block(direction, layout, null_space::densityValue, densityDirection);
|
|
|
|
null_space::assign_value_block(
|
|
direction, layout, null_space::surfaceDeformationValue,
|
|
null_space::const_value_view(uniformRadialMode.direction, layout, null_space::surfaceDeformationValue)
|
|
);
|
|
|
|
/*
|
|
* A physical homology scales rho and h, while the power-law mesh
|
|
* extension moves interior coordinates non-affinely. The scalar
|
|
* tangents therefore contain the coordinate-composition term
|
|
* grad(f) dot (v - (X-Xc)) in addition to their physical scaling.
|
|
* Gravity is completed through the discrete mixed subsystem below,
|
|
* which also supplies the correct fixed-infinity exterior response.
|
|
*/
|
|
const mfem::Vector enthalpyDirection = project_extension_aware_homology_scalar(
|
|
*fixture.fem().enthalpyFes, enthalpyMap,
|
|
null_space::const_value_view(state, layout, null_space::enthalpyValue), coordinateVelocity, referenceCenter,
|
|
1.0
|
|
);
|
|
null_space::assign_value_block(direction, layout, null_space::enthalpyValue, enthalpyDirection);
|
|
|
|
null_space::value_view(direction, layout, null_space::bernoulliValue)(0) =
|
|
-null_space::const_value_view(state, layout, null_space::bernoulliValue)(0);
|
|
|
|
return {
|
|
.name = "n3_homology",
|
|
.family = "homology",
|
|
.axis = -1,
|
|
.requiresGravityCompletion = true,
|
|
.direction = std::move(direction)
|
|
};
|
|
}
|
|
|
|
[[nodiscard]] std::vector<GaugeMode> make_gauge_modes(null_space::N3Equilibrium &fixture) {
|
|
auto surfaceModes = null_space::make_surface_modes(fixture);
|
|
const auto homologyMode = std::ranges::find_if(surfaceModes, [](const null_space::SurfaceMode &mode) {
|
|
return mode.kind == null_space::SurfaceModeKind::uniform_radial;
|
|
});
|
|
MFEM_VERIFY(homologyMode != surfaceModes.end(), "The reduced surface modes do not contain homology.");
|
|
|
|
std::vector<GaugeMode> modes;
|
|
modes.reserve(6);
|
|
modes.push_back(make_homology_mode(fixture, *homologyMode));
|
|
|
|
for (auto &surfaceMode : surfaceModes) {
|
|
if (surfaceMode.kind == null_space::SurfaceModeKind::uniform_radial) {
|
|
continue;
|
|
}
|
|
modes.push_back(
|
|
{.name = surfaceMode.name,
|
|
.family = null_space::surface_mode_kind_name(surfaceMode.kind),
|
|
.axis = surfaceMode.axis,
|
|
.requiresGravityCompletion = true,
|
|
.direction = std::move(surfaceMode.direction)}
|
|
);
|
|
}
|
|
return modes;
|
|
}
|
|
|
|
[[nodiscard]] long long global_nonzero_count(
|
|
const mfem::Vector &vector,
|
|
const MPI_Comm communicator
|
|
) {
|
|
long long localCount = 0;
|
|
for (int index = 0; index < vector.Size(); ++index) {
|
|
if (vector(index) != 0.0) {
|
|
++localCount;
|
|
}
|
|
}
|
|
|
|
long long globalCount = 0;
|
|
MPI_Allreduce(&localCount, &globalCount, 1, MPI_LONG_LONG, MPI_SUM, communicator);
|
|
return globalCount;
|
|
}
|
|
|
|
[[nodiscard]] mean_field::models::structure::StructureSeed make_n3_seed(null_space::Model &model) {
|
|
constexpr double surfaceCoordinate = 6.8968486193769603755;
|
|
constexpr int radialSampleCount = 8192;
|
|
const double pi = std::acos(-1.0);
|
|
const double radius = mean_field::utils::RADIUS;
|
|
const double targetMass = mean_field::utils::MASS;
|
|
constexpr double dimensionlessMass = 2.0182359509662283534;
|
|
const double polytropicConstant =
|
|
pi * mean_field::utils::G * std::pow(targetMass / (4.0 * pi * dimensionlessMass), 2.0 / 3.0);
|
|
const double centralDensity =
|
|
std::pow(surfaceCoordinate * std::sqrt(polytropicConstant / (pi * mean_field::utils::G)) / radius, 3.0);
|
|
return model.makeInitialSeed({.centralDensity = centralDensity, .radialSampleCount = radialSampleCount});
|
|
}
|
|
|
|
[[nodiscard]] mfem::Vector project_n3_density(
|
|
const mean_field::fem::FEM &fem,
|
|
const mean_field::models::structure::StructureSeed &seed
|
|
) {
|
|
const auto interpolate = [](const mfem::Vector &radii, const mfem::Vector &values, const double radius) {
|
|
if (radius <= radii(0)) {
|
|
return values(0);
|
|
}
|
|
const int finalIndex = radii.Size() - 1;
|
|
if (radius >= radii(finalIndex)) {
|
|
return values(finalIndex);
|
|
}
|
|
|
|
int lower = 0;
|
|
int upper = finalIndex;
|
|
while (upper - lower > 1) {
|
|
const int middle = lower + (upper - lower) / 2;
|
|
if (radii(middle) <= radius) {
|
|
lower = middle;
|
|
} else {
|
|
upper = middle;
|
|
}
|
|
}
|
|
const double fraction = (radius - radii(lower)) / (radii(upper) - radii(lower));
|
|
return (1.0 - fraction) * values(lower) + fraction * values(upper);
|
|
};
|
|
|
|
mfem::FunctionCoefficient densityCoefficient([&seed, &interpolate](const mfem::Vector &position) {
|
|
const double radius = position.Norml2();
|
|
return radius >= seed.stellarRadius ? 0.0 : interpolate(seed.radius, seed.density, radius);
|
|
});
|
|
mfem::ParGridFunction densityField(fem.densityFes.get());
|
|
densityField.ProjectCoefficient(densityCoefficient);
|
|
|
|
mfem::Vector densityTrue;
|
|
densityField.GetTrueDofs(densityTrue);
|
|
return densityTrue;
|
|
}
|
|
|
|
void run_homology_mass_cancellation_experiment(const int hRefinementLevel) {
|
|
REQUIRE(hRefinementLevel >= 0);
|
|
|
|
mean_field::utils::Args args = test_utils::setup_args();
|
|
mean_field::fem::FEM fem = mean_field::fem::setup_fem(args.mesh_file, args, hRefinementLevel);
|
|
REQUIRE(fem.okay());
|
|
|
|
null_space::Model model = null_space::make_model();
|
|
const mean_field::models::structure::StructureSeed seed = make_n3_seed(model);
|
|
const mfem::Vector densityTrue = project_n3_density(fem, seed);
|
|
|
|
auto deformation = model.compileDomainDeformation(fem);
|
|
const auto &surface = deformation.surfaceDeformationPrescription();
|
|
mfem::Vector zeroSurfaceParameters(surface.parameterCount());
|
|
mfem::Vector homologySurfaceDirection(surface.parameterCount());
|
|
zeroSurfaceParameters = 0.0;
|
|
for (int parameter = 0; parameter < homologySurfaceDirection.Size(); ++parameter) {
|
|
homologySurfaceDirection(parameter) = surface.referenceRadius(parameter);
|
|
}
|
|
|
|
mfem::Vector volumeDirection(deformation.volumeDisplacementSize());
|
|
deformation.applyJacobian(zeroSurfaceParameters, homologySurfaceDirection, volumeDirection);
|
|
mfem::ParGridFunction coordinateVelocity(fem.displacementFes.get());
|
|
coordinateVelocity.SetFromTrueDofs(volumeDirection);
|
|
|
|
mean_field::operators::context::gravity_field::GravityFieldLinearizationContext gravityContext(
|
|
fem, *fem.domainMapperStateless
|
|
);
|
|
const mean_field::field::FieldDofMap &densityMap = gravityContext.GetDensityMap();
|
|
const mfem::Vector reducedDensity = densityMap.gather(densityTrue);
|
|
const mfem::Vector densityDirection = project_extension_aware_homology_scalar(
|
|
*fem.densityFes, densityMap, reducedDensity, coordinateVelocity, surface.referenceCenter(), 3.0
|
|
);
|
|
|
|
mfem::Vector zeroDisplacement(fem.displacementFes->GetTrueVSize());
|
|
mfem::Vector zeroGravityGradient(fem.gravityFluxFes->GetTrueVSize());
|
|
mfem::Vector zeroGravityPotential(fem.gravityPotentialFes->GetTrueVSize());
|
|
zeroDisplacement = 0.0;
|
|
zeroGravityGradient = 0.0;
|
|
zeroGravityPotential = 0.0;
|
|
|
|
const mean_field::operators::MassNormalizationDependencies dependencies{
|
|
.discretization = {.identity = 9101, .revision = 1},
|
|
.density = {.identity = 9103, .revision = 1},
|
|
.displacement = {.identity = 9109, .revision = 1},
|
|
.targetMass = {.identity = 9127, .revision = 1}
|
|
};
|
|
gravityContext.Prepare(
|
|
{.density = reducedDensity,
|
|
.displacement = gravityContext.GetDisplacementMap().gather(zeroDisplacement),
|
|
.gravity_gradient = gravityContext.GetGravityGradientMap().gather(zeroGravityGradient),
|
|
.gravity_potential = gravityContext.GetGravityPotentialMap().gather(zeroGravityPotential)},
|
|
{.discretization = {.value = dependencies.discretization.revision},
|
|
.displacement = {.value = dependencies.displacement.revision},
|
|
.density = {.value = dependencies.density.revision},
|
|
.gravity_gradient = {.value = 1},
|
|
.gravity_potential = {.value = 1}}
|
|
);
|
|
|
|
mean_field::operators::PreparedMassNormalizationOperator massOperator(
|
|
fem, *fem.domainMapperStateless, gravityContext
|
|
);
|
|
massOperator.Prepare({.targetMass = mean_field::utils::MASS}, dependencies);
|
|
|
|
const mfem::Vector reducedVolumeDirection = gravityContext.GetDisplacementMap().gather(volumeDirection);
|
|
const HomologyMassCancellation cancellation =
|
|
measure_homology_mass_cancellation(massOperator, densityDirection, reducedVolumeDirection);
|
|
|
|
std::map<std::string, double> metrics{
|
|
{"density_direction_norm", null_space::global_norm(densityDirection, fem.mesh->GetComm())},
|
|
{"surface_direction_norm", null_space::global_norm(homologySurfaceDirection, fem.mesh->GetComm())},
|
|
{"volume_direction_norm", null_space::global_norm(volumeDirection, fem.mesh->GetComm())},
|
|
{"global_element_count", static_cast<double>(fem.mesh->GetGlobalNE())},
|
|
{"global_density_true_dof_count", static_cast<double>(fem.densityFes->GlobalTrueVSize())},
|
|
{"global_displacement_true_dof_count", static_cast<double>(fem.displacementFes->GlobalTrueVSize())},
|
|
{"global_surface_parameter_count", static_cast<double>(surface.globalParameterCount())}
|
|
};
|
|
add_homology_mass_metrics(metrics, cancellation);
|
|
|
|
int rank = 0;
|
|
MPI_Comm_rank(fem.mesh->GetComm(), &rank);
|
|
if (rank == 0) {
|
|
const int pRefinementLevel = mean_field::field::uniformPolynomialOrderIncrement;
|
|
experiment::record_experiment_result(
|
|
"n3_homology_mass_cancellation",
|
|
"h" + std::to_string(hRefinementLevel) + "_p" + std::to_string(pRefinementLevel),
|
|
{{"h_refinement_level", std::to_string(hRefinementLevel)},
|
|
{"p_refinement_level", std::to_string(pRefinementLevel)},
|
|
{"density_polynomial_order", std::to_string(mean_field::field::Density::Scalar::familyOrder)},
|
|
{"enthalpy_polynomial_order", std::to_string(mean_field::field::Enthalpy::Scalar::familyOrder)},
|
|
{"displacement_polynomial_order",
|
|
std::to_string(mean_field::field::Displacement::Vector::familyOrder)},
|
|
{"gravity_polynomial_order", std::to_string(mean_field::field::Gravity::Potential::familyOrder)},
|
|
{"mesh_file", test_utils::setup_args().mesh_file}},
|
|
std::move(metrics)
|
|
);
|
|
}
|
|
}
|
|
} // namespace
|
|
|
|
TEST_CASE(
|
|
"Coupled Stellar Equilibrium Homology And Reduced Surface Mode Responses",
|
|
"[null_space][surface_modes][conditioning][homology]"
|
|
) {
|
|
mean_field::utils::Args args = test_utils::setup_args();
|
|
args.p.rtol = std::min(args.p.rtol, 1.0e-12);
|
|
args.p.atol = std::min(args.p.atol, 1.0e-13);
|
|
args.p.max_iters = std::max(args.p.max_iters, 1500);
|
|
|
|
null_space::N3Equilibrium fixture(std::move(args));
|
|
const MPI_Comm communicator = fixture.fem().mesh->GetComm();
|
|
int rank = 0;
|
|
MPI_Comm_rank(communicator, &rank);
|
|
|
|
const std::vector<GaugeMode> modes = make_gauge_modes(fixture);
|
|
const auto &stellarOperator = fixture.stellar_operator();
|
|
const auto &layout = stellarOperator.GetLayout();
|
|
|
|
GravityUnknownJacobian gravityJacobian(stellarOperator);
|
|
mean_field::operators::ReducedGravityFieldPreconditioner gravityPreconditioner(
|
|
fixture.fem(), stellarOperator.GetGravityContext().GetGeometryContext()
|
|
);
|
|
|
|
mfem::MINRESSolver gravitySolver(communicator);
|
|
gravitySolver.SetOperator(gravityJacobian);
|
|
gravitySolver.SetPreconditioner(gravityPreconditioner);
|
|
gravitySolver.SetRelTol(1.0e-10);
|
|
gravitySolver.SetAbsTol(1.0e-12);
|
|
gravitySolver.SetMaxIter(1500);
|
|
gravitySolver.SetPrintLevel(0);
|
|
|
|
for (std::size_t modeIndex = 0; modeIndex < modes.size(); ++modeIndex) {
|
|
const GaugeMode &mode = modes[modeIndex];
|
|
null_space::report_progress(
|
|
communicator,
|
|
"evaluating " + mode.name + " (" + std::to_string(modeIndex + 1) + "/" + std::to_string(modes.size()) + ")"
|
|
);
|
|
|
|
const double prescribedInputNorm = null_space::global_norm(mode.direction, communicator);
|
|
REQUIRE(std::isfinite(prescribedInputNorm));
|
|
REQUIRE(prescribedInputNorm > 0.0);
|
|
|
|
const mfem::Vector prescribedAction = fixture.jacobian_action(mode.direction);
|
|
const double prescribedActionNorm = null_space::global_norm(prescribedAction, communicator);
|
|
|
|
GravityCompletionResult completion;
|
|
completion.direction.SetSize(gravityJacobian.Width());
|
|
completion.direction = 0.0;
|
|
|
|
mfem::Vector completedDirection(mode.direction);
|
|
|
|
if (mode.requiresGravityCompletion) {
|
|
completion =
|
|
solve_gravity_completion(prescribedAction, layout, communicator, gravityJacobian, gravitySolver);
|
|
|
|
assign_gravity_completion(
|
|
completedDirection, layout, completion.direction, gravityJacobian.gravity_gradient_size()
|
|
);
|
|
}
|
|
|
|
const mfem::Vector completedAction = fixture.jacobian_action(completedDirection);
|
|
|
|
const double completedInputNorm = null_space::global_norm(completedDirection, communicator);
|
|
const double completedActionNorm = null_space::global_norm(completedAction, communicator);
|
|
const double completionNorm = null_space::global_norm(completion.direction, communicator);
|
|
|
|
REQUIRE(std::isfinite(prescribedActionNorm));
|
|
REQUIRE(std::isfinite(completedInputNorm));
|
|
REQUIRE(std::isfinite(completedActionNorm));
|
|
REQUIRE(completedInputNorm > 0.0);
|
|
|
|
std::map<std::string, double> metrics{
|
|
{"prescribed_input_norm", prescribedInputNorm},
|
|
{"prescribed_action_norm", prescribedActionNorm},
|
|
{"completed_input_norm", completedInputNorm},
|
|
{"completed_action_norm", completedActionNorm},
|
|
{"normalized_completed_response", completedActionNorm / completedInputNorm},
|
|
{"gravity_completion_norm", completionNorm},
|
|
{"gravity_solve_rhs_norm", completion.rightHandSideNorm},
|
|
{"gravity_solve_residual_norm", completion.residualNorm},
|
|
{"gravity_solve_relative_residual", completion.relativeResidual},
|
|
{"gravity_solve_final_norm", completion.finalNorm},
|
|
{"gravity_solve_iterations", static_cast<double>(completion.iterations)},
|
|
{"global_nonzero_input_dofs", static_cast<double>(global_nonzero_count(mode.direction, communicator))}
|
|
};
|
|
|
|
if (mode.family == "homology") {
|
|
const mfem::Vector densityDirection =
|
|
null_space::const_value_view(mode.direction, layout, null_space::densityValue);
|
|
const mfem::Vector volumeDirection = fixture.lifted_surface_direction(mode.direction);
|
|
const HomologyMassCancellation massCancellation = measure_homology_mass_cancellation(
|
|
stellarOperator.GetMassNormalizationOperator(), densityDirection, volumeDirection
|
|
);
|
|
add_homology_mass_metrics(metrics, massCancellation);
|
|
}
|
|
|
|
if (completedActionNorm > 0.0) {
|
|
metrics.emplace("gravity_completion_reduction", prescribedActionNorm / completedActionNorm);
|
|
}
|
|
|
|
add_block_metrics(
|
|
metrics, "prescribed_", null_space::residual_block_norms(prescribedAction, layout, communicator)
|
|
);
|
|
add_block_metrics(
|
|
metrics, "completed_", null_space::residual_block_norms(completedAction, layout, communicator)
|
|
);
|
|
|
|
if (rank == 0) {
|
|
experiment::record_experiment_result(
|
|
"coupled_reduced_surface_mode_conditioning", mode.name,
|
|
{{"mode_family", mode.family},
|
|
{"axis", std::to_string(mode.axis)},
|
|
{"gravity_completion_requested", mode.requiresGravityCompletion ? "true" : "false"},
|
|
{"gravity_solve_performed", completion.solvePerformed ? "true" : "false"},
|
|
{"rotation_fraction_of_keplerian", "0.0"},
|
|
{"mesh_file", test_utils::setup_args().mesh_file},
|
|
{"local_state_dofs", std::to_string(stellarOperator.Width())}},
|
|
std::move(metrics)
|
|
);
|
|
}
|
|
}
|
|
|
|
null_space::report_progress(communicator, "coupled reduced surface-mode probe complete; writing CSV output");
|
|
}
|
|
|
|
TEST_CASE(
|
|
"Fixed Central Density Phase Couples To The N3 Homology Tangent",
|
|
"[null_space][homology][central_density][phase]"
|
|
) {
|
|
mean_field::utils::Args args = test_utils::setup_args();
|
|
null_space::N3Equilibrium fixture(std::move(args));
|
|
const MPI_Comm communicator = fixture.fem().mesh->GetComm();
|
|
|
|
const std::vector<GaugeMode> modes = make_gauge_modes(fixture);
|
|
const auto homology = std::ranges::find_if(modes, [](const GaugeMode &mode) { return mode.family == "homology"; });
|
|
REQUIRE(homology != modes.end());
|
|
|
|
const auto &layout = fixture.stellar_operator().GetLayout();
|
|
const mfem::Vector enthalpy = null_space::const_value_view(fixture.state(), layout, null_space::enthalpyValue);
|
|
const mfem::Vector enthalpyDirection =
|
|
null_space::const_value_view(homology->direction, layout, null_space::enthalpyValue);
|
|
|
|
const mean_field::field::FieldDofMap enthalpyMap =
|
|
mean_field::field::make_field_dof_map<mean_field::field::Enthalpy, null_space::DomainSchema>(
|
|
*fixture.fem().enthalpyFes
|
|
);
|
|
mfem::Vector origin(fixture.fem().mesh->SpaceDimension());
|
|
origin = 0.0;
|
|
mean_field::field::FieldPointDofMap centerDof =
|
|
mean_field::field::make_field_point_dof_map<mean_field::field::Enthalpy>(
|
|
*fixture.fem().enthalpyFes, enthalpyMap, origin, 1.0e-12
|
|
);
|
|
|
|
double localCentralEnthalpy = 0.0;
|
|
for (const int reducedDof : centerDof.reduced_dofs()) {
|
|
localCentralEnthalpy += enthalpy(reducedDof);
|
|
}
|
|
double centralEnthalpy = 0.0;
|
|
MPI_Allreduce(&localCentralEnthalpy, ¢ralEnthalpy, 1, MPI_DOUBLE, MPI_SUM, communicator);
|
|
REQUIRE(std::isfinite(centralEnthalpy));
|
|
REQUIRE(centralEnthalpy > 0.0);
|
|
|
|
const auto &equationOfState = fixture.model().equationOfState();
|
|
const mean_field::eos::DensityValue targetDensity = mean_field::eos::evaluate<mean_field::eos::quantity::Density>(
|
|
equationOfState, mean_field::eos::SpecificEnthalpyValue{centralEnthalpy}
|
|
);
|
|
const mean_field::models::CompiledFixedCentralDensity compiled =
|
|
mean_field::models::compileConstraint(mean_field::models::FixedCentralDensity{targetDensity}, equationOfState);
|
|
mean_field::operators::PreparedCentralDensityConstraint phase(std::move(centerDof), communicator);
|
|
phase.Prepare(compiled, enthalpy, 0.0, {.enthalpy = {.identity = 3251, .revision = 1}});
|
|
|
|
mfem::Vector enthalpyAction(enthalpy.Size());
|
|
mfem::Vector phaseAction(1);
|
|
enthalpyAction = 0.0;
|
|
phaseAction = 0.0;
|
|
phase.ApplyJacobian(
|
|
{.enthalpyVariation = enthalpyDirection, .borderVariation = 0.0},
|
|
{.enthalpyAction = enthalpyAction, .phaseAction = phaseAction}
|
|
);
|
|
|
|
const double couplingScale = std::max(1.0, std::abs(compiled.targetEnthalpy().value()));
|
|
const double enthalpyDirectionNorm = null_space::global_norm(enthalpyDirection, communicator);
|
|
const double homologyDirectionNorm = null_space::global_norm(homology->direction, communicator);
|
|
const double absolutePhaseCoupling = std::abs(phaseAction(0));
|
|
INFO("Central enthalpy = " << centralEnthalpy);
|
|
INFO("N3 homology phase coupling = " << phaseAction(0));
|
|
REQUIRE(std::isfinite(phaseAction(0)));
|
|
REQUIRE(enthalpyDirectionNorm > 0.0);
|
|
REQUIRE(homologyDirectionNorm > 0.0);
|
|
CHECK(absolutePhaseCoupling > 100.0 * std::numeric_limits<double>::epsilon() * couplingScale);
|
|
|
|
int rank = 0;
|
|
MPI_Comm_rank(communicator, &rank);
|
|
if (rank == 0) {
|
|
experiment::record_experiment_result(
|
|
"fixed_central_density_homology_coupling", "n3_homology",
|
|
{{"mesh_file", test_utils::setup_args().mesh_file},
|
|
{"local_state_dofs", std::to_string(fixture.stellar_operator().Width())}},
|
|
{{"target_density", compiled.targetDensity().value()},
|
|
{"target_enthalpy", compiled.targetEnthalpy().value()},
|
|
{"central_enthalpy", centralEnthalpy},
|
|
{"homology_phase_action", phaseAction(0)},
|
|
{"absolute_phase_coupling", absolutePhaseCoupling},
|
|
{"target_scaled_phase_coupling", absolutePhaseCoupling / couplingScale},
|
|
{"enthalpy_direction_norm", enthalpyDirectionNorm},
|
|
{"enthalpy_normalized_phase_coupling", absolutePhaseCoupling / enthalpyDirectionNorm},
|
|
{"homology_direction_norm", homologyDirectionNorm},
|
|
{"state_normalized_phase_coupling", absolutePhaseCoupling / homologyDirectionNorm}}
|
|
);
|
|
}
|
|
}
|
|
|
|
TEST_CASE(
|
|
"N3 Homology Mass Cancellation At The Registered Polynomial Order",
|
|
"[null_space][homology][mass_normalization][convergence][p_refinement]"
|
|
) {
|
|
run_homology_mass_cancellation_experiment(0);
|
|
}
|
|
|
|
TEST_CASE(
|
|
"N3 Homology Mass Cancellation Under Uniform Spatial Refinement",
|
|
"[null_space][homology][mass_normalization][convergence][h_refinement]"
|
|
) {
|
|
run_homology_mass_cancellation_experiment(1);
|
|
}
|