MeanField/tests/physics/gravity.cpp

#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_floating_point.hpp>

#include <algorithm>
#include <array>
#include <cassert>
#include <cmath>
#include <limits>
#include <memory>

#include <boost/math/quadrature/gauss_kronrod.hpp>
#include <mfem.hpp>
#include <mpi.h>

import mean_field;
import test_helpers;

using namespace mean_field;

namespace {
    double global_vector_norm(
        const mfem::Vector &vector,
        MPI_Comm communicator
    ) {
        const double local_norm_squared = vector * vector;
        double global_norm_squared      = 0.0;
        MPI_Allreduce(
            &local_norm_squared, &global_norm_squared, 1, MPI_DOUBLE, MPI_SUM,
            communicator
        );
        return std::sqrt(global_norm_squared);
    }

    double global_vector_dot(
        const mfem::Vector &lhs,
        const mfem::Vector &rhs,
        MPI_Comm communicator
    ) {
        const double local_dot = lhs * rhs;
        double global_dot      = 0.0;
        MPI_Allreduce(
            &local_dot, &global_dot, 1, MPI_DOUBLE, MPI_SUM, communicator
        );
        return global_dot;
    }

    double global_relative_vector_error(
        const mfem::Vector &computed,
        const mfem::Vector &reference,
        MPI_Comm communicator
    ) {
        mfem::Vector difference(computed);
        difference -= reference;
        return global_vector_norm(difference, communicator) /
               std::max(
                   global_vector_norm(reference, communicator),
                   std::numeric_limits<double>::epsilon()
               );
    }

    struct GravitationalEnergies {
        double binding;
        double virial;
    };

    struct HomogeneousEllipsoidAnalytic {
        double coefficient_x;
        double coefficient_y;
        double coefficient_z;
        double energy_kernel;
    };

    class HomogeneousEllipsoidHDivCoefficient : public mfem::VectorCoefficient {
    public:
        HomogeneousEllipsoidHDivCoefficient(
            const mapping::DomainMapper &domain_mapping,
            const double density,
            const HomogeneousEllipsoidAnalytic &analytic
        )
            : VectorCoefficient(3),
              domain_mapping(domain_mapping),
              density(density),
              coefficient_x(analytic.coefficient_x),
              coefficient_y(analytic.coefficient_y),
              coefficient_z(analytic.coefficient_z) {
        }

        void Eval(
            mfem::Vector &value,
            mfem::ElementTransformation &transformation,
            const mfem::IntegrationPoint &integration_point
        ) override {
            transformation.SetIntPoint(&integration_point);

            mfem::Vector x_physical(3);
            mfem::Vector field_physical(3);
            mfem::DenseMatrix map_jacobian(3, 3);
            mfem::DenseMatrix inverse_map_jacobian(3, 3);

            domain_mapping.GetPhysicalPoint(
                transformation, integration_point, x_physical
            );

            field_physical(0) =
                2.0 * M_PI * utils::G * density * coefficient_x * x_physical(0);
            field_physical(1) =
                2.0 * M_PI * utils::G * density * coefficient_y * x_physical(1);
            field_physical(2) =
                2.0 * M_PI * utils::G * density * coefficient_z * x_physical(2);

            domain_mapping.ComputeJacobian(transformation, map_jacobian);
            const double map_determinant = map_jacobian.Det();
            MFEM_VERIFY(
                map_determinant > 0.0,
                "Domain mapping has a non-positive Jacobian determinant."
            );

            mfem::CalcInverse(map_jacobian, inverse_map_jacobian);
            inverse_map_jacobian.Mult(field_physical, value);
            value *= map_determinant;
        }

    private:
        const mapping::DomainMapper &domain_mapping;
        double density;
        double coefficient_x;
        double coefficient_y;
        double coefficient_z;
    };

    template <typename GravitySolutionType>
    GravitationalEnergies compute_gravitational_energies(
        fem::FEM &f,
        const mfem::GridFunction &rho,
        const GravitySolutionType &gravity_solution,
        const int quadrature_order
    ) {
        const int dim                = f.mesh->Dimension();
        double local_bind_integral   = 0.0;
        double local_virial_integral = 0.0;

        mfem::Vector x_physical(dim);
        mfem::Vector grad_phi_element(dim);
        mfem::Vector grad_phi_physical(dim);
        mfem::DenseMatrix map_jacobian(dim, dim);

        for (int elem_id = 0; elem_id < f.mesh->GetNE(); ++elem_id) {
            if (f.mesh->GetAttribute(elem_id) == 3) {
                continue;
            }

            mfem::ElementTransformation *transformation =
                f.mesh->GetElementTransformation(elem_id);
            const mfem::IntegrationRule &integration_rule = mfem::IntRules.Get(
                transformation->GetGeometryType(), quadrature_order
            );

            for (int q = 0; q < integration_rule.GetNPoints(); ++q) {
                const mfem::IntegrationPoint &integration_point =
                    integration_rule.IntPoint(q);
                transformation->SetIntPoint(&integration_point);

                double weight =
                    transformation->Weight() * integration_point.weight;

                if (f.has_mapping()) {
                    const double map_determinant = f.mapping->ComputeDetJ(
                        *transformation, integration_point
                    );
                    MFEM_VERIFY(
                        map_determinant > 0.0,
                        "Domain mapping has a non-positive Jacobian "
                        "determinant."
                    );

                    weight *= map_determinant;
                    f.mapping->GetPhysicalPoint(
                        *transformation, integration_point, x_physical
                    );
                    gravity_solution.gradPhi.GetVectorValue(
                        elem_id, integration_point, grad_phi_element
                    );
                    f.mapping->ComputeJacobian(*transformation, map_jacobian);
                    map_jacobian.Mult(grad_phi_element, grad_phi_physical);
                    grad_phi_physical /= map_determinant;
                } else {
                    transformation->Transform(integration_point, x_physical);
                    gravity_solution.gradPhi.GetVectorValue(
                        elem_id, integration_point, grad_phi_physical
                    );
                }

                const double rho_value =
                    rho.GetValue(elem_id, integration_point);
                const double phi_value =
                    gravity_solution.phi.GetValue(elem_id, integration_point);
                double radius_dot_gradient = 0.0;

                for (int d = 0; d < dim; ++d) {
                    radius_dot_gradient +=
                        (x_physical(d) - f.com(d)) * grad_phi_physical(d);
                }

                local_bind_integral += rho_value * phi_value * weight;
                local_virial_integral +=
                    rho_value * radius_dot_gradient * weight;
            }
        }

        const double local_w_bind = 0.5 * local_bind_integral;
        const double local_w_vir  = -local_virial_integral;
        double global_w_bind      = 0.0;
        double global_w_vir       = 0.0;
        MPI_Comm communicator     = f.densityFes->GetComm();

        MPI_Allreduce(
            &local_w_bind, &global_w_bind, 1, MPI_DOUBLE, MPI_SUM, communicator
        );
        MPI_Allreduce(
            &local_w_vir, &global_w_vir, 1, MPI_DOUBLE, MPI_SUM, communicator
        );

        return {.binding = global_w_bind, .virial = global_w_vir};
    }

    void zero_vacuum_density(
        const fem::FEM &f,
        mfem::GridFunction &rho
    ) {
        for (int i = 0; i < f.vacuumDensityTdofs.Size(); ++i) {
            rho(f.vacuumDensityTdofs[i]) = 0.0;
        }
    }

    int get_gravity_quadrature_order(const fem::FEM &f) {
        return 2 * std::max(
                       f.gravityPotentialFes->GetMaxElementOrder(),
                       f.gravityFluxFes->GetMaxElementOrder()
                   ) +
               8;
    }

    double compute_ellipsoid_coefficient(
        const double normalized_axis_x,
        const double normalized_axis_y,
        const double normalized_axis_z,
        const double target_axis_squared
    ) {
        auto integrand = [=](const double t) {
            if (t <= 0.0 || t >= 1.0) {
                return 0.0;
            }

            const double one_minus_t   = 1.0 - t;
            const double s             = t / one_minus_t;
            const double s_squared     = s * s;
            const double ds_squared_dt = 2.0 * s / (one_minus_t * one_minus_t);
            const double delta         = std::sqrt(
                (normalized_axis_x * normalized_axis_x + s_squared) *
                (normalized_axis_y * normalized_axis_y + s_squared) *
                (normalized_axis_z * normalized_axis_z + s_squared)
            );

            return normalized_axis_x * normalized_axis_y * normalized_axis_z *
                   ds_squared_dt / ((target_axis_squared + s_squared) * delta);
        };

        double integration_error = 0.0;
        return boost::math::quadrature::gauss_kronrod<double, 61>::integrate(
            integrand, 0.0, 1.0, 15, 1.0e-13, &integration_error
        );
    }

    double compute_ellipsoid_energy_kernel(
        const double normalized_axis_x,
        const double normalized_axis_y,
        const double normalized_axis_z,
        const double length_scale
    ) {
        auto integrand = [=](const double t) {
            if (t <= 0.0) {
                return 0.0;
            }
            if (t >= 1.0) {
                return 2.0;
            }

            const double one_minus_t   = 1.0 - t;
            const double s             = t / one_minus_t;
            const double s_squared     = s * s;
            const double ds_squared_dt = 2.0 * s / (one_minus_t * one_minus_t);
            const double delta         = std::sqrt(
                (normalized_axis_x * normalized_axis_x + s_squared) *
                (normalized_axis_y * normalized_axis_y + s_squared) *
                (normalized_axis_z * normalized_axis_z + s_squared)
            );

            return ds_squared_dt / delta;
        };

        double integration_error = 0.0;
        const double dimensionless_integral =
            boost::math::quadrature::gauss_kronrod<double, 61>::integrate(
                integrand, 0.0, 1.0, 15, 1.0e-13, &integration_error
            );

        return dimensionless_integral / length_scale;
    }

    HomogeneousEllipsoidAnalytic compute_homogeneous_ellipsoid_analytic(
        const double semi_axis_x,
        const double semi_axis_y,
        const double semi_axis_z
    ) {
        const double length_scale =
            std::cbrt(semi_axis_x * semi_axis_y * semi_axis_z);
        const double normalized_axis_x = semi_axis_x / length_scale;
        const double normalized_axis_y = semi_axis_y / length_scale;
        const double normalized_axis_z = semi_axis_z / length_scale;
        const double coefficient_x     = compute_ellipsoid_coefficient(
            normalized_axis_x, normalized_axis_y, normalized_axis_z,
            normalized_axis_x * normalized_axis_x
        );
        const double coefficient_y = compute_ellipsoid_coefficient(
            normalized_axis_x, normalized_axis_y, normalized_axis_z,
            normalized_axis_y * normalized_axis_y
        );
        const double coefficient_z = compute_ellipsoid_coefficient(
            normalized_axis_x, normalized_axis_y, normalized_axis_z,
            normalized_axis_z * normalized_axis_z
        );
        const double energy_kernel = compute_ellipsoid_energy_kernel(
            normalized_axis_x, normalized_axis_y, normalized_axis_z,
            length_scale
        );

        return {
            .coefficient_x = coefficient_x,
            .coefficient_y = coefficient_y,
            .coefficient_z = coefficient_z,
            .energy_kernel = energy_kernel
        };
    }

    struct GravitySolutionComparison {
        double relative_gradient_difference;
        double relative_potential_difference;
    };

    template <
        typename LeftGravitySolution,
        typename RightGravitySolution>
    GravitySolutionComparison compare_gravity_solutions(
        fem::FEM &f,
        const LeftGravitySolution &left_solution,
        const RightGravitySolution &right_solution,
        const int quadrature_order
    ) {
        const int dimension                       = f.mesh->Dimension();

        double local_gradient_difference_squared  = 0.0;
        double local_left_gradient_norm_squared   = 0.0;
        double local_right_gradient_norm_squared  = 0.0;
        double local_potential_difference_squared = 0.0;
        double local_left_potential_norm_squared  = 0.0;
        double local_right_potential_norm_squared = 0.0;

        mfem::Vector left_gradient_element(dimension);
        mfem::Vector right_gradient_element(dimension);
        mfem::Vector left_gradient_physical(dimension);
        mfem::Vector right_gradient_physical(dimension);
        mfem::Vector gradient_difference(dimension);
        mfem::DenseMatrix map_jacobian(dimension, dimension);

        for (int element_id = 0; element_id < f.mesh->GetNE(); ++element_id) {
            if (f.mesh->GetAttribute(element_id) == 3) {
                continue;
            }

            mfem::ElementTransformation *transformation =
                f.mesh->GetElementTransformation(element_id);
            const mfem::IntegrationRule &integration_rule = mfem::IntRules.Get(
                transformation->GetGeometryType(), quadrature_order
            );

            for (int q = 0; q < integration_rule.GetNPoints(); ++q) {
                const mfem::IntegrationPoint &integration_point =
                    integration_rule.IntPoint(q);
                transformation->SetIntPoint(&integration_point);

                double weight =
                    transformation->Weight() * integration_point.weight;

                left_solution.gradPhi.GetVectorValue(
                    element_id, integration_point, left_gradient_element
                );
                right_solution.gradPhi.GetVectorValue(
                    element_id, integration_point, right_gradient_element
                );

                if (f.has_mapping()) {
                    const double map_determinant = f.mapping->ComputeDetJ(
                        *transformation, integration_point
                    );
                    MFEM_VERIFY(
                        map_determinant > 0.0,
                        "Domain mapping has a non-positive Jacobian "
                        "determinant."
                    );

                    weight *= map_determinant;
                    f.mapping->ComputeJacobian(*transformation, map_jacobian);

                    map_jacobian.Mult(
                        left_gradient_element, left_gradient_physical
                    );
                    map_jacobian.Mult(
                        right_gradient_element, right_gradient_physical
                    );

                    left_gradient_physical /= map_determinant;
                    right_gradient_physical /= map_determinant;
                } else {
                    left_gradient_physical  = left_gradient_element;
                    right_gradient_physical = right_gradient_element;
                }

                gradient_difference = left_gradient_physical;
                gradient_difference -= right_gradient_physical;

                local_gradient_difference_squared +=
                    (gradient_difference * gradient_difference) * weight;
                local_left_gradient_norm_squared +=
                    (left_gradient_physical * left_gradient_physical) * weight;
                local_right_gradient_norm_squared +=
                    (right_gradient_physical * right_gradient_physical) *
                    weight;

                const double left_potential =
                    left_solution.phi.GetValue(element_id, integration_point);
                const double right_potential =
                    right_solution.phi.GetValue(element_id, integration_point);
                const double potential_difference =
                    left_potential - right_potential;

                local_potential_difference_squared +=
                    potential_difference * potential_difference * weight;
                local_left_potential_norm_squared +=
                    left_potential * left_potential * weight;
                local_right_potential_norm_squared +=
                    right_potential * right_potential * weight;
            }
        }

        const std::array<double, 6> local_values{
            local_gradient_difference_squared,
            local_left_gradient_norm_squared,
            local_right_gradient_norm_squared,
            local_potential_difference_squared,
            local_left_potential_norm_squared,
            local_right_potential_norm_squared
        };

        std::array<double, 6> global_values{};

        MPI_Allreduce(
            local_values.data(), global_values.data(),
            static_cast<int>(local_values.size()), MPI_DOUBLE, MPI_SUM,
            f.densityFes->GetComm()
        );

        const double gradient_scale_squared =
            0.5 * (global_values[1] + global_values[2]);
        const double potential_scale_squared =
            0.5 * (global_values[4] + global_values[5]);

        MFEM_VERIFY(
            gradient_scale_squared > 0.0,
            "Cannot compare gravity solutions with zero gradient norm."
        );
        MFEM_VERIFY(
            potential_scale_squared > 0.0,
            "Cannot compare gravity solutions with zero potential norm."
        );

        return {
            .relative_gradient_difference =
                std::sqrt(global_values[0] / gradient_scale_squared),
            .relative_potential_difference =
                std::sqrt(global_values[3] / potential_scale_squared)
        };
    }

    struct GravityResidualMetrics {
        double relative_total;
        double relative_gradient;
        double relative_poisson;
    };

    double gravity_test_global_norm(
        const mfem::Vector &vector,
        MPI_Comm communicator
    ) {
        const double local_norm_squared = vector * vector;
        double global_norm_squared      = 0.0;
        MPI_Allreduce(
            &local_norm_squared, &global_norm_squared, 1, MPI_DOUBLE, MPI_SUM,
            communicator
        );
        return std::sqrt(global_norm_squared);
    }

    void gravity_test_get_true_dofs(
        const mfem::ParFiniteElementSpace &finite_element_space,
        const mfem::GridFunction &grid_function,
        mfem::Vector &true_dofs
    ) {
        MFEM_VERIFY(
            grid_function.Size() == finite_element_space.GetVSize(),
            "Grid function has the wrong local size."
        );

        true_dofs.SetSize(finite_element_space.GetTrueVSize());

        const mfem::Operator *restriction =
            finite_element_space.GetRestrictionMatrix();

        if (restriction != nullptr) {
            restriction->Mult(grid_function, true_dofs);
        } else {
            MFEM_VERIFY(
                grid_function.Size() == true_dofs.Size(),
                "Local and true sizes do not match."
            );
            true_dofs = grid_function;
        }
    }

