tboudreaux/MeanField
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Compile and verify radial coarse basis through existing FEM surface map

Browse Source
This commit is contained in:
Emily Boudreaux
2026-09-17 16:55:20 -04:00
parent b3bd244024
commit dc8d4b0501
2 changed files with 104 additions and 41 deletions
Download Patch File Download Diff File Expand all files Collapse all files

123
experiments/plateau_coupled_coarse.hpp
View File

@@ -89,31 +89,97 @@ namespace experiment::plateau_diagnostics {
template <class Inverse>
using MaterialBlock = typename RootBlock<Inverse>::StructureComponent::MaterialSurfaceComponent;
template <mean_field::equilibrium::DiscretizedStellarEquilibriumProblem Problem>
[[nodiscard]] mean_field::deformation::PreparedNodalRadialSurface
MakeDefaultRadialSurface(const Problem &problem) {
using namespace mean_field;
const auto &physical = problem.GetPhysicalOperator();
const auto &domain = physical.GetDomainDeformation();
const auto &fem = physical.GetHydrostaticOperator().GetFEM();
const auto dependencies = domain.discretizationDependencies();
const auto communicator = problem.GetCommunicator();
const auto all = [communicator](bool condition) {
int local = condition, global = 0;
if (MPI_Allreduce(&local, &global, 1, MPI_INT, MPI_MIN, communicator) != MPI_SUCCESS)
throw std::runtime_error("Could not validate the radial shape basis geometry.");
return global != 0;
};
if (!all(
physical.IsPrepared() && fem.mesh && fem.surfaceDeformationFes && fem.displacementFes &&
dependencies.isCurrent() && dependencies.physicalMeshIdentity == fem.mesh.get() &&
dependencies.surfaceScalarSpaceIdentity == fem.surfaceDeformationFes.get() &&
dependencies.volumeDisplacementSpaceIdentity == fem.displacementFes.get()
))
throw std::logic_error("Radial shape modes require the current physical deformation spaces.");
const auto &scalar = *fem.surfaceDeformationFes;
const auto &vector = *fem.displacementFes;
mfem::Vector center(3);
center = 0.0;
const deformation::NodalRadialSurface prescription(std::move(center));
if (!all(
fem.mesh->SpaceDimension() == 3 && scalar.GetVDim() == 1 && vector.GetVDim() == 3 &&
vector.GetOrdering() == mfem::Ordering::byNODES &&
vector.GetTrueVSize() == 3 * scalar.GetTrueVSize() && vector.FEColl() == scalar.FEColl() &&
domain.descriptor().surfaceDeformation == prescription.descriptor()
))
throw std::invalid_argument(
"The shape benchmark requires the default radial surface and its scalar companion."
);
// StellarEquilibriumProblem currently compiles this zero-center
// prescription, then erases its type. Reuse the same compiler and
// exact FEM boundary map instead of reproducing its geometry math.
auto radial = deformation::compileSurfaceDeformationPrescription(
prescription,
deformation::SurfaceDeformationCompilationContext{
*fem.surfaceDeformationFes, field::make_stellar_surface_scalar_dof_map(*fem.surfaceDeformationFes)
}
);
if (!all(radial.parameterCount() == domain.parameterCount()))
throw std::invalid_argument("The reconstructed radial map does not match the physical parameters.");
// The runtime descriptor does not retain the center. Verify actual
// surface motion against the reconstructed directions; metadata
// alone must not admit a future shifted-center prescription.
mfem::Vector unit(domain.parameterCount()), action(domain.volumeDisplacementSize());
unit = 1.0;
domain.applyJacobian(physical.GetSurfaceDeformationParameters(), unit, action);
bool matches = true;
const int scalarSize = scalar.GetTrueVSize();
for (int i = 0; i < radial.parameterCount(); ++i)
for (int component = 0; component < 3; ++component) {
// MFEM byNODES stores component-sized blocks, as required
// by the existing radial extension's true-DOF mapping.
const int dof = radial.surfaceDofMap().volume_true_dof(i) + component * scalarSize;
const double observed = action(dof), expected = radial.radialDirection(i, component);
matches = matches && std::isfinite(observed) && std::isfinite(expected) &&
std::abs(observed - expected) <= 512.0 * std::numeric_limits<double>::epsilon();
}
if (!all(matches && dependencies.isCurrent()))
throw std::invalid_argument("The actual surface motion disagrees with the zero-center radial basis.");
return radial;
}
template <class Problem, class Inverse>
concept SupportsShapeModes = requires(const Problem &problem, const Inverse &inverse) {
typename MaterialBlock<Inverse>::Factorization;
requires std::same_as<
typename MaterialBlock<Inverse>::Factorization,
mean_field::preconditioning::SurfaceThenMaterialTriangular>;
{ inverse.GetProblem() } -> std::same_as<const Problem &>;
{ inverse.IsCurrent() } -> std::same_as<bool>;
inverse.GetCoordinateMap().GetCorrectionRanges();
inverse.GetGroupedPreconditioner().GetStructurePreconditioner().GetCrossCouplings();
inverse.GetGroupedPreconditioner()
.GetStructurePreconditioner()
.GetMaterialSurfacePreconditioner()
.GetCoupledOperator();
inverse.GetGroupedPreconditioner().GetStructurePreconditioner().GetGravityPreconditioner();
{ problem.GetPreparationGeneration() } -> std::convertible_to<std::uint64_t>;
{
problem.GetPhysicalOperator().GetDomainDeformation().surfaceDeformationPrescription().radialDirection(
0, 0
)
} -> std::same_as<double>;
{
problem.GetPhysicalOperator().GetDomainDeformation().surfaceDeformationPrescription().referenceRadius(0)
} -> std::same_as<double>;
};
concept SupportsShapeModes =
mean_field::equilibrium::DiscretizedStellarEquilibriumProblem<Problem> &&
requires(const Problem &problem, const Inverse &inverse) {
typename MaterialBlock<Inverse>::Factorization;
requires std::same_as<
typename MaterialBlock<Inverse>::Factorization,
mean_field::preconditioning::SurfaceThenMaterialTriangular>;
{ inverse.GetProblem() } -> std::same_as<const Problem &>;
{ inverse.IsCurrent() } -> std::same_as<bool>;
inverse.GetCoordinateMap().GetCorrectionRanges();
inverse.GetGroupedPreconditioner().GetStructurePreconditioner().GetCrossCouplings();
inverse.GetGroupedPreconditioner()
.GetStructurePreconditioner()
.GetMaterialSurfacePreconditioner()
.GetCoupledOperator();
inverse.GetGroupedPreconditioner().GetStructurePreconditioner().GetGravityPreconditioner();
{ problem.GetPreparationGeneration() } -> std::convertible_to<std::uint64_t>;
{
MakeDefaultRadialSurface(problem)
} -> std::same_as<mean_field::deformation::PreparedNodalRadialSurface>;
};
} // namespace coupled_shape_detail
// These are geometry modes, not Newton-RHS samples. Complete their material
@@ -123,6 +189,7 @@ namespace experiment::plateau_diagnostics {
// Construction performs nine explicit surface-to-material actions, material
// inverse applications, material-to-gravity actions, and gravity applications.
// The material inverse also performs its usual internal coupling actions.
// One extra domain Jacobian action verifies the reconstructed radial basis.
template <
class Problem,
class PhysicalInverse>
@@ -188,10 +255,10 @@ namespace experiment::plateau_diagnostics {
const auto &offsets = couplings.GetOffsets();
const auto &coordinates = inverse.GetCoordinateMap();
const auto &ranges = coordinates.GetCorrectionRanges();
const auto &radial = problem.GetPhysicalOperator().GetDomainDeformation().surfaceDeformationPrescription();
const auto materialOffset = ranges[materialFirst].preconditionerOffset;
const auto gravityOffset = ranges[gravityFirst].preconditionerOffset;
const auto shapeSize = offsets[shapeColumn + 1] - offsets[shapeColumn];
const auto radial = coupled_shape_detail::MakeDefaultRadialSurface(problem);
const auto materialOffset = ranges[materialFirst].preconditionerOffset;
const auto gravityOffset = ranges[gravityFirst].preconditionerOffset;
const auto shapeSize = offsets[shapeColumn + 1] - offsets[shapeColumn];
if (!all(
radial.spatialDimension() == 3 && radial.parameterCount() == shapeSize &&
material.Height() == material.Width() && gravityInverse.Height() == gravityInverse.Width() &&

22
tests/solver/plateau_coupled_coarse.cpp
View File

@@ -163,22 +163,18 @@ namespace coarse_test {
template <class Model>
using Coarse = probe::PreparedCoupledShapeCoarse<NormalizedOperation<Model>, NormalizedInverse<Model>>;
// The experiment borrows the existing radial geometry rather than inventing
// a second nodal-coordinate map. Instantiation below checks actual solver
// types with optional generated borders, without a costly stellar fixture.
// The production problem erases its prepared deformation. Recompile only
// the existing surface prescription, then verify its actual runtime action.
// Explicit instantiations below keep both complete model paths checked.
template <class Model>
concept BorrowsRadialSurface = requires(const Problem<Model> &problem) {
concept ReusesRadialCompiler = requires(const Problem<Model> &problem) {
{
problem.GetPhysicalOperator().GetDomainDeformation().surfaceDeformationPrescription().radialDirection(0, 0)
} -> std::same_as<double>;
{
problem.GetPhysicalOperator().GetDomainDeformation().surfaceDeformationPrescription().referenceRadius(0)
} -> std::same_as<double>;
requires std::is_lvalue_reference_v<
decltype(problem.GetPhysicalOperator().GetDomainDeformation().surfaceDeformationPrescription())>;
probe::coupled_shape_detail::MakeDefaultRadialSurface(problem)
} -> std::same_as<deformation::PreparedNodalRadialSurface>;
requires probe::coupled_shape_detail::SupportsShapeModes<Problem<Model>, PhysicalInverse<Model>>;
};
static_assert(BorrowsRadialSurface<OrdinaryModel>);
static_assert(BorrowsRadialSurface<CentralModel>);
static_assert(ReusesRadialCompiler<OrdinaryModel>);
static_assert(ReusesRadialCompiler<CentralModel>);
static_assert(solver::LinearPreconditionerApplicationContractAvailable<Coarse<OrdinaryModel>>);
static_assert(solver::LinearPreconditionerApplicationContractAvailable<Coarse<CentralModel>>);
} // namespace coarse_test