diff --git a/experiments/plateau_coupled_coarse.hpp b/experiments/plateau_coupled_coarse.hpp index e854c3d..a6bfa98 100644 --- a/experiments/plateau_coupled_coarse.hpp +++ b/experiments/plateau_coupled_coarse.hpp @@ -89,31 +89,97 @@ namespace experiment::plateau_diagnostics { template using MaterialBlock = typename RootBlock::StructureComponent::MaterialSurfaceComponent; + template + [[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::epsilon(); + } + if (!all(matches && dependencies.isCurrent())) + throw std::invalid_argument("The actual surface motion disagrees with the zero-center radial basis."); + return radial; + } + template - concept SupportsShapeModes = requires(const Problem &problem, const Inverse &inverse) { - typename MaterialBlock::Factorization; - requires std::same_as< - typename MaterialBlock::Factorization, - mean_field::preconditioning::SurfaceThenMaterialTriangular>; - { inverse.GetProblem() } -> std::same_as; - { inverse.IsCurrent() } -> std::same_as; - inverse.GetCoordinateMap().GetCorrectionRanges(); - inverse.GetGroupedPreconditioner().GetStructurePreconditioner().GetCrossCouplings(); - inverse.GetGroupedPreconditioner() - .GetStructurePreconditioner() - .GetMaterialSurfacePreconditioner() - .GetCoupledOperator(); - inverse.GetGroupedPreconditioner().GetStructurePreconditioner().GetGravityPreconditioner(); - { problem.GetPreparationGeneration() } -> std::convertible_to; - { - problem.GetPhysicalOperator().GetDomainDeformation().surfaceDeformationPrescription().radialDirection( - 0, 0 - ) - } -> std::same_as; - { - problem.GetPhysicalOperator().GetDomainDeformation().surfaceDeformationPrescription().referenceRadius(0) - } -> std::same_as; - }; + concept SupportsShapeModes = + mean_field::equilibrium::DiscretizedStellarEquilibriumProblem && + requires(const Problem &problem, const Inverse &inverse) { + typename MaterialBlock::Factorization; + requires std::same_as< + typename MaterialBlock::Factorization, + mean_field::preconditioning::SurfaceThenMaterialTriangular>; + { inverse.GetProblem() } -> std::same_as; + { inverse.IsCurrent() } -> std::same_as; + inverse.GetCoordinateMap().GetCorrectionRanges(); + inverse.GetGroupedPreconditioner().GetStructurePreconditioner().GetCrossCouplings(); + inverse.GetGroupedPreconditioner() + .GetStructurePreconditioner() + .GetMaterialSurfacePreconditioner() + .GetCoupledOperator(); + inverse.GetGroupedPreconditioner().GetStructurePreconditioner().GetGravityPreconditioner(); + { problem.GetPreparationGeneration() } -> std::convertible_to; + { + MakeDefaultRadialSurface(problem) + } -> std::same_as; + }; } // 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() && diff --git a/tests/solver/plateau_coupled_coarse.cpp b/tests/solver/plateau_coupled_coarse.cpp index f8eb02a..902a341 100644 --- a/tests/solver/plateau_coupled_coarse.cpp +++ b/tests/solver/plateau_coupled_coarse.cpp @@ -163,22 +163,18 @@ namespace coarse_test { template using Coarse = probe::PreparedCoupledShapeCoarse, NormalizedInverse>; - // 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 - concept BorrowsRadialSurface = requires(const Problem &problem) { + concept ReusesRadialCompiler = requires(const Problem &problem) { { - problem.GetPhysicalOperator().GetDomainDeformation().surfaceDeformationPrescription().radialDirection(0, 0) - } -> std::same_as; - { - problem.GetPhysicalOperator().GetDomainDeformation().surfaceDeformationPrescription().referenceRadius(0) - } -> std::same_as; - requires std::is_lvalue_reference_v< - decltype(problem.GetPhysicalOperator().GetDomainDeformation().surfaceDeformationPrescription())>; + probe::coupled_shape_detail::MakeDefaultRadialSurface(problem) + } -> std::same_as; + requires probe::coupled_shape_detail::SupportsShapeModes, PhysicalInverse>; }; - static_assert(BorrowsRadialSurface); - static_assert(BorrowsRadialSurface); + static_assert(ReusesRadialCompiler); + static_assert(ReusesRadialCompiler); static_assert(solver::LinearPreconditionerApplicationContractAvailable>); static_assert(solver::LinearPreconditionerApplicationContractAvailable>); } // namespace coarse_test