    GravityResidualMetrics measure_gravity_residual(
        const mfem::Vector &residual,
        const mfem::Array<int> &offsets,
        const double right_hand_side_norm,
        MPI_Comm communicator
    ) {
        MFEM_VERIFY(
            offsets.Size() == 3, "Gravity residual must contain two blocks."
        );
        MFEM_VERIFY(
            residual.Size() == offsets.Last(),
            "Gravity residual has the wrong size."
        );
        MFEM_VERIFY(
            right_hand_side_norm > 0.0,
            "Gravity right-hand side must be nonzero."
        );

        mfem::Vector gradient_residual(offsets[1] - offsets[0]);
        mfem::Vector poisson_residual(offsets[2] - offsets[1]);

        for (int i = 0; i < gradient_residual.Size(); ++i) {
            gradient_residual(i) = residual(offsets[0] + i);
        }

        for (int i = 0; i < poisson_residual.Size(); ++i) {
            poisson_residual(i) = residual(offsets[1] + i);
        }

        return {
            .relative_total = gravity_test_global_norm(residual, communicator) /
                              right_hand_side_norm,
            .relative_gradient =
                gravity_test_global_norm(gradient_residual, communicator) /
                right_hand_side_norm,
            .relative_poisson =
                gravity_test_global_norm(poisson_residual, communicator) /
                right_hand_side_norm
        };
    }

    enum class ExteriorMonopoleMapping { legacy, stateless };

    struct ExteriorMonopoleShellMetrics {
        long long quadrature_points{0};
        double minimum_radius{std::numeric_limits<double>::infinity()};
        double maximum_radius{0.0};
        double potential_rms_error{0.0};
        double radial_field_rms_error{0.0};
        double tangential_field_rms{0.0};
    };

    struct ExteriorMonopoleShellAccumulator {
        long long quadrature_points{0};
        double minimum_radius{std::numeric_limits<double>::infinity()};
        double maximum_radius{0.0};
        double reference_weight{0.0};
        double potential_error_squared{0.0};
        double radial_field_error_squared{0.0};
        double tangential_field_squared{0.0};
    };

    constexpr std::array<double, 6> exterior_shell_boundaries{0.0,  0.25, 0.50,
                                                              0.75, 0.90, 1.0};

    int get_exterior_shell(const double compactification_coordinate) {
        REQUIRE(std::isfinite(compactification_coordinate));
        REQUIRE(compactification_coordinate >= -1.0e-12);
        REQUIRE(compactification_coordinate <= 1.0 + 1.0e-12);

        const double coordinate = std::clamp(
            compactification_coordinate, 0.0, std::nextafter(1.0, 0.0)
        );

        for (int shell = 0;
             shell < static_cast<int>(exterior_shell_boundaries.size()) - 1;
             ++shell) {
            if (coordinate < exterior_shell_boundaries[shell + 1]) {
                return shell;
            }
        }

        return static_cast<int>(exterior_shell_boundaries.size()) - 2;
    }

    std::array<
        ExteriorMonopoleShellMetrics,
        5>
    measure_exterior_monopole_shells(
        fem::FEM &f,
        const physics::GravitySolution &solution,
        const mfem::GridFunction &displacement,
        const ExteriorMonopoleMapping mapping_path,
        const double mass
    ) {
        REQUIRE(f.mesh != nullptr);
        REQUIRE(f.gravityFluxFes != nullptr);
        REQUIRE(f.displacementFes != nullptr);
        REQUIRE(f.compactificationFes != nullptr);
        REQUIRE(f.compactificationCoordinate != nullptr);
        REQUIRE(f.mapping != nullptr);
        REQUIRE(f.domainMapperStateless != nullptr);

        constexpr int shell_count =
            static_cast<int>(exterior_shell_boundaries.size()) - 1;

        std::array<ExteriorMonopoleShellAccumulator, shell_count>
            local_shells{};

        mapping::DomainMapperStateless::Workspace workspace(
            f.mesh->Dimension()
        );

        const int vacuum_attribute =
            f.domainMapperStateless->GetVacuumElementAttribute();

        const int quadrature_order = get_gravity_quadrature_order(f);

        for (int element_id = 0; element_id < f.mesh->GetNE(); ++element_id) {
            mfem::ElementTransformation *transformation =
                f.mesh->GetElementTransformation(element_id);

            REQUIRE(transformation != nullptr);

            if (transformation->Attribute != vacuum_attribute) {
                continue;
            }

            const mfem::FiniteElement &displacement_element =
                *f.displacementFes->GetFE(element_id);

            const mfem::FiniteElement &compactification_element =
                *f.compactificationFes->GetFE(element_id);

            mfem::Array<int> displacement_dofs;
            mfem::Array<int> compactification_dofs;

            mfem::DofTransformation *displacement_dof_transformation =
                f.displacementFes->GetElementVDofs(
                    element_id, displacement_dofs
                );

            mfem::DofTransformation *compactification_dof_transformation =
                f.compactificationFes->GetElementDofs(
                    element_id, compactification_dofs
                );

            mfem::Vector element_displacement;
            mfem::Vector element_compactification;

            displacement.GetSubVector(displacement_dofs, element_displacement);

            f.compactificationCoordinate->GetSubVector(
                compactification_dofs, element_compactification
            );

            if (displacement_dof_transformation != nullptr) {
                displacement_dof_transformation->InvTransformPrimal(
                    element_displacement
                );
            }

            if (compactification_dof_transformation != nullptr) {
                compactification_dof_transformation->InvTransformPrimal(
                    element_compactification
                );
            }

            const mapping::ElementDisplacementData displacement_data(
                displacement_element, element_displacement,
                f.displacementFes->GetOrdering()
            );

            const mapping::ElementCompactificationData compactification_data(
                compactification_element, element_compactification
            );

            const mapping::ElementMappingData mapping_data{
                .displacement     = displacement_data,
                .compactification = compactification_data
            };

            mfem::Vector compactification_shape(
                compactification_element.GetDof()
            );

            const mfem::IntegrationRule &integration_rule = mfem::IntRules.Get(
                transformation->GetGeometryType(), quadrature_order
            );

            for (int q = 0; q < integration_rule.GetNPoints(); ++q) {
                const mfem::IntegrationPoint &integration_point =
                    integration_rule.IntPoint(q);

                transformation->SetIntPoint(&integration_point);

                compactification_element.CalcShape(
                    integration_point, compactification_shape
                );

                const double compactification_coordinate =
                    element_compactification * compactification_shape;

                const int shell =
                    get_exterior_shell(compactification_coordinate);

                mfem::Vector reference_field(3);
                mfem::Vector physical_field(3);
                mfem::Vector physical_position(3);

                solution.gradPhi.GetVectorValue(
                    element_id, integration_point, reference_field
                );

                if (mapping_path == ExteriorMonopoleMapping::stateless) {
                    mapping::VolumeMappingContext mapping_context;

                    const mapping::MappingStatus status =
                        f.domainMapperStateless->EvaluateVolume(
                            mapping_data, *transformation, integration_point,
                            workspace, mapping_context
                        );

                    CAPTURE(
                        element_id, q, compactification_coordinate,
                        static_cast<int>(status)
                    );

                    REQUIRE(
                        status == mean_field::mapping::MappingStatus::valid
                    );

                    physical_position =
                        mapping_context.mapping.physical_position;

                    mean_field::mapping::MapHDivFluxToPhysical(
                        mapping_context.mapping, reference_field, physical_field
                    );
                } else {
                    /*
                     * The legacy path intentionally does not use the exterior
                     * coordinate to construct its physical mapping.
                     */
                    f.mapping->GetPhysicalPoint(
                        *transformation, integration_point, physical_position
                    );

                    mfem::DenseMatrix mapping_jacobian(3);
                    f.mapping->ComputeJacobian(
                        *transformation, mapping_jacobian
                    );

                    const double mapping_determinant = mapping_jacobian.Det();

                    REQUIRE(std::isfinite(mapping_determinant));
                    REQUIRE(mapping_determinant > 0.0);

                    mapping_jacobian.Mult(reference_field, physical_field);

                    physical_field /= mapping_determinant;
                }

                const double radius = physical_position.Norml2();

                CAPTURE(
                    element_id, q, shell, compactification_coordinate, radius
                );

                REQUIRE(std::isfinite(radius));
                REQUIRE(radius > 0.0);

                mfem::Vector radial_unit_vector(physical_position);
                radial_unit_vector /= radius;

                const double numerical_radial_field =
                    physical_field * radial_unit_vector;

                mfem::Vector tangential_field(physical_field);
                tangential_field.Add(
                    -numerical_radial_field, radial_unit_vector
                );

                const double numerical_potential =
                    solution.phi.GetValue(element_id, integration_point);

                /*
                 * For an exterior monopole:
                 *
                 *     phi       = -GM/r
                 *     grad(phi) =  GM r_hat/r^2
                 *
                 * These scaled quantities should therefore be one, one, and
                 * zero respectively. They remain well-conditioned as r -> inf.
                 */
                const double scaled_potential =
                    -radius * numerical_potential / (utils::G * mass);

                const double scaled_radial_field = radius * radius *
                                                   numerical_radial_field /
                                                   (utils::G * mass);

                const double scaled_tangential_field =
                    radius * radius * tangential_field.Norml2() /
                    (utils::G * mass);

                REQUIRE(std::isfinite(scaled_potential));
                REQUIRE(std::isfinite(scaled_radial_field));
                REQUIRE(std::isfinite(scaled_tangential_field));

                /*
                 * Use the finite reference-domain measure for averaging. A
                 * physical L2 norm of phi over an infinite three-dimensional
                 * exterior domain is not finite.
                 */
                const double reference_weight =
                    integration_point.weight * transformation->Weight();

                ExteriorMonopoleShellAccumulator &accumulator =
                    local_shells[shell];

                ++accumulator.quadrature_points;

                accumulator.minimum_radius =
                    std::min(accumulator.minimum_radius, radius);

                accumulator.maximum_radius =
                    std::max(accumulator.maximum_radius, radius);

                accumulator.reference_weight += reference_weight;

                accumulator.potential_error_squared +=
                    reference_weight * std::pow(scaled_potential - 1.0, 2);

                accumulator.radial_field_error_squared +=
                    reference_weight * std::pow(scaled_radial_field - 1.0, 2);

                accumulator.tangential_field_squared +=
                    reference_weight * scaled_tangential_field *
                    scaled_tangential_field;
            }
        }

        MPI_Comm communicator = f.gravityFluxFes->GetComm();

        std::array<ExteriorMonopoleShellMetrics, shell_count> metrics{};

        for (int shell = 0; shell < shell_count; ++shell) {
            long long global_points = 0;

            MPI_Allreduce(
                &local_shells[shell].quadrature_points, &global_points, 1,
                MPI_LONG_LONG, MPI_SUM, communicator
            );

            double local_sums[4]{
                local_shells[shell].reference_weight,
                local_shells[shell].potential_error_squared,
                local_shells[shell].radial_field_error_squared,
                local_shells[shell].tangential_field_squared
            };

            double global_sums[4]{};

            MPI_Allreduce(
                local_sums, global_sums, 4, MPI_DOUBLE, MPI_SUM, communicator
            );

            double global_minimum_radius = 0.0;
            double global_maximum_radius = 0.0;

            MPI_Allreduce(
                &local_shells[shell].minimum_radius, &global_minimum_radius, 1,
                MPI_DOUBLE, MPI_MIN, communicator
            );

            MPI_Allreduce(
                &local_shells[shell].maximum_radius, &global_maximum_radius, 1,
                MPI_DOUBLE, MPI_MAX, communicator
            );

            REQUIRE(global_points > 0);
            REQUIRE(global_sums[0] > 0.0);

            metrics[shell] = {
                .quadrature_points = global_points,
                .minimum_radius    = global_minimum_radius,
                .maximum_radius    = global_maximum_radius,
                .potential_rms_error =
                    std::sqrt(global_sums[1] / global_sums[0]),
                .radial_field_rms_error =
                    std::sqrt(global_sums[2] / global_sums[0]),
                .tangential_field_rms =
                    std::sqrt(global_sums[3] / global_sums[0])
            };
        }

        return metrics;
    }

    class StatelessProjectionGeometry {
    public:
        StatelessProjectionGeometry(
            const fem::FEM &f,
            const mapping::DomainMapperStateless &domain_mapper,
            const mfem::GridFunction &displacement
        )
            : m_fem(f),
              m_domain_mapper(domain_mapper),
              m_displacement(displacement),
              m_workspace(f.mesh->Dimension()) {
        }

        mapping::MappingStatus Evaluate(
            mfem::ElementTransformation &transformation,
            const mfem::IntegrationPoint &integration_point,
            mapping::MappingPointContext &context,
            const bool permit_infinity_limit
        ) {
            m_last_evaluation_used_infinity_limit = false;

            const int element_id                  = transformation.ElementNo;

            MFEM_VERIFY(
                element_id >= 0,
                "Projection coefficient received an invalid element number."
            );

            const mfem::FiniteElement &displacement_element =
                *m_fem.displacementFes->GetFE(element_id);

            const mfem::FiniteElement &compactification_element =
                *m_fem.compactificationFes->GetFE(element_id);

            mfem::Array<int> displacement_dofs;
            mfem::Array<int> compactification_dofs;

            mfem::DofTransformation *displacement_dof_transformation =
                m_fem.displacementFes->GetElementVDofs(
                    element_id, displacement_dofs
                );

            mfem::DofTransformation *compactification_dof_transformation =
                m_fem.compactificationFes->GetElementDofs(
                    element_id, compactification_dofs
                );

            mfem::Vector element_displacement;
            mfem::Vector element_compactification;

            m_displacement.GetSubVector(
                displacement_dofs, element_displacement
            );

            m_fem.compactificationCoordinate->GetSubVector(
                compactification_dofs, element_compactification
            );

            if (displacement_dof_transformation != nullptr) {
                displacement_dof_transformation->InvTransformPrimal(
                    element_displacement
                );
            }

            if (compactification_dof_transformation != nullptr) {
                compactification_dof_transformation->InvTransformPrimal(
                    element_compactification
                );
            }

            const mapping::ElementDisplacementData displacement_data(
                displacement_element, element_displacement,
                m_fem.displacementFes->GetOrdering()
            );

            const mapping::ElementCompactificationData compactification_data(
                compactification_element, element_compactification
            );

            mfem::Vector requested_compactification_shape(
                compactification_element.GetDof()
            );

            compactification_element.CalcShape(
                integration_point, requested_compactification_shape
            );

            const double requested_compactification_coordinate =
                element_compactification * requested_compactification_shape;

            constexpr double infinity_candidate_tolerance = 1.0e-8;

            const bool requested_infinity_limit =
                std::isfinite(requested_compactification_coordinate) &&
                requested_compactification_coordinate >=
                    1.0 - infinity_candidate_tolerance;

            const mapping::ElementMappingData mapping_data{
                .displacement     = displacement_data,
                .compactification = compactification_data
            };

            mapping::MappingStatus status = m_domain_mapper.EvaluatePoint(
                mapping_data, transformation, integration_point, m_workspace,
                context
            );

            if (status == mapping::MappingStatus::valid) {
                transformation.SetIntPoint(&integration_point);
                return status;
            }

            if (!permit_infinity_limit ||
                !m_domain_mapper.IsCompactifiedElement(transformation)) {
                transformation.SetIntPoint(&integration_point);
                return status;
            }

            const bool retryable_boundary_status =
                status == mapping::MappingStatus::at_compactified_infinity ||
                status == mapping::MappingStatus::outside_reference_domain ||
                status == mapping::MappingStatus::non_finite_result ||
                (requested_infinity_limit &&
                 status == mapping::MappingStatus::non_positive_determinant);

            if (!retryable_boundary_status) {
                transformation.SetIntPoint(&integration_point);
                return status;
            }

            const mfem::IntegrationPoint &element_center =
                mfem::Geometries.GetCenter(transformation.GetGeometryType());

            /*
             * Use the nearest admissible point. Starting extremely close to
             * the requested point preserves the limiting RT trace, while the
             * larger fallbacks accommodate the mapper's infinity guard.
             */
            constexpr std::array<double, 9> inward_fractions{
                1.0e-12, 1.0e-11, 1.0e-10, 1.0e-9, 1.0e-8,
                1.0e-7,  1.0e-6,  1.0e-5,  1.0e-4
            };

            for (const double inward_fraction : inward_fractions) {
                mfem::IntegrationPoint inward_point;

                inward_point.x = (1.0 - inward_fraction) * integration_point.x +
                                 inward_fraction * element_center.x;

                inward_point.y = (1.0 - inward_fraction) * integration_point.y +
                                 inward_fraction * element_center.y;

                inward_point.z = (1.0 - inward_fraction) * integration_point.z +
                                 inward_fraction * element_center.z;

                inward_point.weight = integration_point.weight;

                status              = m_domain_mapper.EvaluatePoint(
                    mapping_data, transformation, inward_point, m_workspace,
                    context
                );

                if (status == mapping::MappingStatus::valid) {
                    m_last_evaluation_used_infinity_limit = true;
                    transformation.SetIntPoint(&integration_point);
                    return status;
                }

                const bool still_retryable =
                    status ==
                        mapping::MappingStatus::at_compactified_infinity ||
                    status ==
                        mapping::MappingStatus::outside_reference_domain ||
                    status == mapping::MappingStatus::non_finite_result ||
                    (requested_infinity_limit &&
                     status ==
                         mapping::MappingStatus::non_positive_determinant);
                if (!still_retryable) {
                    break;
                }
            }

            transformation.SetIntPoint(&integration_point);
            return status;
        }

        [[nodiscard]]
        bool LastEvaluationUsedInfinityLimit() const noexcept {
            return m_last_evaluation_used_infinity_limit;
        }

    private:
        const fem::FEM &m_fem;
        const mapping::DomainMapperStateless &m_domain_mapper;
        const mfem::GridFunction &m_displacement;
        mapping::DomainMapperStateless::Workspace m_workspace;
        bool m_last_evaluation_used_infinity_limit{false};
    };
    class StatelessMonopolePotentialCoefficient final
        : public mfem::Coefficient {
    public:
        StatelessMonopolePotentialCoefficient(
            const fem::FEM &f,
            const mapping::DomainMapperStateless &domain_mapper,
            const mfem::GridFunction &displacement,
            const double mass,
            const double stellar_radius
        )
            : m_geometry(
                  f,
                  domain_mapper,
                  displacement
              ),
              m_vacuum_attribute(domain_mapper.GetVacuumElementAttribute()),
              m_mass(mass),
              m_stellar_radius(stellar_radius) {
        }

        double Eval(
            mfem::ElementTransformation &transformation,
            const mfem::IntegrationPoint &integration_point
        ) override {
            mapping::MappingPointContext context;

            const mapping::MappingStatus status = m_geometry.Evaluate(
                transformation, integration_point, context, true
            );

            MFEM_VERIFY(
                status == mean_field::mapping::MappingStatus::valid,
                "Stateless monopole-potential projection failed."
                    << "\nMapping status = " << static_cast<int>(status)
                    << "\nElement ID = " << transformation.ElementNo
                    << "\nElement attribute = " << transformation.Attribute
                    << "\nIntegration point = <" << integration_point.x << ", "
                    << integration_point.y << ", " << integration_point.z << ">"
            );

            if (m_geometry.LastEvaluationUsedInfinityLimit()) {
                return 0.0;
            }

            const double radius = context.physical_position.Norml2();

            MFEM_VERIFY(
                std::isfinite(radius) && radius > 0.0,
                "Monopole projection encountered an invalid physical radius."
            );

            if (transformation.Attribute == m_vacuum_attribute) {
                return -utils::G * m_mass / radius;
            }

            return -utils::G * m_mass /
                   (2.0 * m_stellar_radius * m_stellar_radius *
                    m_stellar_radius) *
                   (3.0 * m_stellar_radius * m_stellar_radius -
                    radius * radius);
        }

    private:
        StatelessProjectionGeometry m_geometry;
        int m_vacuum_attribute;
        double m_mass;
        double m_stellar_radius;
    };

    class StatelessMonopoleHDivCoefficient final
        : public mfem::VectorCoefficient {
    public:
        StatelessMonopoleHDivCoefficient(
            const fem::FEM &f,
            const mapping::DomainMapperStateless &domain_mapper,
            const mfem::GridFunction &displacement,
            const double mass,
            const double stellar_radius
        )
            : VectorCoefficient(f.mesh->Dimension()),
              m_geometry(
                  f,
                  domain_mapper,
                  displacement
              ),
              m_vacuum_attribute(domain_mapper.GetVacuumElementAttribute()),
              m_mass(mass),
              m_stellar_radius(stellar_radius) {
        }

        void Eval(
            mfem::Vector &value,
            mfem::ElementTransformation &transformation,
            const mfem::IntegrationPoint &integration_point
        ) override {
            mapping::MappingPointContext context;

            const mapping::MappingStatus status = m_geometry.Evaluate(
                transformation, integration_point, context, true
            );

            MFEM_VERIFY(
                status == mean_field::mapping::MappingStatus::valid,
                "Stateless monopole H(div) projection failed."
                    << "\nMapping status = " << static_cast<int>(status)
                    << "\nElement ID = " << transformation.ElementNo
                    << "\nElement attribute = " << transformation.Attribute
                    << "\nIntegration point = <" << integration_point.x << ", "
                    << integration_point.y << ", " << integration_point.z << ">"
                    << "\nInfinity-limit evaluation attempted = "
                    << m_geometry.LastEvaluationUsedInfinityLimit()
            );
            const mfem::Vector &displaced_position = context.displaced_position;

            const double computational_radius = displaced_position.Norml2();

            MFEM_VERIFY(
                std::isfinite(computational_radius) &&
                    computational_radius > 0.0,
                "Monopole H(div) projection encountered an invalid displaced "
                "computational radius."
            );

            const double displacement_determinant =
                context.displacement_jacobian.Det();

            MFEM_VERIFY(
                std::isfinite(displacement_determinant) &&
                    displacement_determinant > 0.0,
                "Monopole H(div) projection encountered an invalid "
                "displacement "
                "Jacobian determinant."
            );

            mfem::DenseMatrix inverse_displacement_jacobian;

            inverse_displacement_jacobian.SetSize(
                context.displacement_jacobian.Height(),
                context.displacement_jacobian.Width()
            );

            mfem::CalcInverse(
                context.displacement_jacobian, inverse_displacement_jacobian
            );

            /*
             * Pull the radial field back only through the regular displacement
             * map.
             *
             * In the compactified vacuum, the Kelvin scale and its radial
             * derivative cancel exactly from the three-dimensional H(div) Piola
             * pullback of the inverse-square monopole field:
             *
             *   det(J) J^{-1} (GM x / |x|^3)
             *     = GM det(A) A^{-1} y / |y|^3.
             *
             * This is also the finite reference-space limit at compactified
             * infinity.
             */
            inverse_displacement_jacobian.Mult(displaced_position, value);

            double radial_denominator = 0.0;

            if (transformation.Attribute == m_vacuum_attribute) {
                radial_denominator = computational_radius *
                                     computational_radius *
                                     computational_radius;
            } else {
                radial_denominator =
                    m_stellar_radius * m_stellar_radius * m_stellar_radius;
            }

            value *= mean_field::utils::G * m_mass * displacement_determinant /
                     radial_denominator;

            for (int component = 0; component < value.Size(); ++component) {
                MFEM_VERIFY(
                    std::isfinite(value(component)),
                    "Monopole H(div) projection produced a non-finite "
                    "reference flux."
                );
            }
        }

    private:
        StatelessProjectionGeometry m_geometry;
        int m_vacuum_attribute;
        double m_mass;
        double m_stellar_radius;
    };

} // namespace

TEST_CASE(
    "Uniform Potential Matches Analytic",
    tags::gravity &tags::analytic_comparison &tags::initialization
) {
    auto args  = test_utils::setup_args();
    fem::FEM f = fem::setup_fem(args.mesh_file, args, 0);
    f.mapping->ResetDisplacement();
    physics::update_stiffness_matrix(f);

    const double radius          = utils::RADIUS;
    const double mass            = utils::MASS;
    const double analytic_volume = (4.0 / 3.0) * M_PI * std::pow(radius, 3.0);
    const double density         = mass / analytic_volume;

    mfem::GridFunction rho_uniform(f.densityFes.get());
    rho_uniform = density;
    zero_vacuum_density(f, rho_uniform);
    analysis::conserve_mass(f, rho_uniform, mass);

    f.com = analysis::get_com(f, rho_uniform);
    f.Q   = physics::compute_quadrupole_moment_tensor(f, rho_uniform, f.com);

    const auto gravity_solution = physics::grav_potential(f, args, rho_uniform);
    constexpr double potential_tolerance =
        utils::APPROX_MAX_ACCEPTABLE_POTENTIAL_ERROR_SI_BURNING;
    double local_max_abs_error     = 0.0;
    double local_max_rel_error     = 0.0;

    const int num_elements_to_test = std::min(30, f.mesh->GetNE());
    for (int elem_id = 0; elem_id < num_elements_to_test; ++elem_id) {
        mfem::ElementTransformation *transformation =
            f.mesh->GetElementTransformation(elem_id);
        const mfem::IntegrationRule &integration_rule =
            mfem::IntRules.Get(transformation->GetGeometryType(), 2);
        const mfem::IntegrationPoint &integration_point =
            integration_rule.IntPoint(0);
        transformation->SetIntPoint(&integration_point);

        mfem::Vector x_physical;
        f.mapping->GetPhysicalPoint(
            *transformation, integration_point, x_physical
        );

        const double radial_coordinate = x_physical.Norml2();
        if (radial_coordinate < 1.0e-9) {
            continue;
        }

        const double phi_analytic =
            -(utils::G * mass / (2.0 * std::pow(radius, 3.0))) *
            (3.0 * radius * radius - radial_coordinate * radial_coordinate);
        const double phi_fem =
            gravity_solution.phi.GetValue(elem_id, integration_point);
        const double absolute_error = std::abs(phi_fem - phi_analytic);
        const double relative_error = absolute_error / std::abs(phi_analytic);

        local_max_abs_error = std::max(local_max_abs_error, absolute_error);
        local_max_rel_error = std::max(local_max_rel_error, relative_error);
        CHECK_THAT(
            relative_error,
            Catch::Matchers::WithinAbs(0.0, 0.1 * potential_tolerance)
        );
    }

    double global_max_abs_error = 0.0;
    double global_max_rel_error = 0.0;
    MPI_Comm communicator       = f.densityFes->GetComm();
    MPI_Allreduce(
        &local_max_abs_error, &global_max_abs_error, 1, MPI_DOUBLE, MPI_MAX,
        communicator
    );
    MPI_Allreduce(
        &local_max_rel_error, &global_max_rel_error, 1, MPI_DOUBLE, MPI_MAX,
        communicator
    );

    const int quadrature_order           = get_gravity_quadrature_order(f);
    const GravitationalEnergies energies = compute_gravitational_energies(
        f, rho_uniform, gravity_solution, quadrature_order
    );
    const double analytic_binding_energy =
        -(3.0 / 5.0) * utils::G * mass * mass / radius;
    const double relative_binding_error =
        std::abs(energies.binding - analytic_binding_energy) /
        std::abs(analytic_binding_energy);
    const double relative_virial_error =
        std::abs(energies.virial - analytic_binding_energy) /
        std::abs(analytic_binding_energy);
    const double relative_consistency_error =
        std::abs(energies.binding - energies.virial) /
        std::abs(energies.binding);

    INFO("Analytic binding energy = " << analytic_binding_energy);
    INFO("Computed binding energy = " << energies.binding);
    INFO("Computed virial energy = " << energies.virial);
    INFO("Relative virial consistency error = " << relative_consistency_error);

    constexpr double energy_tolerance      = 1.0e-5;
    constexpr double consistency_tolerance = 1.0e-6;

    CHECK_THAT(
        global_max_rel_error,
        Catch::Matchers::WithinAbs(0.0, 0.1 * potential_tolerance)
    );
    CHECK_THAT(
        global_max_abs_error,
        Catch::Matchers::WithinAbs(0.0, 0.1 * potential_tolerance)
    );
    CHECK_THAT(
        relative_binding_error,
        Catch::Matchers::WithinAbs(0.0, energy_tolerance)
    );
    CHECK_THAT(
        relative_virial_error, Catch::Matchers::WithinAbs(0.0, energy_tolerance)
    );
    CHECK_THAT(
        relative_consistency_error,
        Catch::Matchers::WithinAbs(0.0, consistency_tolerance)
    );
}

TEST_CASE(
    "Parabolic Density Virial Self-Consistency",
    tags::gravity &tags::analytic_comparison &tags::initialization
) {
    auto args  = test_utils::setup_args();
    fem::FEM f = fem::setup_fem(args.mesh_file, args, 0);
    f.mapping->ResetDisplacement();
    physics::update_stiffness_matrix(f);

    const double radius = utils::RADIUS;
    const double mass   = utils::MASS;
    const double central_density =
        (15.0 * mass) / (8.0 * M_PI * std::pow(radius, 3.0));

    auto parabolic_rho = [central_density, radius](const mfem::Vector &x) {
        const double radial_coordinate = x.Norml2();
        return central_density * (1.0 - radial_coordinate * radial_coordinate /
                                            (radius * radius));
    };

    std::unique_ptr<mfem::Coefficient> rho_coeff;
    if (f.has_mapping()) {
        rho_coeff =
            std::make_unique<mapping::PhysicalPositionFunctionCoefficient>(
                *f.mapping, parabolic_rho
            );
    } else {
        rho_coeff = std::make_unique<mfem::FunctionCoefficient>(parabolic_rho);
    }

    mfem::GridFunction rho_grid(f.densityFes.get());
    rho_grid.ProjectCoefficient(*rho_coeff);
    zero_vacuum_density(f, rho_grid);
    analysis::conserve_mass(f, rho_grid, mass);

    f.com = analysis::get_com(f, rho_grid);
    f.Q   = physics::compute_quadrupole_moment_tensor(f, rho_grid, f.com);

    const auto gravity_solution = physics::grav_potential(f, args, rho_grid);
    const int quadrature_order  = get_gravity_quadrature_order(f);
    const GravitationalEnergies energies = compute_gravitational_energies(
        f, rho_grid, gravity_solution, quadrature_order
    );
    const double analytic_binding_energy =
        -(5.0 / 7.0) * utils::G * mass * mass / radius;
    const double relative_binding_error =
        std::abs(energies.binding - analytic_binding_energy) /
        std::abs(analytic_binding_energy);
    const double relative_virial_error =
        std::abs(energies.virial - analytic_binding_energy) /
        std::abs(analytic_binding_energy);
    const double relative_consistency_error =
        std::abs(energies.binding - energies.virial) /
        std::abs(energies.binding);

    INFO("Analytic binding energy = " << analytic_binding_energy);
    INFO("Computed binding energy = " << energies.binding);
    INFO("Computed virial energy = " << energies.virial);
    INFO("Relative virial consistency error = " << relative_consistency_error);

    constexpr double analytic_tolerance    = 1.0e-5;
    constexpr double consistency_tolerance = 1.0e-6;

    CHECK_THAT(
        relative_binding_error,
        Catch::Matchers::WithinAbs(0.0, analytic_tolerance)
    );
    CHECK_THAT(
        relative_virial_error,
        Catch::Matchers::WithinAbs(0.0, analytic_tolerance)
    );
    CHECK_THAT(
        relative_consistency_error,
        Catch::Matchers::WithinAbs(0.0, consistency_tolerance)
    );
}

TEST_CASE(
    "Rational Density Virial Self-Consistency",
    tags::gravity &tags::self_consistency &tags::initialization
) {
    auto args  = test_utils::setup_args();
    fem::FEM f = fem::setup_fem(args.mesh_file, args, 0);
    f.mapping->ResetDisplacement();
    physics::update_stiffness_matrix(f);

    const double radius            = utils::RADIUS;
    const double mass              = utils::MASS;

    // Larger values are more centrally concentrated and generally harder for a
    // polynomial to represent. A regression would be considered if this test
    // does not pass for concentrations <= 16.0.
    constexpr double concentration = 16.0;
    const double density_scale     = mass / std::pow(radius, 3.0);

    auto rational_rho = [radius, density_scale](const mfem::Vector &x) {
        const double normalized_radius_squared = (x * x) / (radius * radius);
        if (normalized_radius_squared >= 1.0) {
            return 0.0;
        }

        const double denominator =
            1.0 + concentration * normalized_radius_squared;
        return density_scale * (1.0 - normalized_radius_squared) /
               (denominator * denominator);
    };

    std::unique_ptr<mfem::Coefficient> rho_coeff;
    if (f.has_mapping()) {
        rho_coeff =
            std::make_unique<mapping::PhysicalPositionFunctionCoefficient>(
                *f.mapping, rational_rho
            );
    } else {
        rho_coeff = std::make_unique<mfem::FunctionCoefficient>(rational_rho);
    }

    mfem::GridFunction rho_grid(f.densityFes.get());
    rho_grid.ProjectCoefficient(*rho_coeff);
    zero_vacuum_density(f, rho_grid);
    analysis::conserve_mass(f, rho_grid, mass);

    f.com = analysis::get_com(f, rho_grid);
    f.Q   = physics::compute_quadrupole_moment_tensor(f, rho_grid, f.com);

    const auto gravity_solution = physics::grav_potential(f, args, rho_grid);
    const int quadrature_order  = get_gravity_quadrature_order(f);
    const GravitationalEnergies energies = compute_gravitational_energies(
        f, rho_grid, gravity_solution, quadrature_order
    );

    REQUIRE(energies.binding < 0.0);
    REQUIRE(energies.virial < 0.0);

    const double relative_consistency_error =
        std::abs(energies.binding - energies.virial) /
        std::abs(energies.binding);
    INFO("W_bind = " << energies.binding);
    INFO("W_vir = " << energies.virial);
    INFO("Relative virial consistency error = " << relative_consistency_error);

    constexpr double virial_tolerance = 1.0e-5;
    CHECK_THAT(
        relative_consistency_error,
        Catch::Matchers::WithinAbs(0.0, virial_tolerance)
    );
}

TEST_CASE(
    "Homogeneous Ellipsoid Analytic Gravity",
    tags::gravity &tags::analytic_comparison &tags::initialization
) {
    auto args                = test_utils::setup_args();
    fem::FEM f               = fem::setup_fem(args.mesh_file, args, 0);

    const double radius      = utils::RADIUS;
    const double mass        = utils::MASS;
    constexpr double x_scale = 1.15;
    constexpr double y_scale = 0.95;
    constexpr double z_scale = 1.0 / (x_scale * y_scale);
    assert(std::abs(x_scale * y_scale * z_scale - 1.0) < 1.0e-14);

    const double semi_axis_x = x_scale * radius;
    const double semi_axis_y = y_scale * radius;
    const double semi_axis_z = z_scale * radius;

    auto affine_displacement = [](const mfem::Vector &x,
                                  mfem::Vector &displacement_value) {
        displacement_value.SetSize(3);
        displacement_value(0) = (x_scale - 1.0) * x(0);
        displacement_value(1) = (y_scale - 1.0) * x(1);
        displacement_value(2) = (z_scale - 1.0) * x(2);
    };

    mfem::VectorFunctionCoefficient displacement_coeff(3, affine_displacement);
    mfem::ParGridFunction displacement(f.displacementFes.get());
    displacement.ProjectCoefficient(displacement_coeff);
    f.mapping->SetDisplacement(displacement);
    physics::update_stiffness_matrix(f);

    const double analytic_volume =
        (4.0 / 3.0) * M_PI * semi_axis_x * semi_axis_y * semi_axis_z;
    const double density = mass / analytic_volume;

    mfem::GridFunction rho_grid(f.densityFes.get());
    rho_grid = density;
    zero_vacuum_density(f, rho_grid);

    const double projected_mass = analysis::domain_integrate_grid_function(
        f, rho_grid, utils::DOMAINS::STELLAR
    );
    const double numerical_density = density * mass / projected_mass;
    analysis::conserve_mass(f, rho_grid, mass);

    f.com = analysis::get_com(f, rho_grid);
    f.Q   = physics::compute_quadrupole_moment_tensor(f, rho_grid, f.com);

    const HomogeneousEllipsoidAnalytic analytic =
        compute_homogeneous_ellipsoid_analytic(
            semi_axis_x, semi_axis_y, semi_axis_z
        );
    const double coefficient_sum = analytic.coefficient_x +
                                   analytic.coefficient_y +
                                   analytic.coefficient_z;

    INFO("A_x = " << analytic.coefficient_x);
    INFO("A_y = " << analytic.coefficient_y);
    INFO("A_z = " << analytic.coefficient_z);
    INFO("A_x + A_y + A_z = " << coefficient_sum);
    REQUIRE_THAT(coefficient_sum, Catch::Matchers::WithinAbs(2.0, 1.0e-11));

    mfem::DenseMatrix analytic_quadrupole(3, 3);
    analytic_quadrupole = 0.0;
    analytic_quadrupole(0, 0) =
        (mass / 5.0) * (2.0 * semi_axis_x * semi_axis_x -
                        semi_axis_y * semi_axis_y - semi_axis_z * semi_axis_z);
    analytic_quadrupole(1, 1) =
        (mass / 5.0) * (2.0 * semi_axis_y * semi_axis_y -
                        semi_axis_x * semi_axis_x - semi_axis_z * semi_axis_z);
    analytic_quadrupole(2, 2) =
        (mass / 5.0) * (2.0 * semi_axis_z * semi_axis_z -
                        semi_axis_x * semi_axis_x - semi_axis_y * semi_axis_y);

    mfem::DenseMatrix quadrupole_difference(f.Q);
    quadrupole_difference -= analytic_quadrupole;
    const double relative_quadrupole_error =
        quadrupole_difference.FNorm() / analytic_quadrupole.FNorm();
    INFO("Relative quadrupole error = " << relative_quadrupole_error);

    HomogeneousEllipsoidHDivCoefficient analytic_field_coefficient(
        *f.mapping, numerical_density, analytic
    );
    mfem::ParGridFunction analytic_field_projection(f.gravityFluxFes.get());
    analytic_field_projection = 0.0;

    for (int i = 0; i < f.mesh->attributes.Size(); ++i) {
        const int attribute = f.mesh->attributes[i];

        if (attribute != 3) {
            analytic_field_projection.ProjectCoefficient(
                analytic_field_coefficient, attribute
            );
        }
    }

    const auto gravity_solution = physics::grav_potential(f, args, rho_grid);
    const int quadrature_order  = get_gravity_quadrature_order(f);
    double local_field_error_squared                   = 0.0;
    double local_field_norm_squared                    = 0.0;
    double local_projection_error_squared              = 0.0;
    double local_gravity_projection_difference_squared = 0.0;

    mfem::Vector x_physical(3);
    mfem::Vector grad_phi_element(3);
    mfem::Vector grad_phi_physical(3);
    mfem::Vector grad_phi_analytic(3);
    mfem::Vector grad_phi_difference(3);
    mfem::Vector projected_field_element(3);
    mfem::Vector projected_field_physical(3);
    mfem::Vector projected_field_difference(3);
    mfem::DenseMatrix map_jacobian(3, 3);
    mfem::Vector gravity_projection_difference(3);

    for (int elem_id = 0; elem_id < f.mesh->GetNE(); ++elem_id) {
        if (f.mesh->GetAttribute(elem_id) == 3) {
            continue;
        }

        mfem::ElementTransformation *transformation =
            f.mesh->GetElementTransformation(elem_id);
        const mfem::IntegrationRule &integration_rule = mfem::IntRules.Get(
            transformation->GetGeometryType(), quadrature_order
        );

        for (int q = 0; q < integration_rule.GetNPoints(); ++q) {
            const mfem::IntegrationPoint &integration_point =
                integration_rule.IntPoint(q);
            transformation->SetIntPoint(&integration_point);

            const double map_determinant =
                f.mapping->ComputeDetJ(*transformation, integration_point);
            MFEM_VERIFY(
                map_determinant > 0.0,
                "Domain mapping has a non-positive Jacobian determinant."
            );

            const double weight = transformation->Weight() *
                                  integration_point.weight * map_determinant;
            f.mapping->GetPhysicalPoint(
                *transformation, integration_point, x_physical
            );
            gravity_solution.gradPhi.GetVectorValue(
                elem_id, integration_point, grad_phi_element
            );
            f.mapping->ComputeJacobian(*transformation, map_jacobian);
            map_jacobian.Mult(grad_phi_element, grad_phi_physical);
            grad_phi_physical /= map_determinant;

            analytic_field_projection.GetVectorValue(
                elem_id, integration_point, projected_field_element
            );
            map_jacobian.Mult(
                projected_field_element, projected_field_physical
            );
            projected_field_physical /= map_determinant;

            grad_phi_analytic(0) = 2.0 * M_PI * utils::G * numerical_density *
                                   analytic.coefficient_x * x_physical(0);
            grad_phi_analytic(1) = 2.0 * M_PI * utils::G * numerical_density *
                                   analytic.coefficient_y * x_physical(1);
            grad_phi_analytic(2) = 2.0 * M_PI * utils::G * numerical_density *
                                   analytic.coefficient_z * x_physical(2);

            projected_field_difference = projected_field_physical;
            projected_field_difference -= grad_phi_analytic;
            local_projection_error_squared +=
                (projected_field_difference * projected_field_difference) *
                weight;

            gravity_projection_difference = grad_phi_physical;
            gravity_projection_difference -= projected_field_physical;
            local_gravity_projection_difference_squared +=
                (gravity_projection_difference *
                 gravity_projection_difference) *
                weight;

            grad_phi_difference = grad_phi_physical;
            grad_phi_difference -= grad_phi_analytic;
            local_field_error_squared +=
                (grad_phi_difference * grad_phi_difference) * weight;
            local_field_norm_squared +=
                (grad_phi_analytic * grad_phi_analytic) * weight;
        }
    }

    double global_field_error_squared                   = 0.0;
    double global_field_norm_squared                    = 0.0;
    double global_projection_error_squared              = 0.0;
    double global_gravity_projection_difference_squared = 0.0;
    MPI_Comm communicator = f.densityFes->GetComm();
    MPI_Allreduce(
        &local_field_error_squared, &global_field_error_squared, 1, MPI_DOUBLE,
        MPI_SUM, communicator
    );
    MPI_Allreduce(
        &local_field_norm_squared, &global_field_norm_squared, 1, MPI_DOUBLE,
        MPI_SUM, communicator
    );
    MPI_Allreduce(
        &local_projection_error_squared, &global_projection_error_squared, 1,
        MPI_DOUBLE, MPI_SUM, communicator
    );
    MPI_Allreduce(
        &local_gravity_projection_difference_squared,
        &global_gravity_projection_difference_squared, 1, MPI_DOUBLE, MPI_SUM,
        communicator
    );

    const double relative_field_error =
        std::sqrt(global_field_error_squared / global_field_norm_squared);
    const GravitationalEnergies energies = compute_gravitational_energies(
        f, rho_grid, gravity_solution, quadrature_order
    );
    const double analytic_binding_energy =
        -(3.0 / 10.0) * utils::G * mass * mass * analytic.energy_kernel;
    const double relative_binding_energy_error =
        std::abs(energies.binding - analytic_binding_energy) /
        std::abs(analytic_binding_energy);
    const double relative_virial_energy_error =
        std::abs(energies.virial - analytic_binding_energy) /
        std::abs(analytic_binding_energy);
    const double relative_consistency_error =
        std::abs(energies.binding - energies.virial) /
        std::abs(energies.binding);
    const double relative_projection_error =
        std::sqrt(global_projection_error_squared / global_field_norm_squared);
    const double gravity_to_projection_error_ratio =
        relative_projection_error > 0.0
            ? relative_field_error / relative_projection_error
            : std::numeric_limits<double>::infinity();
    const double relative_gravity_projection_difference = std::sqrt(
        global_gravity_projection_difference_squared / global_field_norm_squared
    );
    const double projection_gap_ratio =
        relative_gravity_projection_difference / relative_projection_error;

    INFO("Analytic binding energy = " << analytic_binding_energy);
    INFO("Computed binding energy = " << energies.binding);
    INFO("Computed virial energy = " << energies.virial);
    INFO("Relative field L2 error = " << relative_field_error);
    INFO("Relative binding energy error = " << relative_binding_energy_error);
    INFO("Relative virial energy error = " << relative_virial_energy_error);
    INFO("Relative virial consistency error = " << relative_consistency_error);
    INFO(
        "Relative gravity-to-RT-projection difference = "
        << relative_gravity_projection_difference
    );
    INFO("Gravity-to-projection gap ratio = " << projection_gap_ratio);

    INFO("Relative RT projection L2 error = " << relative_projection_error);
    INFO(
        "Gravity-to-projection error ratio = "
        << gravity_to_projection_error_ratio
    );
    REQUIRE(std::isfinite(relative_projection_error));

    constexpr double quadrupole_tolerance  = 2.0e-4;
    constexpr double field_tolerance       = 1.0e-5;
    constexpr double energy_tolerance      = 1.0e-5;
    constexpr double consistency_tolerance = 1.0e-5;

    CHECK_THAT(
        relative_quadrupole_error,
        Catch::Matchers::WithinAbs(0.0, quadrupole_tolerance)
    );
    CHECK_THAT(
        relative_field_error, Catch::Matchers::WithinAbs(0.0, field_tolerance)
    );
    CHECK_THAT(
        relative_binding_energy_error,
        Catch::Matchers::WithinAbs(0.0, energy_tolerance)
    );
    CHECK_THAT(
        relative_virial_energy_error,
        Catch::Matchers::WithinAbs(0.0, energy_tolerance)
    );
    CHECK_THAT(
        relative_consistency_error,
        Catch::Matchers::WithinAbs(0.0, consistency_tolerance)
    );
}

TEST_CASE(
    "Deformed Rational Density Virial Self-Consistency",
    tags::gravity &tags::self_consistency &tags::initialization
) {
    auto args                = test_utils::setup_args();
    fem::FEM f               = fem::setup_fem(args.mesh_file, args, 0);

    const double radius      = utils::RADIUS;
    const double mass        = utils::MASS;
    constexpr double x_scale = 1.15;
    constexpr double y_scale = 0.95;
    constexpr double z_scale = 1.0 / (x_scale * y_scale);
    assert(std::abs(x_scale * y_scale * z_scale - 1.0) < 1.0e-14);

    const double semi_axis_x = x_scale * radius;
    const double semi_axis_y = y_scale * radius;
    const double semi_axis_z = z_scale * radius;

    auto affine_displacement = [](const mfem::Vector &x,
                                  mfem::Vector &displacement_value) {
        displacement_value.SetSize(3);
        displacement_value(0) = (x_scale - 1.0) * x(0);
        displacement_value(1) = (y_scale - 1.0) * x(1);
        displacement_value(2) = (z_scale - 1.0) * x(2);
    };

    mfem::VectorFunctionCoefficient displacement_coeff(3, affine_displacement);
    mfem::ParGridFunction displacement(f.displacementFes.get());
    displacement.ProjectCoefficient(displacement_coeff);
    f.mapping->SetDisplacement(displacement);
    physics::update_stiffness_matrix(f);

    constexpr double concentration = 16.0;
    const double density_scale     = mass / std::pow(radius, 3.0);

    auto ellipsoidal_rho           = [semi_axis_x, semi_axis_y, semi_axis_z,
                                      density_scale](const mfem::Vector &x) {
        const double ellipsoidal_radius_squared =
            x(0) * x(0) / (semi_axis_x * semi_axis_x) +
            x(1) * x(1) / (semi_axis_y * semi_axis_y) +
            x(2) * x(2) / (semi_axis_z * semi_axis_z);

        if (ellipsoidal_radius_squared >= 1.0) {
            return 0.0;
        }

        const double denominator =
            1.0 + concentration * ellipsoidal_radius_squared;
        return density_scale * (1.0 - ellipsoidal_radius_squared) /
               (denominator * denominator);
    };

    std::unique_ptr<mfem::Coefficient> rho_coeff;
    if (f.has_mapping()) {
        rho_coeff =
            std::make_unique<mapping::PhysicalPositionFunctionCoefficient>(
                *f.mapping, ellipsoidal_rho
            );
    } else {
        rho_coeff =
            std::make_unique<mfem::FunctionCoefficient>(ellipsoidal_rho);
    }

    mfem::GridFunction rho_grid(f.densityFes.get());
    rho_grid.ProjectCoefficient(*rho_coeff);
    zero_vacuum_density(f, rho_grid);
    analysis::conserve_mass(f, rho_grid, mass);

    f.com = analysis::get_com(f, rho_grid);
    f.Q   = physics::compute_quadrupole_moment_tensor(f, rho_grid, f.com);

    const double normalized_quadrupole = f.Q.FNorm() / (mass * radius * radius);
    INFO("Normalized quadrupole = " << normalized_quadrupole);
    REQUIRE(normalized_quadrupole > 1.0e-3);

    const auto gravity_solution = physics::grav_potential(f, args, rho_grid);
    const int quadrature_order  = get_gravity_quadrature_order(f);
    const GravitationalEnergies energies = compute_gravitational_energies(
        f, rho_grid, gravity_solution, quadrature_order
    );

    REQUIRE(energies.binding < 0.0);
    REQUIRE(energies.virial < 0.0);

    const double relative_consistency_error =
        std::abs(energies.binding - energies.virial) /
        std::abs(energies.binding);
    INFO("W_bind = " << energies.binding);
    INFO("W_vir = " << energies.virial);
    INFO("Relative virial consistency error = " << relative_consistency_error);

    constexpr double virial_tolerance = 1.0e-5;
    CHECK_THAT(
        relative_consistency_error,
        Catch::Matchers::WithinAbs(0.0, virial_tolerance)
    );
}

TEST_CASE(
    "New Gravity Potential Matches Uniform Sphere Analytic",
    tags::gravity &tags::analytic_comparison &tags::initialization
        &tags::integration
) {
    auto args        = test_utils::setup_args();
    args.p.rtol      = 1.0e-13;
    args.p.max_iters = std::max(args.p.max_iters, 1000);

    fem::FEM f       = fem::setup_fem(args.mesh_file, args, 0);

    f.mapping->ResetDisplacement();
    physics::update_stiffness_matrix(f);

    mfem::ParGridFunction displacement(f.displacementFes.get());
    displacement                 = 0.0;

    const double radius          = utils::RADIUS;
    const double mass            = utils::MASS;
    const double analytic_volume = (4.0 / 3.0) * M_PI * std::pow(radius, 3.0);
    const double density         = mass / analytic_volume;

    mfem::GridFunction rho_uniform(f.densityFes.get());
    rho_uniform = density;
    zero_vacuum_density(f, rho_uniform);
    analysis::conserve_mass(f, rho_uniform, mass);

    f.com = analysis::get_com(f, rho_uniform);
    f.Q   = physics::compute_quadrupole_moment_tensor(f, rho_uniform, f.com);

    const physics::GravitySolution gravity_solution =
        physics::grav_potential_new(f, args, rho_uniform, displacement);

    constexpr double potential_tolerance =
        utils::APPROX_MAX_ACCEPTABLE_POTENTIAL_ERROR_SI_BURNING;
    double local_maximum_absolute_error = 0.0;
    double local_maximum_relative_error = 0.0;

    const int elements_to_test          = std::min(30, f.mesh->GetNE());

    for (int element_id = 0; element_id < elements_to_test; ++element_id) {
        if (f.mesh->GetAttribute(element_id) == 3) {
            continue;
        }

        mfem::ElementTransformation *transformation =
            f.mesh->GetElementTransformation(element_id);
        const mfem::IntegrationRule &integration_rule =
            mfem::IntRules.Get(transformation->GetGeometryType(), 2);
        const mfem::IntegrationPoint &integration_point =
            integration_rule.IntPoint(0);
        transformation->SetIntPoint(&integration_point);

        mfem::Vector physical_position;
        f.mapping->GetPhysicalPoint(
            *transformation, integration_point, physical_position
        );

        const double radial_coordinate = physical_position.Norml2();

        if (radial_coordinate < 1.0e-9) {
            continue;
        }

        const double analytic_potential =
            -(utils::G * mass / (2.0 * std::pow(radius, 3.0))) *
            (3.0 * radius * radius - radial_coordinate * radial_coordinate);

        const double computed_potential =
            gravity_solution.phi.GetValue(element_id, integration_point);
        const double absolute_error =
            std::abs(computed_potential - analytic_potential);
        const double relative_error =
            absolute_error / std::abs(analytic_potential);

        local_maximum_absolute_error =
            std::max(local_maximum_absolute_error, absolute_error);
        local_maximum_relative_error =
            std::max(local_maximum_relative_error, relative_error);
    }

    double global_maximum_absolute_error = 0.0;
    double global_maximum_relative_error = 0.0;

    MPI_Allreduce(
        &local_maximum_absolute_error, &global_maximum_absolute_error, 1,
        MPI_DOUBLE, MPI_MAX, f.densityFes->GetComm()
    );

    MPI_Allreduce(
        &local_maximum_relative_error, &global_maximum_relative_error, 1,
        MPI_DOUBLE, MPI_MAX, f.densityFes->GetComm()
    );

    const int quadrature_order           = get_gravity_quadrature_order(f);
    const GravitationalEnergies energies = compute_gravitational_energies(
        f, rho_uniform, gravity_solution, quadrature_order
    );
    const double analytic_binding_energy =
        -(3.0 / 5.0) * utils::G * mass * mass / radius;
    const double relative_binding_error =
        std::abs(energies.binding - analytic_binding_energy) /
        std::abs(analytic_binding_energy);
    const double relative_virial_error =
        std::abs(energies.virial - analytic_binding_energy) /
        std::abs(analytic_binding_energy);
    const double relative_consistency_error =
        std::abs(energies.binding - energies.virial) /
        std::abs(energies.binding);

    INFO(
        "Global maximum absolute potential error = "
        << global_maximum_absolute_error
    );
    INFO(
        "Global maximum relative potential error = "
        << global_maximum_relative_error
    );
    INFO("Analytic binding energy = " << analytic_binding_energy);
    INFO("New-solver binding energy = " << energies.binding);
    INFO("New-solver virial energy = " << energies.virial);
    INFO("Relative binding-energy error = " << relative_binding_error);
    INFO("Relative virial-energy error = " << relative_virial_error);
    INFO("Relative virial consistency error = " << relative_consistency_error);

    REQUIRE(energies.binding < 0.0);
    REQUIRE(energies.virial < 0.0);

    constexpr double energy_tolerance      = 1.0e-5;
    constexpr double consistency_tolerance = 1.0e-6;

    CHECK_THAT(
        global_maximum_relative_error,
        Catch::Matchers::WithinAbs(0.0, 0.1 * potential_tolerance)
    );
    CHECK_THAT(
        global_maximum_absolute_error,
        Catch::Matchers::WithinAbs(0.0, 0.1 * potential_tolerance)
    );
    CHECK_THAT(
        relative_binding_error,
        Catch::Matchers::WithinAbs(0.0, energy_tolerance)
    );
    CHECK_THAT(
        relative_virial_error, Catch::Matchers::WithinAbs(0.0, energy_tolerance)
    );
    CHECK_THAT(
        relative_consistency_error,
        Catch::Matchers::WithinAbs(0.0, consistency_tolerance)
    );
}

TEST_CASE(
    "New Gravity Potential Matches Legacy Solver On Homogeneous Ellipsoid",
    tags::gravity &tags::analytic_comparison &tags::initialization
        &tags::integration &tags::legacy_comparison
) {
    auto args                = test_utils::setup_args();
    args.p.rtol              = 1.0e-13;
    args.p.max_iters         = std::max(args.p.max_iters, 1000);

    fem::FEM f               = fem::setup_fem(args.mesh_file, args, 0);

    const double radius      = utils::RADIUS;
    const double mass        = utils::MASS;

    constexpr double x_scale = 1.15;
    constexpr double y_scale = 0.95;
    constexpr double z_scale = 1.0 / (x_scale * y_scale);

    REQUIRE_THAT(
        x_scale * y_scale * z_scale, Catch::Matchers::WithinAbs(1.0, 1.0e-14)
    );

    const double semi_axis_x = x_scale * radius;
    const double semi_axis_y = y_scale * radius;
    const double semi_axis_z = z_scale * radius;

    auto affine_displacement = [](const mfem::Vector &position,
                                  mfem::Vector &value) {
        value.SetSize(3);
        value(0) = (x_scale - 1.0) * position(0);
        value(1) = (y_scale - 1.0) * position(1);
        value(2) = (z_scale - 1.0) * position(2);
    };

    mfem::VectorFunctionCoefficient displacement_coefficient(
        3, affine_displacement
    );
    mfem::ParGridFunction displacement(f.displacementFes.get());
    displacement.ProjectCoefficient(displacement_coefficient);

    f.mapping->SetDisplacement(displacement);
    physics::update_stiffness_matrix(f);

    REQUIRE(f.gravityContext.block_A != nullptr);
    REQUIRE(f.gravityContext.source_form != nullptr);
    REQUIRE(f.domainMapperStateless != nullptr);

    const double analytic_volume =
        (4.0 / 3.0) * M_PI * semi_axis_x * semi_axis_y * semi_axis_z;
    const double density = mass / analytic_volume;

    mfem::GridFunction rho_grid(f.densityFes.get());
    rho_grid = density;
    zero_vacuum_density(f, rho_grid);

    const double projected_mass = analysis::domain_integrate_grid_function(
        f, rho_grid, utils::DOMAINS::STELLAR
    );
    const double numerical_density = density * mass / projected_mass;
    analysis::conserve_mass(f, rho_grid, mass);

    f.com = analysis::get_com(f, rho_grid);
    f.Q   = physics::compute_quadrupole_moment_tensor(f, rho_grid, f.com);

    const physics::GravitySolution legacy_solution =
        physics::grav_potential(f, args, rho_grid);
    const physics::GravitySolution new_solution =
        physics::grav_potential_new(f, args, rho_grid, displacement);

    using gravity_form = utils::blocks::gravity_field_form;

    constexpr auto gravity_gradient_residual_block =
        utils::blocks::get_residual_block<gravity_form>(
            utils::blocks::gravity_field.gradient_term
        );

    constexpr auto gravity_poisson_residual_block =
        utils::blocks::get_residual_block<gravity_form>(
            utils::blocks::gravity_field.poisson_term
        );

    const std::array<int, gravity_form::value_block_count> value_sizes{
        f.densityFes->GetTrueVSize(), f.displacementFes->GetTrueVSize(),
        f.gravityFluxFes->GetTrueVSize(), f.gravityPotentialFes->GetTrueVSize()
    };

    const std::array<int, gravity_form::residual_block_count> residual_sizes{
        f.gravityFluxFes->GetTrueVSize(), f.gravityPotentialFes->GetTrueVSize()
    };

    const utils::blocks::form_layout<gravity_form> gravity_layout(
        value_sizes, residual_sizes
    );
    const int gravity_system_size = gravity_layout.residual_offsets().Last();
    const int gravity_gradient_size =
        gravity_layout.size(gravity_gradient_residual_block);
    const int gravity_poisson_size =
        gravity_layout.size(gravity_poisson_residual_block);

    mfem::Vector density_true;
    mfem::Vector displacement_true;
    gravity_test_get_true_dofs(*f.densityFes, rho_grid, density_true);
    gravity_test_get_true_dofs(
        *f.displacementFes, displacement, displacement_true
    );

    REQUIRE(density_true.Size() == f.densityFes->GetTrueVSize());
    REQUIRE(displacement_true.Size() == f.displacementFes->GetTrueVSize());

    operators::context::gravity_field::GravityFieldLinearizationContext
        residual_linearization_context(f, *f.domainMapperStateless);

    operators::GravityFieldJacobianOperator residual_jacobian(
        f, *f.domainMapperStateless, residual_linearization_context,
        gravity_layout.value_offsets(), gravity_layout.residual_offsets()
    );

    operators::GravityFieldOperator residual_gravity_operator(
        f, *f.domainMapperStateless, residual_linearization_context,
        gravity_layout.value_offsets(), residual_jacobian
    );

    operators::context::gravity_field::GravityFieldGeometryContext
        residual_geometry_context(f, *f.domainMapperStateless);

    operators::ReducedGravityFieldOperator new_reduced_operator(
        residual_gravity_operator, residual_geometry_context, displacement_true
    );

    REQUIRE(new_reduced_operator.Width() == gravity_system_size);
    REQUIRE(new_reduced_operator.Height() == gravity_system_size);

    mfem::Vector new_right_hand_side;
    new_reduced_operator.BuildRightHandSide(density_true, new_right_hand_side);

    REQUIRE(new_right_hand_side.Size() == gravity_system_size);

    REQUIRE(f.gravityContext.source_form->Width() == density_true.Size());
    REQUIRE(f.gravityContext.source_form->Height() == gravity_poisson_size);

    mfem::Vector legacy_source_action(f.gravityContext.source_form->Height());
    legacy_source_action = 0.0;

    f.gravityContext.source_form->Mult(density_true, legacy_source_action);

    REQUIRE(legacy_source_action.Size() == gravity_poisson_size);

    mfem::BlockVector legacy_right_hand_side(gravity_layout.residual_offsets());
    legacy_right_hand_side = 0.0;
    legacy_right_hand_side.GetBlock(gravity_poisson_residual_block) =
        legacy_source_action;

    REQUIRE(legacy_right_hand_side.Size() == gravity_system_size);

    MPI_Comm communicator = f.densityFes->GetComm();

    const double new_right_hand_side_norm =
        gravity_test_global_norm(new_right_hand_side, communicator);
    const double legacy_right_hand_side_norm =
        gravity_test_global_norm(legacy_right_hand_side, communicator);

    REQUIRE(new_right_hand_side_norm > 0.0);
    REQUIRE(legacy_right_hand_side_norm > 0.0);

    mfem::Vector right_hand_side_difference(new_right_hand_side);
    right_hand_side_difference -= legacy_right_hand_side;

    const double right_hand_side_scale = std::max(
        0.5 * (new_right_hand_side_norm + legacy_right_hand_side_norm),
        std::numeric_limits<double>::min()
    );

    const double relative_right_hand_side_difference =
        gravity_test_global_norm(right_hand_side_difference, communicator) /
        right_hand_side_scale;

    mfem::BlockVector legacy_gravity_state(gravity_layout.residual_offsets());
    mfem::BlockVector new_gravity_state(gravity_layout.residual_offsets());

    legacy_gravity_state = 0.0;
    new_gravity_state    = 0.0;

    {
        mfem::Vector gradient_true;
        mfem::Vector potential_true;

        legacy_solution.gradPhi.GetTrueDofs(gradient_true);
        legacy_solution.phi.GetTrueDofs(potential_true);

        REQUIRE(gradient_true.Size() == gravity_gradient_size);
        REQUIRE(potential_true.Size() == gravity_poisson_size);

        legacy_gravity_state.GetBlock(gravity_gradient_residual_block) =
            gradient_true;
        legacy_gravity_state.GetBlock(gravity_poisson_residual_block) =
            potential_true;

        new_solution.gradPhi.GetTrueDofs(gradient_true);
        new_solution.phi.GetTrueDofs(potential_true);

        REQUIRE(gradient_true.Size() == gravity_gradient_size);
        REQUIRE(potential_true.Size() == gravity_poisson_size);

        new_gravity_state.GetBlock(gravity_gradient_residual_block) =
            gradient_true;
        new_gravity_state.GetBlock(gravity_poisson_residual_block) =
            potential_true;
    }

    REQUIRE(legacy_gravity_state.Size() == gravity_system_size);
    REQUIRE(new_gravity_state.Size() == gravity_system_size);

    REQUIRE(f.gravityContext.block_A->Width() == gravity_system_size);
    REQUIRE(f.gravityContext.block_A->Height() == gravity_system_size);

    mfem::Vector legacy_action_at_legacy_solution(gravity_system_size);
    mfem::Vector legacy_action_at_new_solution(gravity_system_size);
    mfem::Vector new_action_at_legacy_solution(gravity_system_size);
    mfem::Vector new_action_at_new_solution(gravity_system_size);

    legacy_action_at_legacy_solution = 0.0;
    legacy_action_at_new_solution    = 0.0;
    new_action_at_legacy_solution    = 0.0;
    new_action_at_new_solution       = 0.0;

    f.gravityContext.block_A->Mult(
        legacy_gravity_state, legacy_action_at_legacy_solution
    );

    f.gravityContext.block_A->Mult(
        new_gravity_state, legacy_action_at_new_solution
    );

    new_reduced_operator.Mult(
        legacy_gravity_state, new_action_at_legacy_solution
    );

    new_reduced_operator.Mult(new_gravity_state, new_action_at_new_solution);

    REQUIRE(legacy_action_at_legacy_solution.Size() == gravity_system_size);
    REQUIRE(legacy_action_at_new_solution.Size() == gravity_system_size);
    REQUIRE(new_action_at_legacy_solution.Size() == gravity_system_size);
    REQUIRE(new_action_at_new_solution.Size() == gravity_system_size);

    mfem::Vector legacy_residual_at_legacy_solution(
        legacy_action_at_legacy_solution
    );
    mfem::Vector legacy_residual_at_new_solution(legacy_action_at_new_solution);
    mfem::Vector new_residual_at_legacy_solution(new_action_at_legacy_solution);
    mfem::Vector new_residual_at_new_solution(new_action_at_new_solution);

    legacy_residual_at_legacy_solution -= legacy_right_hand_side;
    legacy_residual_at_new_solution -= legacy_right_hand_side;
    new_residual_at_legacy_solution -= new_right_hand_side;
    new_residual_at_new_solution -= new_right_hand_side;

    const GravityResidualMetrics legacy_at_legacy = measure_gravity_residual(
        legacy_residual_at_legacy_solution, gravity_layout.residual_offsets(),
        legacy_right_hand_side_norm, communicator
    );

    const GravityResidualMetrics legacy_at_new = measure_gravity_residual(
        legacy_residual_at_new_solution, gravity_layout.residual_offsets(),
        legacy_right_hand_side_norm, communicator
    );

    const GravityResidualMetrics new_at_legacy = measure_gravity_residual(
        new_residual_at_legacy_solution, gravity_layout.residual_offsets(),
        new_right_hand_side_norm, communicator
    );

    const GravityResidualMetrics new_at_new = measure_gravity_residual(
        new_residual_at_new_solution, gravity_layout.residual_offsets(),
        new_right_hand_side_norm, communicator
    );

    mfem::Vector operator_gap_at_legacy_solution(
        new_residual_at_legacy_solution
    );
    operator_gap_at_legacy_solution -= legacy_residual_at_legacy_solution;

    mfem::Vector operator_gap_at_new_solution(new_residual_at_new_solution);
    operator_gap_at_new_solution -= legacy_residual_at_new_solution;

    const GravityResidualMetrics gap_at_legacy = measure_gravity_residual(
        operator_gap_at_legacy_solution, gravity_layout.residual_offsets(),
        right_hand_side_scale, communicator
    );

    const GravityResidualMetrics gap_at_new = measure_gravity_residual(
        operator_gap_at_new_solution, gravity_layout.residual_offsets(),
        right_hand_side_scale, communicator
    );

    INFO(
        "New/legacy right-hand-side difference = "
        << relative_right_hand_side_difference
    );
    INFO(
        "Legacy operator at legacy solution: total = "
        << legacy_at_legacy.relative_total
        << ", gradient = " << legacy_at_legacy.relative_gradient
        << ", Poisson = " << legacy_at_legacy.relative_poisson
    );
    INFO(
        "Legacy operator at new solution: total = "
        << legacy_at_new.relative_total
        << ", gradient = " << legacy_at_new.relative_gradient
        << ", Poisson = " << legacy_at_new.relative_poisson
    );
    INFO(
        "New operator at legacy solution: total = "
        << new_at_legacy.relative_total
        << ", gradient = " << new_at_legacy.relative_gradient
        << ", Poisson = " << new_at_legacy.relative_poisson
    );
    INFO(
        "New operator at new solution: total = "
        << new_at_new.relative_total
        << ", gradient = " << new_at_new.relative_gradient
        << ", Poisson = " << new_at_new.relative_poisson
    );
    INFO(
        "Operator gap at legacy solution: total = "
        << gap_at_legacy.relative_total
        << ", gradient = " << gap_at_legacy.relative_gradient
        << ", Poisson = " << gap_at_legacy.relative_poisson
    );
    INFO(
        "Operator gap at new solution: total = "
        << gap_at_new.relative_total
        << ", gradient = " << gap_at_new.relative_gradient
        << ", Poisson = " << gap_at_new.relative_poisson
    );

    REQUIRE(std::isfinite(legacy_at_legacy.relative_total));
    REQUIRE(std::isfinite(legacy_at_new.relative_total));
    REQUIRE(std::isfinite(new_at_legacy.relative_total));
    REQUIRE(std::isfinite(new_at_new.relative_total));
    REQUIRE(std::isfinite(gap_at_legacy.relative_total));
    REQUIRE(std::isfinite(gap_at_new.relative_total));

    constexpr double source_parity_tolerance     = 1.0e-12;
    constexpr double diagonal_residual_tolerance = 1.0e-8;
    constexpr double poisson_gap_tolerance       = 1.0e-11;

    CHECK_THAT(
        relative_right_hand_side_difference,
        Catch::Matchers::WithinAbs(0.0, source_parity_tolerance)
    );

    CHECK_THAT(
        legacy_at_legacy.relative_total,
        Catch::Matchers::WithinAbs(0.0, diagonal_residual_tolerance)
    );

    CHECK_THAT(
        new_at_new.relative_total,
        Catch::Matchers::WithinAbs(0.0, diagonal_residual_tolerance)
    );

    CHECK_THAT(
        gap_at_legacy.relative_poisson,
        Catch::Matchers::WithinAbs(0.0, poisson_gap_tolerance)
    );

    CHECK_THAT(
        gap_at_new.relative_poisson,
        Catch::Matchers::WithinAbs(0.0, poisson_gap_tolerance)
    );

    const int quadrature_order = get_gravity_quadrature_order(f);

    const GravitySolutionComparison comparison = compare_gravity_solutions(
        f, legacy_solution, new_solution, quadrature_order
    );

    const GravitationalEnergies legacy_energies =
        compute_gravitational_energies(
            f, rho_grid, legacy_solution, quadrature_order
        );

    const GravitationalEnergies new_energies = compute_gravitational_energies(
        f, rho_grid, new_solution, quadrature_order
    );

    const HomogeneousEllipsoidAnalytic analytic =
        compute_homogeneous_ellipsoid_analytic(
            semi_axis_x, semi_axis_y, semi_axis_z
        );

    const double analytic_binding_energy =
        -(3.0 / 10.0) * utils::G * mass * mass * analytic.energy_kernel;

    const double new_binding_analytic_error =
        std::abs(new_energies.binding - analytic_binding_energy) /
        std::abs(analytic_binding_energy);

    const double new_virial_analytic_error =
        std::abs(new_energies.virial - analytic_binding_energy) /
        std::abs(analytic_binding_energy);

    const double new_virial_consistency_error =
        std::abs(new_energies.binding - new_energies.virial) /
        std::abs(new_energies.binding);

    const double relative_binding_difference =
        std::abs(new_energies.binding - legacy_energies.binding) /
        std::max(
            0.5 * (std::abs(new_energies.binding) +
                   std::abs(legacy_energies.binding)),
            std::numeric_limits<double>::min()
        );

    const double relative_virial_difference =
        std::abs(new_energies.virial - legacy_energies.virial) /
        std::max(
            0.5 * (std::abs(new_energies.virial) +
                   std::abs(legacy_energies.virial)),
            std::numeric_limits<double>::min()
        );

    INFO("Numerical density = " << numerical_density);
    INFO("Analytic binding energy = " << analytic_binding_energy);
    INFO("Legacy binding energy = " << legacy_energies.binding);
    INFO("New binding energy = " << new_energies.binding);
    INFO("Legacy virial energy = " << legacy_energies.virial);
    INFO("New virial energy = " << new_energies.virial);
    INFO("New binding analytic error = " << new_binding_analytic_error);
    INFO("New virial analytic error = " << new_virial_analytic_error);
    INFO("New virial consistency error = " << new_virial_consistency_error);
    INFO(
        "New/legacy physical gradient difference = "
        << comparison.relative_gradient_difference
    );
    INFO(
        "New/legacy physical potential difference = "
        << comparison.relative_potential_difference
    );
    INFO(
        "New/legacy binding-energy difference = " << relative_binding_difference
    );
    INFO(
        "New/legacy virial-energy difference = " << relative_virial_difference
    );

    REQUIRE(legacy_energies.binding < 0.0);
    REQUIRE(legacy_energies.virial < 0.0);
    REQUIRE(new_energies.binding < 0.0);
    REQUIRE(new_energies.virial < 0.0);

    REQUIRE(std::isfinite(comparison.relative_gradient_difference));
    REQUIRE(std::isfinite(comparison.relative_potential_difference));

    constexpr double gradient_parity_tolerance  = 5.0e-4;
    constexpr double potential_parity_tolerance = 1.0e-4;
    constexpr double energy_parity_tolerance    = 1.0e-5;
    constexpr double analytic_energy_tolerance  = 1.0e-5;
    constexpr double consistency_tolerance      = 1.0e-5;

    CHECK_THAT(
        comparison.relative_gradient_difference,
        Catch::Matchers::WithinAbs(0.0, gradient_parity_tolerance)
    );

    CHECK_THAT(
        comparison.relative_potential_difference,
        Catch::Matchers::WithinAbs(0.0, potential_parity_tolerance)
    );

    CHECK_THAT(
        relative_binding_difference,
        Catch::Matchers::WithinAbs(0.0, energy_parity_tolerance)
    );

    CHECK_THAT(
        relative_virial_difference,
        Catch::Matchers::WithinAbs(0.0, energy_parity_tolerance)
    );

    CHECK_THAT(
        new_binding_analytic_error,
        Catch::Matchers::WithinAbs(0.0, analytic_energy_tolerance)
    );

    CHECK_THAT(
        new_virial_analytic_error,
        Catch::Matchers::WithinAbs(0.0, analytic_energy_tolerance)
    );

    CHECK_THAT(
        new_virial_consistency_error,
        Catch::Matchers::WithinAbs(0.0, consistency_tolerance)
    );
}

TEST_CASE(
    "New Gravity Potential Deformed Rational Density Virial Self-Consistency",
    tags::gravity &tags::self_consistency &tags::initialization
        &tags::integration
) {
    auto args                = test_utils::setup_args();
    args.p.rtol              = 1.0e-13;
    args.p.max_iters         = std::max(args.p.max_iters, 1000);

    fem::FEM f               = fem::setup_fem(args.mesh_file, args, 0);

    const double radius      = utils::RADIUS;
    const double mass        = utils::MASS;

    constexpr double x_scale = 1.15;
    constexpr double y_scale = 0.95;
    constexpr double z_scale = 1.0 / (x_scale * y_scale);

    REQUIRE_THAT(
        x_scale * y_scale * z_scale, Catch::Matchers::WithinAbs(1.0, 1.0e-14)
    );

    const double semi_axis_x = x_scale * radius;
    const double semi_axis_y = y_scale * radius;
    const double semi_axis_z = z_scale * radius;

    auto affine_displacement = [](const mfem::Vector &position,
                                  mfem::Vector &value) {
        value.SetSize(3);
        value(0) = (x_scale - 1.0) * position(0);
        value(1) = (y_scale - 1.0) * position(1);
        value(2) = (z_scale - 1.0) * position(2);
    };

    mfem::VectorFunctionCoefficient displacement_coefficient(
        3, affine_displacement
    );

    mfem::ParGridFunction displacement(f.displacementFes.get());

    displacement.ProjectCoefficient(displacement_coefficient);

    f.mapping->SetDisplacement(displacement);
    physics::update_stiffness_matrix(f);

    constexpr double concentration = 16.0;
    const double density_scale     = mass / std::pow(radius, 3.0);

    auto ellipsoidal_density = [semi_axis_x, semi_axis_y, semi_axis_z,
                                density_scale](const mfem::Vector &position) {
        const double ellipsoidal_radius_squared =
            position(0) * position(0) / (semi_axis_x * semi_axis_x) +
            position(1) * position(1) / (semi_axis_y * semi_axis_y) +
            position(2) * position(2) / (semi_axis_z * semi_axis_z);

        if (ellipsoidal_radius_squared >= 1.0) {
            return 0.0;
        }

        const double denominator =
            1.0 + concentration * ellipsoidal_radius_squared;

        return density_scale * (1.0 - ellipsoidal_radius_squared) /
               (denominator * denominator);
    };

    mapping::PhysicalPositionFunctionCoefficient density_coefficient(
        *f.mapping, ellipsoidal_density
    );

    mfem::GridFunction density(f.densityFes.get());

    density.ProjectCoefficient(density_coefficient);

    zero_vacuum_density(f, density);
    analysis::conserve_mass(f, density, mass);

    f.com = analysis::get_com(f, density);
    f.Q   = physics::compute_quadrupole_moment_tensor(f, density, f.com);

    const double normalized_quadrupole = f.Q.FNorm() / (mass * radius * radius);

    INFO("Normalized quadrupole = " << normalized_quadrupole);

    REQUIRE(normalized_quadrupole > 1.0e-3);

    const physics::GravitySolution gravity_solution =
        physics::grav_potential_new(f, args, density, displacement);

    const int base_quadrature_order           = get_gravity_quadrature_order(f);

    const GravitationalEnergies base_energies = compute_gravitational_energies(
        f, density, gravity_solution, base_quadrature_order
    );

    const GravitationalEnergies medium_energies =
        compute_gravitational_energies(
            f, density, gravity_solution, base_quadrature_order + 4
        );

    const GravitationalEnergies fine_energies = compute_gravitational_energies(
        f, density, gravity_solution, base_quadrature_order + 8
    );

    REQUIRE(fine_energies.binding < 0.0);
    REQUIRE(fine_energies.virial < 0.0);

    const double base_consistency_error =
        std::abs(base_energies.binding - base_energies.virial) /
        std::abs(base_energies.binding);

    const double medium_consistency_error =
        std::abs(medium_energies.binding - medium_energies.virial) /
        std::abs(medium_energies.binding);

    const double fine_consistency_error =
        std::abs(fine_energies.binding - fine_energies.virial) /
        std::abs(fine_energies.binding);

    const double binding_quadrature_change =
        std::abs(fine_energies.binding - medium_energies.binding) /
        std::abs(fine_energies.binding);

    const double virial_quadrature_change =
        std::abs(fine_energies.virial - medium_energies.virial) /
        std::abs(fine_energies.virial);

    INFO("Base-order consistency error = " << base_consistency_error);

    INFO("Medium-order consistency error = " << medium_consistency_error);

    INFO("Fine-order consistency error = " << fine_consistency_error);

    INFO(
        "Medium-to-fine binding-energy change = " << binding_quadrature_change
    );

    INFO("Medium-to-fine virial-energy change = " << virial_quadrature_change);

    constexpr double virial_tolerance                = 1.0e-5;
    constexpr double diagnostic_quadrature_tolerance = 1.0e-7;

    CHECK_THAT(
        fine_consistency_error,
        Catch::Matchers::WithinAbs(0.0, virial_tolerance)
    );

    CHECK_THAT(
        binding_quadrature_change,
        Catch::Matchers::WithinAbs(0.0, diagnostic_quadrature_tolerance)
    );

    CHECK_THAT(
        virial_quadrature_change,
        Catch::Matchers::WithinAbs(0.0, diagnostic_quadrature_tolerance)
    );
}

TEST_CASE(
    "New Gravity Potential Resolves Exterior Monopole By Compactification "
    "Shell",
    tags::gravity &tags::analytic_comparison &tags::initialization
        &tags::integration
) {
    auto args        = test_utils::setup_args();
    args.p.rtol      = 1.0e-13;
    args.p.max_iters = std::max(args.p.max_iters, 1000);

    fem::FEM f       = fem::setup_fem(args.mesh_file, args, 0);

    REQUIRE(f.mapping != nullptr);
    REQUIRE(f.domainMapperStateless != nullptr);
    REQUIRE(f.compactificationCoordinate != nullptr);

    const double radius = utils::RADIUS;
    const double mass   = utils::MASS;
    const double analytic_volume =
        (4.0 / 3.0) * M_PI * radius * radius * radius;
    const double density = mass / analytic_volume;

    mfem::ParGridFunction displacement(f.displacementFes.get());
    displacement = 0.0;

    f.mapping->ResetDisplacement();
    physics::update_stiffness_matrix(f);

    mfem::GridFunction rho_uniform(f.densityFes.get());
    rho_uniform = density;

    zero_vacuum_density(f, rho_uniform);
    analysis::conserve_mass(f, rho_uniform, mass);

    f.com = analysis::get_com(f, rho_uniform);
    f.Q   = physics::compute_quadrupole_moment_tensor(f, rho_uniform, f.com);

    const physics::GravitySolution legacy_solution =
        physics::grav_potential(f, args, rho_uniform);

    const physics::GravitySolution new_solution =
        physics::grav_potential_new(f, args, rho_uniform, displacement);

    const std::array<ExteriorMonopoleShellMetrics, 5> legacy_metrics =
        measure_exterior_monopole_shells(
            f, legacy_solution, displacement, ExteriorMonopoleMapping::legacy,
            mass
        );

    const std::array<ExteriorMonopoleShellMetrics, 5> new_metrics =
        measure_exterior_monopole_shells(
            f, new_solution, displacement, ExteriorMonopoleMapping::stateless,
            mass
        );

    double maximum_new_potential_error    = 0.0;
    double maximum_new_radial_field_error = 0.0;
    double maximum_new_tangential_field   = 0.0;

    for (int shell = 0; shell < 5; ++shell) {
        DYNAMIC_SECTION(
            "Exterior coordinate in ["
            << exterior_shell_boundaries[shell] << ", "
            << exterior_shell_boundaries[shell + 1] << ")"
        ) {
            const ExteriorMonopoleShellMetrics &legacy  = legacy_metrics[shell];

            const ExteriorMonopoleShellMetrics &current = new_metrics[shell];

            maximum_new_potential_error                 = std::max(
                maximum_new_potential_error, current.potential_rms_error
            );

            maximum_new_radial_field_error = std::max(
                maximum_new_radial_field_error, current.radial_field_rms_error
            );

            maximum_new_tangential_field = std::max(
                maximum_new_tangential_field, current.tangential_field_rms
            );

            INFO("Shell = " << shell);
            INFO(
                "Exterior-coordinate interval = ["
                << exterior_shell_boundaries[shell] << ", "
                << exterior_shell_boundaries[shell + 1] << ")"
            );

            INFO(
                "Legacy physical-radius range = ["
                << legacy.minimum_radius << ", " << legacy.maximum_radius << "]"
            );

            INFO(
                "New physical-radius range = ["
                << current.minimum_radius << ", " << current.maximum_radius
                << "]"
            );

            INFO(
                "Legacy scaled-potential RMS error = "
                << legacy.potential_rms_error
            );

            INFO(
                "New scaled-potential RMS error = "
                << current.potential_rms_error
            );

            INFO(
                "Legacy scaled-radial-field RMS error = "
                << legacy.radial_field_rms_error
            );

            INFO(
                "New scaled-radial-field RMS error = "
                << current.radial_field_rms_error
            );

            INFO(
                "Legacy scaled-tangential-field RMS = "
                << legacy.tangential_field_rms
            );

            INFO(
                "New scaled-tangential-field RMS = "
                << current.tangential_field_rms
            );

            REQUIRE(legacy.quadrature_points > 0);
            REQUIRE(current.quadrature_points > 0);

            REQUIRE(std::isfinite(current.minimum_radius));
            REQUIRE(std::isfinite(current.maximum_radius));
            REQUIRE(std::isfinite(current.potential_rms_error));
            REQUIRE(std::isfinite(current.radial_field_rms_error));
            REQUIRE(std::isfinite(current.tangential_field_rms));

            CHECK(current.minimum_radius > 0.0);
            CHECK(current.maximum_radius > current.minimum_radius);

            /*
             * These are broad regression bounds, not the virial target.
             * The measured values will determine whether a monopole lift
             * is warranted.
             */
            CHECK(current.potential_rms_error < 1.0e-2);
            CHECK(current.radial_field_rms_error < 1.0e-2);
            CHECK(current.tangential_field_rms < 1.0e-2);
        }
    }

    for (int shell = 1; shell < 5; ++shell) {
        CHECK(
            new_metrics[shell].minimum_radius >=
            new_metrics[shell - 1].minimum_radius
        );

        CHECK(
            new_metrics[shell].maximum_radius >
            new_metrics[shell - 1].maximum_radius
        );
    }

    INFO(
        "Maximum new scaled-potential shell error = "
        << maximum_new_potential_error
    );

    INFO(
        "Maximum new scaled-radial-field shell error = "
        << maximum_new_radial_field_error
    );

    INFO(
        "Maximum new scaled-tangential-field shell amplitude = "
        << maximum_new_tangential_field
    );
}

TEST_CASE(
    "New Exterior Monopole Error Is Separated From Finite Element Projection "
    "Floor",
    tags::gravity &tags::analytic_comparison &tags::initialization
        &tags::integration &tags::accuracy
) {
    auto args        = test_utils::setup_args();
    args.p.rtol      = 1.0e-13;
    args.p.max_iters = std::max(args.p.max_iters, 1000);

    fem::FEM f       = fem::setup_fem(args.mesh_file, args, 0);

    REQUIRE(f.mapping != nullptr);
    REQUIRE(f.domainMapperStateless != nullptr);
    REQUIRE(f.compactificationCoordinate != nullptr);

    const double stellar_radius = utils::RADIUS;

    const double mass           = utils::MASS;

    const double analytic_volume =
        (4.0 / 3.0) * M_PI * stellar_radius * stellar_radius * stellar_radius;

    const double density = mass / analytic_volume;

    mfem::ParGridFunction displacement(f.displacementFes.get());
    displacement = 0.0;

    f.mapping->ResetDisplacement();
    physics::update_stiffness_matrix(f);

    mfem::GridFunction rho_uniform(f.densityFes.get());
    rho_uniform = density;

    zero_vacuum_density(f, rho_uniform);

    analysis::conserve_mass(f, rho_uniform, mass);

    f.com = analysis::get_com(f, rho_uniform);

    f.Q   = physics::compute_quadrupole_moment_tensor(f, rho_uniform, f.com);

    const physics::GravitySolution numerical_solution =
        physics::grav_potential_new(f, args, rho_uniform, displacement);

    StatelessMonopolePotentialCoefficient analytic_potential_coefficient(
        f, *f.domainMapperStateless, displacement, mass, stellar_radius
    );

    StatelessMonopoleHDivCoefficient analytic_field_coefficient(
        f, *f.domainMapperStateless, displacement, mass, stellar_radius
    );

    physics::GravitySolution analytic_projection(f);
    analytic_projection.phi     = 0.0;
    analytic_projection.gradPhi = 0.0;

    analytic_projection.phi.ProjectCoefficient(analytic_potential_coefficient);

    analytic_projection.gradPhi.ProjectCoefficient(analytic_field_coefficient);

    const std::array<ExteriorMonopoleShellMetrics, 5> numerical_metrics =
        measure_exterior_monopole_shells(
            f, numerical_solution, displacement,
            ExteriorMonopoleMapping::stateless, mass
        );

    const std::array<ExteriorMonopoleShellMetrics, 5> projection_metrics =
        measure_exterior_monopole_shells(
            f, analytic_projection, displacement,
            ExteriorMonopoleMapping::stateless, mass
        );

    mfem::Vector numerical_gradient_true;
    mfem::Vector numerical_potential_true;
    mfem::Vector projected_gradient_true;
    mfem::Vector projected_potential_true;

    numerical_solution.gradPhi.GetTrueDofs(numerical_gradient_true);

    numerical_solution.phi.GetTrueDofs(numerical_potential_true);

    analytic_projection.gradPhi.GetTrueDofs(projected_gradient_true);

    analytic_projection.phi.GetTrueDofs(projected_potential_true);

    MPI_Comm communicator                = f.gravityFluxFes->GetComm();

    const double gradient_projection_gap = global_relative_vector_error(
        numerical_gradient_true, projected_gradient_true, communicator
    );

    const double potential_projection_gap = global_relative_vector_error(
        numerical_potential_true, projected_potential_true, communicator
    );

    double maximum_numerical_potential_error  = 0.0;
    double maximum_projected_potential_error  = 0.0;
    double maximum_numerical_radial_error     = 0.0;
    double maximum_projected_radial_error     = 0.0;
    double maximum_numerical_tangential_field = 0.0;
    double maximum_projected_tangential_field = 0.0;

    std::ostringstream report;

    report << "Global numerical/projection gradient DOF gap = "
           << gradient_projection_gap << '\n'
           << "Global numerical/projection potential DOF gap = "
           << potential_projection_gap << '\n';

    for (int shell = 0; shell < 5; ++shell) {
        const ExteriorMonopoleShellMetrics &numerical =
            numerical_metrics[shell];

        const ExteriorMonopoleShellMetrics &projected =
            projection_metrics[shell];

        maximum_numerical_potential_error = std::max(
            maximum_numerical_potential_error, numerical.potential_rms_error
        );

        maximum_projected_potential_error = std::max(
            maximum_projected_potential_error, projected.potential_rms_error
        );

        maximum_numerical_radial_error = std::max(
            maximum_numerical_radial_error, numerical.radial_field_rms_error
        );

        maximum_projected_radial_error = std::max(
            maximum_projected_radial_error, projected.radial_field_rms_error
        );

        maximum_numerical_tangential_field = std::max(
            maximum_numerical_tangential_field, numerical.tangential_field_rms
        );

        maximum_projected_tangential_field = std::max(
            maximum_projected_tangential_field, projected.tangential_field_rms
        );

        report << "Shell " << shell << " xi=["
               << exterior_shell_boundaries[shell] << ", "
               << exterior_shell_boundaries[shell + 1] << "):\n"
               << "  radius range = [" << numerical.minimum_radius << ", "
               << numerical.maximum_radius << "]\n"
               << "  numerical potential error = "
               << numerical.potential_rms_error << '\n'
               << "  projected potential error = "
               << projected.potential_rms_error << '\n'
               << "  numerical radial-field error = "
               << numerical.radial_field_rms_error << '\n'
               << "  projected radial-field error = "
               << projected.radial_field_rms_error << '\n'
               << "  numerical tangential field = "
               << numerical.tangential_field_rms << '\n'
               << "  projected tangential field = "
               << projected.tangential_field_rms << '\n';
    }

    INFO(report.str());

    REQUIRE(std::isfinite(gradient_projection_gap));
    REQUIRE(std::isfinite(potential_projection_gap));

    REQUIRE(maximum_numerical_potential_error > 0.0);
    REQUIRE(maximum_projected_potential_error > 0.0);
    REQUIRE(maximum_numerical_radial_error > 0.0);
    REQUIRE(maximum_projected_radial_error > 0.0);

    /*
     * Broad guards against a broken projection. These are not the final
     * physical acceptance thresholds.
     */
    CHECK(maximum_projected_potential_error < 5.0e-2);
    CHECK(maximum_projected_radial_error < 5.0e-3);
    CHECK(maximum_projected_tangential_field < 5.0e-3);

    CHECK(maximum_numerical_potential_error < 5.0e-2);
    CHECK(maximum_numerical_radial_error < 5.0e-3);
    CHECK(maximum_numerical_tangential_field < 5.0e-3);

    /*
     * These are the decisive comparisons. If either fails, the solved
     * field is farther from the direct FE representation than it is from
     * the continuum monopole, indicating an operator-consistency issue
     * rather than a simple approximation floor.
     */
    CHECK(gradient_projection_gap < maximum_numerical_radial_error);

    CHECK(potential_projection_gap < maximum_numerical_potential_error);
}

TEST_CASE(
    "New Gravity Potential Matches Analytic Interior Potential For A Deformed "
    "Homogeneous Star",
    tags::gravity &tags::analytic_comparison &tags::initialization
        &tags::integration &tags::accuracy
) {
    auto args        = test_utils::setup_args();
    args.p.rtol      = 1.0e-13;
    args.p.atol      = 1.0e-14;
    args.p.max_iters = std::max(args.p.max_iters, 2000);

    fem::FEM fem     = fem::setup_fem(args.mesh_file, args, 0);

    REQUIRE(fem.mapping != nullptr);
    REQUIRE(fem.domainMapperStateless != nullptr);

    const double radius      = utils::RADIUS;
    const double mass        = utils::MASS;

    constexpr double x_scale = 1.15;
    constexpr double y_scale = 0.95;
    constexpr double z_scale = 1.0 / (x_scale * y_scale);

    const double semi_axis_x = x_scale * radius;
    const double semi_axis_y = y_scale * radius;
    const double semi_axis_z = z_scale * radius;

    REQUIRE_THAT(
        x_scale * y_scale * z_scale, Catch::Matchers::WithinAbs(1.0, 1.0e-14)
    );

    auto displacement_function = [](const mfem::Vector &position,
                                    mfem::Vector &value) {
        value.SetSize(3);
        value(0) = (x_scale - 1.0) * position(0);
        value(1) = (y_scale - 1.0) * position(1);
        value(2) = (z_scale - 1.0) * position(2);
    };

    mfem::VectorFunctionCoefficient displacement_coefficient(
        3, displacement_function
    );
    mfem::ParGridFunction displacement(fem.displacementFes.get());
    displacement.ProjectCoefficient(displacement_coefficient);

    fem.mapping->SetDisplacement(displacement);
    physics::update_stiffness_matrix(fem);

    const double analytic_volume =
        (4.0 / 3.0) * M_PI * semi_axis_x * semi_axis_y * semi_axis_z;
    const double density_value = mass / analytic_volume;

    mfem::GridFunction density(fem.densityFes.get());
    density = density_value;
    zero_vacuum_density(fem, density);

    const double projected_mass = analysis::domain_integrate_grid_function(
        fem, density, utils::DOMAINS::STELLAR
    );
    const double numerical_density = density_value * mass / projected_mass;
    analysis::conserve_mass(fem, density, mass);

    fem.com = analysis::get_com(fem, density);
    fem.Q   = physics::compute_quadrupole_moment_tensor(fem, density, fem.com);

    const HomogeneousEllipsoidAnalytic analytic =
        compute_homogeneous_ellipsoid_analytic(
            semi_axis_x, semi_axis_y, semi_axis_z
        );

    auto analytic_potential = [numerical_density, semi_axis_x, semi_axis_y,
                               semi_axis_z,
                               analytic](const mfem::Vector &position) {
        const double potential_kernel =
            semi_axis_x * semi_axis_y * semi_axis_z * analytic.energy_kernel -
            analytic.coefficient_x * position(0) * position(0) -
            analytic.coefficient_y * position(1) * position(1) -
            analytic.coefficient_z * position(2) * position(2);

        return -M_PI * utils::G * numerical_density * potential_kernel;
    };

    mapping::PhysicalPositionFunctionCoefficient analytic_potential_coefficient(
        *fem.mapping, analytic_potential
    );
    const int vacuum_attribute =
        fem.domainMapperStateless->GetVacuumElementAttribute();
    mfem::ParGridFunction projected_potential(fem.gravityPotentialFes.get());
    projected_potential.ProjectCoefficient(analytic_potential_coefficient);
    const physics::GravitySolution solution =
        physics::grav_potential_new(fem, args, density, displacement);

    const int quadrature_order            = get_gravity_quadrature_order(fem);
    double local_solution_error_squared   = 0.0;
    double local_projection_error_squared = 0.0;
    double local_solution_projection_gap_squared = 0.0;
    double local_analytic_norm_squared           = 0.0;
    double local_projected_norm_squared          = 0.0;
    double local_maximum_relative_error          = 0.0;

    mfem::Vector physical_position(3);
    mfem::DenseMatrix mapping_jacobian(3);

    for (int element_id = 0; element_id < fem.mesh->GetNE(); ++element_id) {
        mfem::ElementTransformation *transformation =
            fem.mesh->GetElementTransformation(element_id);
        if (transformation->Attribute == vacuum_attribute) {
            continue;
        }

        const mfem::IntegrationRule &rule = mfem::IntRules.Get(
            transformation->GetGeometryType(), quadrature_order
        );

        for (int quadrature_point_id = 0;
             quadrature_point_id < rule.GetNPoints(); ++quadrature_point_id) {
            const mfem::IntegrationPoint &point =
                rule.IntPoint(quadrature_point_id);
            transformation->SetIntPoint(&point);

            fem.mapping->GetPhysicalPoint(
                *transformation, point, physical_position
            );
            fem.mapping->ComputeJacobian(*transformation, mapping_jacobian);
            const double mapping_determinant = mapping_jacobian.Det();
            MFEM_VERIFY(
                mapping_determinant > 0.0,
                "Deformed potential test encountered a "
                "non-positive mapping determinant."
            );

            const double expected_potential =
                analytic_potential(physical_position);
            const double computed_potential =
                solution.phi.GetValue(element_id, point);
            const double projected_potential_value =
                projected_potential.GetValue(element_id, point);
            const double weight =
                point.weight * transformation->Weight() * mapping_determinant;

            local_solution_error_squared +=
                weight * (computed_potential - expected_potential) *
                (computed_potential - expected_potential);
            local_projection_error_squared +=
                weight * (projected_potential_value - expected_potential) *
                (projected_potential_value - expected_potential);
            local_solution_projection_gap_squared +=
                weight * (computed_potential - projected_potential_value) *
                (computed_potential - projected_potential_value);
            local_analytic_norm_squared +=
                weight * expected_potential * expected_potential;
            local_projected_norm_squared +=
                weight * projected_potential_value * projected_potential_value;
            local_maximum_relative_error = std::max(
                local_maximum_relative_error,
                std::abs(computed_potential - expected_potential) /
                    std::max(
                        std::abs(expected_potential),
                        std::numeric_limits<double>::epsilon()
                    )
            );
        }
    }

    const std::array<double, 5> local_values{
        local_solution_error_squared, local_projection_error_squared,
        local_solution_projection_gap_squared, local_analytic_norm_squared,
        local_projected_norm_squared
    };
    std::array<double, 5> global_values{};
    MPI_Allreduce(
        local_values.data(), global_values.data(),
        static_cast<int>(local_values.size()), MPI_DOUBLE, MPI_SUM,
        fem.densityFes->GetComm()
    );

    double maximum_relative_error = 0.0;
    MPI_Allreduce(
        &local_maximum_relative_error, &maximum_relative_error, 1, MPI_DOUBLE,
        MPI_MAX, fem.densityFes->GetComm()
    );

    const double solution_relative_error =
        std::sqrt(global_values[0] / global_values[3]);
    const double projection_relative_error =
        std::sqrt(global_values[1] / global_values[3]);
    const double solution_projection_gap =
        std::sqrt(global_values[2] / global_values[4]);

    INFO(
        "Ellipsoid coefficients = (" << analytic.coefficient_x << ", "
                                     << analytic.coefficient_y << ", "
                                     << analytic.coefficient_z << ")"
    );
    INFO(
        "Analytic interior-potential L2 relative error = "
        << solution_relative_error
    );
    INFO(
        "Analytic-potential FE projection L2 relative error = "
        << projection_relative_error
    );
    INFO(
        "New-solver / analytic-potential projection relative gap = "
        << solution_projection_gap
    );
    INFO(
        "Maximum interior pointwise relative potential error = "
        << maximum_relative_error
    );

    REQUIRE(std::isfinite(solution_relative_error));
    REQUIRE(std::isfinite(projection_relative_error));
    REQUIRE(std::isfinite(solution_projection_gap));
    REQUIRE(std::isfinite(maximum_relative_error));

    /*
     * The field test establishes an O(1e-3) representation floor on this
     * mesh. These are intentionally accuracy-regression guards, not claims
     * of analytic convergence of the unrefined RT/L2 representation.
     */
    CHECK(solution_relative_error < 1.0e-5);
    CHECK(projection_relative_error < 1.0e-5);
    CHECK(solution_projection_gap < 1.0e-5);
    CHECK(maximum_relative_error < 1.0e-5);
}

struct FerrersN1Analytic {
    double potential_constant;
    std::array<double, 3> first_coefficients;
    std::array<std::array<double, 3>, 3> second_coefficients;
};

class FerrersVacuumMaskedCoefficient final : public mfem::Coefficient {
public:
    FerrersVacuumMaskedCoefficient(
        Coefficient &coefficient,
        const int vacuum_attribute
    )
        : m_coefficient(coefficient),
          m_vacuum_attribute(vacuum_attribute) {
    }

    double Eval(
        mfem::ElementTransformation &transformation,
        const mfem::IntegrationPoint &integration_point
    ) override {
        if (transformation.Attribute == m_vacuum_attribute) {
            return 0.0;
        }

        return m_coefficient.Eval(transformation, integration_point);
    }

private:
    Coefficient &m_coefficient;
    int m_vacuum_attribute;
};

static double compute_ferrers_n1_coefficient(
    const std::array<
        double,
        3> &semi_axes,
    const int first_denominator_axis,
    const int second_denominator_axis
) {
    const double axis_product = semi_axes[0] * semi_axes[1] * semi_axes[2];

    const double length_scale = std::cbrt(axis_product);

    auto integrand            = [semi_axes, axis_product, length_scale,
                                 first_denominator_axis,
                                 second_denominator_axis](const double t) {
        if (t <= 0.0 || t >= 1.0) {
            return 0.0;
        }

        /*
         * Map u in [0, infinity) to t in [0, 1]:
         *
         *     u = L^2 [t / (1 - t)]^2.
         */
        const double one_minus_t = 1.0 - t;
        const double s           = t / one_minus_t;

        const double u           = length_scale * length_scale * s * s;

        const double du_dt =
            length_scale * length_scale * 2.0 * s / (one_minus_t * one_minus_t);

        const double delta = std::sqrt(
            (semi_axes[0] * semi_axes[0] + u) *
            (semi_axes[1] * semi_axes[1] + u) *
            (semi_axes[2] * semi_axes[2] + u)
        );

        double value = axis_product * du_dt / delta;

        if (first_denominator_axis >= 0) {
            value /= semi_axes[first_denominator_axis] *
                         semi_axes[first_denominator_axis] +
                     u;
        }

        if (second_denominator_axis >= 0) {
            value /= semi_axes[second_denominator_axis] *
                         semi_axes[second_denominator_axis] +
                     u;
        }

        return value;
    };

    double integration_error = 0.0;

    return boost::math::quadrature::gauss_kronrod<double, 61>::integrate(
        integrand, 0.0, 1.0, 15, 1.0e-13, &integration_error
    );
}

static FerrersN1Analytic compute_ferrers_n1_analytic(
    const double semi_axis_x,
    const double semi_axis_y,
    const double semi_axis_z
) {
    const std::array<double, 3> semi_axes{
        semi_axis_x, semi_axis_y, semi_axis_z
    };

    FerrersN1Analytic analytic{
        .potential_constant = compute_ferrers_n1_coefficient(semi_axes, -1, -1),
        .first_coefficients = {},
        .second_coefficients = {}
    };

    for (int axis = 0; axis < 3; ++axis) {
        analytic.first_coefficients[axis] =
            compute_ferrers_n1_coefficient(semi_axes, axis, -1);
    }

    for (int first_axis = 0; first_axis < 3; ++first_axis) {
        for (int second_axis = first_axis; second_axis < 3; ++second_axis) {
            const double coefficient = compute_ferrers_n1_coefficient(
                semi_axes, first_axis, second_axis
            );

            analytic.second_coefficients[first_axis][second_axis] = coefficient;

            analytic.second_coefficients[second_axis][first_axis] = coefficient;
        }
    }

    return analytic;
}

static double evaluate_ferrers_n1_potential(
    const mfem::Vector &position,
    const double central_density,
    const FerrersN1Analytic &analytic
) {
    const std::array<double, 3> coordinate_squared{
        position(0) * position(0), position(1) * position(1),
        position(2) * position(2)
    };

    /*
     * Expansion of
     *
     *  -pi G rho_c abc / 2
     *      integral [(1 - m^2(u))^2 / Delta(u)] du.
     */
    double potential_kernel = analytic.potential_constant;

    for (int first_axis = 0; first_axis < 3; ++first_axis) {
        potential_kernel -= 2.0 * analytic.first_coefficients[first_axis] *
                            coordinate_squared[first_axis];

        for (int second_axis = 0; second_axis < 3; ++second_axis) {
            potential_kernel +=
                analytic.second_coefficients[first_axis][second_axis] *
                coordinate_squared[first_axis] *
                coordinate_squared[second_axis];
        }
    }

    return -0.5 * M_PI * utils::G * central_density * potential_kernel;
}

static void evaluate_ferrers_n1_gradient(
    const mfem::Vector &position,
    const double central_density,
    const FerrersN1Analytic &analytic,
    mfem::Vector &gradient
) {
    gradient.SetSize(3);

    const std::array<double, 3> coordinate_squared{
        position(0) * position(0), position(1) * position(1),
        position(2) * position(2)
    };

    for (int axis = 0; axis < 3; ++axis) {
        double coefficient = analytic.first_coefficients[axis];

        for (int other_axis = 0; other_axis < 3; ++other_axis) {
            coefficient -= analytic.second_coefficients[axis][other_axis] *
                           coordinate_squared[other_axis];
        }

        /*
         * The solver stores grad(Phi), which points outward for a
         * negative gravitational potential.
         */
        gradient(axis) = 2.0 * M_PI * utils::G * central_density *
                         position(axis) * coefficient;
    }
}

TEST_CASE(
    "New Gravity Potential Matches Analytic Ferrers Ellipsoid",
    tags::gravity &tags::analytic_comparison &tags::initialization
        &tags::integration &tags::accuracy
) {
    auto args        = test_utils::setup_args();

    args.p.rtol      = 1.0e-13;
    args.p.atol      = 1.0e-14;
    args.p.max_iters = std::max(args.p.max_iters, 2000);

    fem::FEM fem     = fem::setup_fem(args.mesh_file, args, 0);

    REQUIRE(fem.mapping != nullptr);
    REQUIRE(fem.domainMapperStateless != nullptr);

    const double radius      = utils::RADIUS;

    const double mass        = utils::MASS;

    constexpr double x_scale = 1.0;
    constexpr double y_scale = 1.0;
    constexpr double z_scale = 1.0 / (x_scale * y_scale);

    const double semi_axis_x = x_scale * radius;
    const double semi_axis_y = y_scale * radius;
    const double semi_axis_z = z_scale * radius;

    REQUIRE_THAT(
        x_scale * y_scale * z_scale, Catch::Matchers::WithinAbs(1.0, 1.0e-14)
    );

    auto displacement_function = [](const mfem::Vector &position,
                                    mfem::Vector &value) {
        value.SetSize(3);

        value(0) = (x_scale - 1.0) * position(0);
        value(1) = (y_scale - 1.0) * position(1);
        value(2) = (z_scale - 1.0) * position(2);
    };

    mfem::VectorFunctionCoefficient displacement_coefficient(
        3, displacement_function
    );

    mfem::ParGridFunction displacement(fem.displacementFes.get());

    displacement.ProjectCoefficient(displacement_coefficient);

    fem.mapping->SetDisplacement(displacement);

    physics::update_stiffness_matrix(fem);

    const int vacuum_attribute =
        fem.domainMapperStateless->GetVacuumElementAttribute();

    /*
     * For rho = rho_c (1 - m^2), the exact mass is
     *
     *     M = 8 pi a b c rho_c / 15.
     */
    const double central_density =
        15.0 * mass / (8.0 * M_PI * semi_axis_x * semi_axis_y * semi_axis_z);

    auto density_function = [central_density, semi_axis_x, semi_axis_y,
                             semi_axis_z](const mfem::Vector &position) {
        const double ellipsoidal_radius_squared =
            position(0) * position(0) / (semi_axis_x * semi_axis_x) +
            position(1) * position(1) / (semi_axis_y * semi_axis_y) +
            position(2) * position(2) / (semi_axis_z * semi_axis_z);

        return central_density *
               std::max(0.0, 1.0 - ellipsoidal_radius_squared);
    };

    mapping::PhysicalPositionFunctionCoefficient physical_density_coefficient(
        *fem.mapping, density_function
    );

    FerrersVacuumMaskedCoefficient stellar_density_coefficient(
        physical_density_coefficient, vacuum_attribute
    );

    mfem::GridFunction density(fem.densityFes.get());

    density = 0.0;

    density.ProjectCoefficient(stellar_density_coefficient);

    const double projected_mass = analysis::domain_integrate_grid_function(
        fem, density, utils::DOMAINS::STELLAR
    );

    REQUIRE(std::isfinite(projected_mass));
    REQUIRE(projected_mass > 0.0);

    /*
     * Keep the projected source at exactly the requested mass. Because
     * projection and scaling are linear, this also gives the central
     * density appropriate to the represented source.
     */
    const double density_scale = mass / projected_mass;

    density *= density_scale;

    const double represented_central_density = central_density * density_scale;

    fem.com                                  = analysis::get_com(fem, density);

    fem.Q = physics::compute_quadrupole_moment_tensor(fem, density, fem.com);

    const double normalized_quadrupole =
        fem.Q.FNorm() / (mass * radius * radius);

    // REQUIRE(normalized_quadrupole > 1.0e-3);

    const FerrersN1Analytic analytic =
        compute_ferrers_n1_analytic(semi_axis_x, semi_axis_y, semi_axis_z);

    /*
     * Independent analytic consistency checks.
     *
     * Sum(A_i) = 2 supplies the constant part of Poisson's
     * equation. The B_ij identities supply the -m^2 part.
     */
    const double first_coefficient_sum = analytic.first_coefficients[0] +
                                         analytic.first_coefficients[1] +
                                         analytic.first_coefficients[2];

    REQUIRE_THAT(
        first_coefficient_sum, Catch::Matchers::WithinAbs(2.0, 1.0e-11)
    );

    const std::array<double, 3> semi_axes_squared{
        semi_axis_x * semi_axis_x, semi_axis_y * semi_axis_y,
        semi_axis_z * semi_axis_z
    };

    for (int axis = 0; axis < 3; ++axis) {
        double poisson_coefficient =
            3.0 * analytic.second_coefficients[axis][axis];

        for (int other_axis = 0; other_axis < 3; ++other_axis) {
            if (other_axis != axis) {
                poisson_coefficient +=
                    analytic.second_coefficients[axis][other_axis];
            }
        }

        REQUIRE_THAT(
            poisson_coefficient,
            Catch::Matchers::WithinRel(2.0 / semi_axes_squared[axis], 1.0e-10)
        );
    }

    const physics::GravitySolution solution =
        physics::grav_potential_new(fem, args, density, displacement);

    auto analytic_potential_function =
        [represented_central_density, analytic](const mfem::Vector &position) {
            return evaluate_ferrers_n1_potential(
                position, represented_central_density, analytic
            );
        };

    mapping::PhysicalPositionFunctionCoefficient physical_potential_coefficient(
        *fem.mapping, analytic_potential_function
    );

    /*
     * This wrapper is required because scalar ProjectCoefficient has no
     * attribute overload. It also prevents evaluation of the quartic
     * interior formula in the compactified vacuum.
     */
    FerrersVacuumMaskedCoefficient stellar_potential_coefficient(
        physical_potential_coefficient, vacuum_attribute
    );

    mfem::ParGridFunction projected_potential(fem.gravityPotentialFes.get());

    projected_potential = 0.0;

    projected_potential.ProjectCoefficient(stellar_potential_coefficient);

    const int quadrature_order = get_gravity_quadrature_order(fem) + 4;

    double local_potential_error_squared         = 0.0;
    double local_projection_error_squared        = 0.0;
    double local_solution_projection_gap_squared = 0.0;
    double local_potential_norm_squared          = 0.0;
    double local_projection_norm_squared         = 0.0;
    double local_field_error_squared             = 0.0;
    double local_field_norm_squared              = 0.0;
    double local_maximum_potential_error         = 0.0;

    mfem::Vector physical_position(3);
    mfem::Vector reference_field(3);
    mfem::Vector physical_field(3);
    mfem::Vector analytic_field(3);
    mfem::Vector field_difference(3);

    mfem::DenseMatrix mapping_jacobian(3, 3);

    for (int element_id = 0; element_id < fem.mesh->GetNE(); ++element_id) {
        mfem::ElementTransformation *transformation =
            fem.mesh->GetElementTransformation(element_id);

        if (transformation->Attribute == vacuum_attribute) {
            continue;
        }

        const mfem::IntegrationRule &rule = mfem::IntRules.Get(
            transformation->GetGeometryType(), quadrature_order
        );

        for (int quadrature_point_id = 0;
             quadrature_point_id < rule.GetNPoints(); ++quadrature_point_id) {
            const mfem::IntegrationPoint &point =
                rule.IntPoint(quadrature_point_id);

            transformation->SetIntPoint(&point);

            fem.mapping->GetPhysicalPoint(
                *transformation, point, physical_position
            );

            fem.mapping->ComputeJacobian(*transformation, mapping_jacobian);

            const double mapping_determinant = mapping_jacobian.Det();

            MFEM_VERIFY(
                std::isfinite(mapping_determinant) && mapping_determinant > 0.0,
                "Ferrers test encountered an invalid mapping determinant."
            );

            const double expected_potential = evaluate_ferrers_n1_potential(
                physical_position, represented_central_density, analytic
            );

            evaluate_ferrers_n1_gradient(
                physical_position, represented_central_density, analytic,
                analytic_field
            );

            const double computed_potential =
                solution.phi.GetValue(element_id, point);

            const double projected_potential_value =
                projected_potential.GetValue(element_id, point);

            solution.gradPhi.GetVectorValue(element_id, point, reference_field);

            mapping_jacobian.Mult(reference_field, physical_field);

            physical_field /= mapping_determinant;

            field_difference = physical_field;

            field_difference -= analytic_field;

            const double weight =
                point.weight * transformation->Weight() * mapping_determinant;

            const double potential_error =
                computed_potential - expected_potential;

            const double projection_error =
                projected_potential_value - expected_potential;

            const double solution_projection_difference =
                computed_potential - projected_potential_value;

            local_potential_error_squared +=
                weight * potential_error * potential_error;

            local_projection_error_squared +=
                weight * projection_error * projection_error;

            local_solution_projection_gap_squared +=
                weight * solution_projection_difference *
                solution_projection_difference;

            local_potential_norm_squared +=
                weight * expected_potential * expected_potential;

            local_projection_norm_squared +=
                weight * projected_potential_value * projected_potential_value;

            local_field_error_squared +=
                weight * (field_difference * field_difference);

            local_field_norm_squared +=
                weight * (analytic_field * analytic_field);

            local_maximum_potential_error = std::max(
                local_maximum_potential_error,
                std::abs(potential_error) /
                    std::max(
                        std::abs(expected_potential),
                        std::numeric_limits<double>::epsilon()
                    )
            );
        }
    }

    const std::array<double, 7> local_values{
        local_potential_error_squared,
        local_projection_error_squared,
        local_solution_projection_gap_squared,
        local_potential_norm_squared,
        local_projection_norm_squared,
        local_field_error_squared,
        local_field_norm_squared
    };

    std::array<double, 7> global_values{};

    MPI_Allreduce(
        local_values.data(), global_values.data(),
        static_cast<int>(local_values.size()), MPI_DOUBLE, MPI_SUM,
        fem.densityFes->GetComm()
    );

    double maximum_potential_error = 0.0;

    MPI_Allreduce(
        &local_maximum_potential_error, &maximum_potential_error, 1, MPI_DOUBLE,
        MPI_MAX, fem.densityFes->GetComm()
    );

    REQUIRE(global_values[3] > 0.0);
    REQUIRE(global_values[4] > 0.0);
    REQUIRE(global_values[6] > 0.0);

    const double potential_relative_error =
        std::sqrt(global_values[0] / global_values[3]);

    const double projection_relative_error =
        std::sqrt(global_values[1] / global_values[3]);

    const double solution_projection_gap =
        std::sqrt(global_values[2] / global_values[4]);

    const double field_relative_error =
        std::sqrt(global_values[5] / global_values[6]);

    INFO("Projected mass before normalization = " << projected_mass);

    INFO("Density normalization factor = " << density_scale);

    INFO("Normalized quadrupole = " << normalized_quadrupole);

    INFO("Ferrers potential L2 relative error = " << potential_relative_error);

    INFO(
        "Ferrers potential FE-projection relative error = "
        << projection_relative_error
    );

    INFO(
        "New-solver / Ferrers-potential projection gap = "
        << solution_projection_gap
    );

    INFO("Ferrers field L2 relative error = " << field_relative_error);

    INFO(
        "Maximum interior pointwise potential relative error = "
        << maximum_potential_error
    );

    REQUIRE(std::isfinite(potential_relative_error));
    REQUIRE(std::isfinite(projection_relative_error));
    REQUIRE(std::isfinite(solution_projection_gap));
    REQUIRE(std::isfinite(field_relative_error));
    REQUIRE(std::isfinite(maximum_potential_error));

    /*
     * Initial characterization guards. Unlike the homogeneous case,
     * this exact potential is quartic, so the FE representation floor
     * will generally be higher. Record the values before deciding
     * whether tighter regression thresholds are appropriate.
     */
    CHECK(potential_relative_error < 1.0e-5);
    CHECK(projection_relative_error < 1.0e-5);
    CHECK(solution_projection_gap < 1.0e-5);
    CHECK(field_relative_error < 1.0e-5);
    CHECK(maximum_potential_error < 1.0e-5);
}