Emily Boudreaux
85500fef3b
restricted the unknown state vector to surface deformation and implemented one prescription, NodalRadialSurface, while the full volumetric displacment field is reconstructed analytically from that. This reduced the number of degrees of freedom in the system by a factor of 80 while also removing many null vectors from the system.
606 lines
28 KiB
C++
606 lines
28 KiB
C++
#include <catch2/catch_test_macros.hpp>
|
|
|
|
#include <algorithm>
|
|
#include <array>
|
|
#include <cmath>
|
|
#include <limits>
|
|
#include <map>
|
|
#include <string>
|
|
|
|
#include <mfem.hpp>
|
|
#include <mpi.h>
|
|
|
|
import mean_field;
|
|
import test_helpers;
|
|
import experiment;
|
|
|
|
using namespace experiment;
|
|
|
|
struct AccuracyBudgetEnergies {
|
|
double binding{0.0};
|
|
double virial{0.0};
|
|
};
|
|
|
|
struct AccuracyBudgetMetrics {
|
|
double direct_relative_residual{0.0};
|
|
double gradient_relative_error{0.0};
|
|
double gradient_projection_relative_error{0.0};
|
|
double gradient_solution_projection_gap{0.0};
|
|
double potential_relative_error{0.0};
|
|
double potential_projection_relative_error{0.0};
|
|
double potential_solution_projection_gap{0.0};
|
|
double binding_relative_error{0.0};
|
|
double virial_relative_error{0.0};
|
|
double virial_consistency_error{0.0};
|
|
};
|
|
|
|
static double 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);
|
|
}
|
|
|
|
static double global_dot(
|
|
const mfem::Vector &left,
|
|
const mfem::Vector &right,
|
|
MPI_Comm communicator
|
|
) {
|
|
const double local_dot = left * right;
|
|
double global_dot_product = 0.0;
|
|
MPI_Allreduce(&local_dot, &global_dot_product, 1, MPI_DOUBLE, MPI_SUM, communicator);
|
|
return global_dot_product;
|
|
}
|
|
|
|
static void zero_vacuum_density(
|
|
const mean_field::fem::FEM &fem,
|
|
mfem::GridFunction &density
|
|
) {
|
|
using DomainSchema = mean_field::utils::domain::CoreEnvelopeVacuumDomainSchema;
|
|
|
|
const mean_field::field::FieldDofMap density_map =
|
|
mean_field::field::make_field_dof_map<mean_field::field::Density, DomainSchema>(*fem.densityFes);
|
|
|
|
mfem::Vector density_true;
|
|
density.GetTrueDofs(density_true);
|
|
|
|
const mfem::Vector supported_density = density_map.gather(density_true);
|
|
density_map.scatter(supported_density, density_true);
|
|
density.SetFromTrueDofs(density_true);
|
|
}
|
|
|
|
static int diagnostic_quadrature_order(const mean_field::fem::FEM &fem) {
|
|
return 2 * std::max(fem.gravityPotentialFes->GetMaxElementOrder(), fem.gravityFluxFes->GetMaxElementOrder()) + 8;
|
|
}
|
|
|
|
static mfem::Vector assemble_monopole_projection_rhs(
|
|
mean_field::fem::FEM &fem,
|
|
const mfem::GridFunction &displacement,
|
|
const double mass,
|
|
const double stellar_radius
|
|
) {
|
|
*fem.displacement = displacement;
|
|
|
|
mfem::Vector local_rhs(fem.gravityFluxFes->GetVSize());
|
|
local_rhs = 0.0;
|
|
|
|
const int vacuum_attribute = field_dof_test_utils::vacuum_material_attribute;
|
|
const int quadrature_order = diagnostic_quadrature_order(fem);
|
|
mean_field::mapping::GridFunctionMappingEvaluator mapping_evaluator(
|
|
*fem.domainMapperStateless, *fem.displacement, *fem.compactificationCoordinate
|
|
);
|
|
|
|
for (int element_id = 0; element_id < fem.mesh->GetNE(); ++element_id) {
|
|
const mfem::FiniteElement &gravity_element = *fem.gravityFluxFes->GetFE(element_id);
|
|
mfem::ElementTransformation *transformation = fem.mesh->GetElementTransformation(element_id);
|
|
|
|
mfem::Array<int> gravity_dofs;
|
|
mfem::DofTransformation *gravity_transform = fem.gravityFluxFes->GetElementVDofs(element_id, gravity_dofs);
|
|
|
|
const int dof_count = gravity_element.GetDof();
|
|
const int dimension = transformation->GetSpaceDim();
|
|
mfem::Vector element_rhs(dof_count);
|
|
mfem::Vector physical_position(dimension);
|
|
mfem::Vector analytic_field(dimension);
|
|
mfem::Vector pulled_field(dimension);
|
|
mfem::DenseMatrix vector_shape(dof_count, dimension);
|
|
element_rhs = 0.0;
|
|
|
|
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);
|
|
mean_field::mapping::MappingPointContext mapping_context;
|
|
MFEM_VERIFY(
|
|
mapping_evaluator.EvaluatePoint(*transformation, point, mapping_context) ==
|
|
mean_field::mapping::MappingStatus::valid,
|
|
"Invalid mapping in monopole projection RHS."
|
|
);
|
|
physical_position = mapping_context.physical_position;
|
|
|
|
const double radius = physical_position.Norml2();
|
|
MFEM_VERIFY(std::isfinite(radius) && radius > 0.0, "Invalid radius in monopole projection RHS.");
|
|
|
|
analytic_field = physical_position;
|
|
if (transformation->Attribute == vacuum_attribute) {
|
|
analytic_field *= mean_field::utils::G * mass / (radius * radius * radius);
|
|
} else {
|
|
analytic_field *= mean_field::utils::G * mass / (stellar_radius * stellar_radius * stellar_radius);
|
|
}
|
|
|
|
mapping_context.mapping_jacobian.MultTranspose(analytic_field, pulled_field);
|
|
|
|
transformation->SetIntPoint(&point);
|
|
gravity_element.CalcVShape(*transformation, vector_shape);
|
|
const double reference_weight = point.weight * transformation->Weight();
|
|
|
|
for (int dof = 0; dof < dof_count; ++dof) {
|
|
for (int component = 0; component < dimension; ++component) {
|
|
element_rhs(dof) += reference_weight * vector_shape(dof, component) * pulled_field(component);
|
|
}
|
|
}
|
|
}
|
|
|
|
if (gravity_transform != nullptr) {
|
|
gravity_transform->TransformDual(element_rhs);
|
|
}
|
|
local_rhs.AddElementVector(gravity_dofs, element_rhs);
|
|
}
|
|
|
|
mfem::Vector true_rhs(fem.gravityFluxFes->GetTrueVSize());
|
|
true_rhs = 0.0;
|
|
const mfem::Operator *prolongation = fem.gravityFluxFes->GetProlongationMatrix();
|
|
if (prolongation != nullptr) {
|
|
prolongation->MultTranspose(local_rhs, true_rhs);
|
|
} else {
|
|
true_rhs = local_rhs;
|
|
}
|
|
|
|
return true_rhs;
|
|
}
|
|
|
|
static mfem::Vector project_monopole_gradient(
|
|
mean_field::fem::FEM &fem,
|
|
const mfem::GridFunction &displacement,
|
|
const double mass,
|
|
const double stellar_radius
|
|
) {
|
|
mfem::Vector displacement_true;
|
|
displacement.GetTrueDofs(displacement_true);
|
|
|
|
const mfem::Vector projection_rhs_true = assemble_monopole_projection_rhs(fem, displacement, mass, stellar_radius);
|
|
|
|
mean_field::operators::PreparedMappedHDivMassOperator mass_operator(fem, *fem.domainMapperStateless);
|
|
mass_operator.Prepare(mass_operator.GetDisplacementMap().gather(displacement_true));
|
|
|
|
const mfem::Vector projection_rhs = mass_operator.GetFluxMap().gather(projection_rhs_true);
|
|
|
|
mfem::CGSolver solver(fem.gravityFluxFes->GetComm());
|
|
solver.SetOperator(mass_operator);
|
|
solver.SetRelTol(1.0e-11);
|
|
solver.SetAbsTol(1.0e-13);
|
|
solver.SetMaxIter(4000);
|
|
solver.SetPrintLevel(0);
|
|
|
|
mfem::Vector projected_gradient_reduced(mass_operator.GetFluxMap().reduced_size());
|
|
projected_gradient_reduced = 0.0;
|
|
solver.Mult(projection_rhs, projected_gradient_reduced);
|
|
|
|
mfem::Vector residual;
|
|
mass_operator.Mult(projected_gradient_reduced, residual);
|
|
residual -= projection_rhs;
|
|
|
|
const double relative_residual =
|
|
global_norm(residual, fem.gravityFluxFes->GetComm()) /
|
|
std::max(global_norm(projection_rhs, fem.gravityFluxFes->GetComm()), std::numeric_limits<double>::epsilon());
|
|
|
|
REQUIRE(std::isfinite(relative_residual));
|
|
REQUIRE(relative_residual < 1.0e-8);
|
|
return mass_operator.GetFluxMap().scatter(projected_gradient_reduced);
|
|
}
|
|
|
|
static double mapped_hdiv_relative_gap(
|
|
mean_field::fem::FEM &fem,
|
|
const mfem::GridFunction &displacement,
|
|
const mfem::Vector &calculated,
|
|
const mfem::Vector &reference
|
|
) {
|
|
mfem::Vector displacement_true;
|
|
displacement.GetTrueDofs(displacement_true);
|
|
|
|
mean_field::operators::PreparedMappedHDivMassOperator mass_operator(fem, *fem.domainMapperStateless);
|
|
mass_operator.Prepare(mass_operator.GetDisplacementMap().gather(displacement_true));
|
|
|
|
mfem::Vector difference(calculated);
|
|
difference -= reference;
|
|
mfem::Vector difference_action;
|
|
mfem::Vector reference_action;
|
|
const mfem::Vector reduced_difference = mass_operator.GetFluxMap().gather(difference);
|
|
const mfem::Vector reduced_reference = mass_operator.GetFluxMap().gather(reference);
|
|
mass_operator.Mult(reduced_difference, difference_action);
|
|
mass_operator.Mult(reduced_reference, reference_action);
|
|
|
|
const double difference_energy = global_dot(reduced_difference, difference_action, fem.gravityFluxFes->GetComm());
|
|
const double reference_energy = global_dot(reduced_reference, reference_action, fem.gravityFluxFes->GetComm());
|
|
MFEM_VERIFY(reference_energy > 0.0, "Projected monopole field has zero mapped H(div) norm.");
|
|
|
|
return std::sqrt(std::max(0.0, difference_energy) / reference_energy);
|
|
}
|
|
|
|
static AccuracyBudgetEnergies measure_stellar_energies(
|
|
mean_field::fem::FEM &fem,
|
|
const mfem::GridFunction &density,
|
|
const mean_field::physics::GravitySolution &solution
|
|
) {
|
|
const int vacuum_attribute = field_dof_test_utils::vacuum_material_attribute;
|
|
const int quadrature_order = diagnostic_quadrature_order(fem);
|
|
mean_field::mapping::GridFunctionMappingEvaluator mapping_evaluator(
|
|
*fem.domainMapperStateless, *fem.displacement, *fem.compactificationCoordinate
|
|
);
|
|
double local_binding = 0.0;
|
|
double local_virial = 0.0;
|
|
|
|
mfem::Vector physical_position(3);
|
|
mfem::Vector reference_field(3);
|
|
mfem::Vector physical_field(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);
|
|
|
|
mean_field::mapping::MappingPointContext mapping_context;
|
|
MFEM_VERIFY(
|
|
mapping_evaluator.EvaluatePoint(*transformation, point, mapping_context) ==
|
|
mean_field::mapping::MappingStatus::valid,
|
|
"Invalid mapping in energy diagnostic."
|
|
);
|
|
physical_position = mapping_context.physical_position;
|
|
const mfem::DenseMatrix &mapping_jacobian = mapping_context.mapping_jacobian;
|
|
const double mapping_determinant = mapping_context.mapping_determinant;
|
|
MFEM_VERIFY(mapping_determinant > 0.0, "Non-positive mapping determinant in energy diagnostic.");
|
|
|
|
solution.gradPhi.GetVectorValue(element_id, point, reference_field);
|
|
mapping_jacobian.Mult(reference_field, physical_field);
|
|
physical_field /= mapping_determinant;
|
|
|
|
const double weight = point.weight * transformation->Weight() * mapping_determinant;
|
|
const double rho = density.GetValue(element_id, point);
|
|
const double phi = solution.phi.GetValue(element_id, point);
|
|
local_binding += 0.5 * rho * phi * weight;
|
|
local_virial -= rho * (physical_position * physical_field) * weight;
|
|
}
|
|
}
|
|
|
|
AccuracyBudgetEnergies energies;
|
|
MPI_Allreduce(&local_binding, &energies.binding, 1, MPI_DOUBLE, MPI_SUM, fem.densityFes->GetComm());
|
|
MPI_Allreduce(&local_virial, &energies.virial, 1, MPI_DOUBLE, MPI_SUM, fem.densityFes->GetComm());
|
|
return energies;
|
|
}
|
|
|
|
static double reduced_gravity_relative_residual(
|
|
mean_field::fem::FEM &fem,
|
|
const mfem::GridFunction &density,
|
|
const mfem::GridFunction &displacement,
|
|
const mean_field::physics::GravitySolution &solution
|
|
) {
|
|
using GravityFieldForm = mean_field::utils::blocks::gravity_field_form;
|
|
|
|
constexpr auto gradient_block = mean_field::utils::blocks::get_residual_block<GravityFieldForm>(
|
|
mean_field::utils::blocks::gravity_field.gradient_term
|
|
);
|
|
constexpr auto poisson_block = mean_field::utils::blocks::get_residual_block<GravityFieldForm>(
|
|
mean_field::utils::blocks::gravity_field.poisson_term
|
|
);
|
|
|
|
using DomainSchema = mean_field::utils::domain::CoreEnvelopeVacuumDomainSchema;
|
|
const mean_field::field::FieldDofMap density_map =
|
|
mean_field::field::make_field_dof_map<mean_field::field::Density, DomainSchema>(*fem.densityFes);
|
|
const mean_field::field::FieldDofMap displacement_map =
|
|
mean_field::field::make_field_dof_map<mean_field::field::Displacement, DomainSchema>(*fem.displacementFes);
|
|
const mean_field::field::FieldDofMap gravity_flux_map =
|
|
mean_field::field::make_field_dof_map<mean_field::field::Gravity, DomainSchema>(*fem.gravityFluxFes);
|
|
const mean_field::field::FieldDofMap gravity_potential_map =
|
|
mean_field::field::make_field_dof_map<mean_field::field::Gravity, DomainSchema>(*fem.gravityPotentialFes);
|
|
|
|
const std::array<int, GravityFieldForm::value_block_count> value_sizes{
|
|
density_map.reduced_size(), displacement_map.reduced_size(), gravity_flux_map.reduced_size(),
|
|
gravity_potential_map.reduced_size()
|
|
};
|
|
const std::array<int, GravityFieldForm::residual_block_count> residual_sizes{
|
|
gravity_flux_map.reduced_size(), gravity_potential_map.reduced_size()
|
|
};
|
|
const mean_field::utils::blocks::form_layout<GravityFieldForm> layout(value_sizes, residual_sizes);
|
|
|
|
mfem::Vector density_true;
|
|
mfem::Vector displacement_true;
|
|
mfem::Vector gradient_true;
|
|
mfem::Vector potential_true;
|
|
density.GetTrueDofs(density_true);
|
|
displacement.GetTrueDofs(displacement_true);
|
|
solution.gradPhi.GetTrueDofs(gradient_true);
|
|
solution.phi.GetTrueDofs(potential_true);
|
|
|
|
mean_field::operators::context::gravity_field::GravityFieldLinearizationContext linearization_context(
|
|
fem, *fem.domainMapperStateless
|
|
);
|
|
mean_field::operators::GravityFieldJacobianOperator jacobian(
|
|
fem, *fem.domainMapperStateless, linearization_context, layout.value_offsets(), layout.residual_offsets()
|
|
);
|
|
mean_field::operators::GravityFieldOperator field_operator(
|
|
fem, *fem.domainMapperStateless, linearization_context, layout.value_offsets(), jacobian
|
|
);
|
|
mean_field::operators::context::gravity_field::GravityFieldGeometryContext geometry_context(
|
|
fem, *fem.domainMapperStateless
|
|
);
|
|
mean_field::operators::ReducedGravityFieldOperator reduced_operator(
|
|
field_operator, geometry_context, displacement_map.gather(displacement_true)
|
|
);
|
|
|
|
mfem::Vector right_hand_side;
|
|
reduced_operator.BuildRightHandSide(density_map.gather(density_true), right_hand_side);
|
|
|
|
mfem::BlockVector state(layout.residual_offsets());
|
|
state = 0.0;
|
|
state.GetBlock(gradient_block) = gravity_flux_map.gather(gradient_true);
|
|
state.GetBlock(poisson_block) = gravity_potential_map.gather(potential_true);
|
|
|
|
mfem::Vector residual;
|
|
reduced_operator.Mult(state, residual);
|
|
residual -= right_hand_side;
|
|
|
|
return global_norm(residual, fem.mesh->GetComm()) /
|
|
std::max(global_norm(right_hand_side, fem.mesh->GetComm()), std::numeric_limits<double>::epsilon());
|
|
}
|
|
|
|
static AccuracyBudgetMetrics measure_monopole_accuracy(
|
|
mean_field::fem::FEM &fem,
|
|
const mfem::GridFunction &density,
|
|
const mfem::GridFunction &displacement,
|
|
const mean_field::physics::GravitySolution &solution,
|
|
const mfem::ParGridFunction &projected_potential,
|
|
const mfem::Vector &projected_gradient,
|
|
const double mass,
|
|
const double stellar_radius
|
|
) {
|
|
mfem::Vector solution_gradient;
|
|
solution.gradPhi.GetTrueDofs(solution_gradient);
|
|
|
|
mfem::Vector solution_potential;
|
|
mfem::Vector projection_potential;
|
|
solution.phi.GetTrueDofs(solution_potential);
|
|
projected_potential.GetTrueDofs(projection_potential);
|
|
|
|
mfem::ParGridFunction projected_gradient_grid_function(fem.gravityFluxFes.get());
|
|
projected_gradient_grid_function.SetFromTrueDofs(projected_gradient);
|
|
|
|
double local_solution_gradient_error = 0.0;
|
|
double local_projection_gradient_error = 0.0;
|
|
double local_gradient_norm = 0.0;
|
|
double local_solution_potential_error = 0.0;
|
|
double local_projection_potential_error = 0.0;
|
|
double local_potential_norm = 0.0;
|
|
|
|
const int vacuum_attribute = field_dof_test_utils::vacuum_material_attribute;
|
|
const int quadrature_order = diagnostic_quadrature_order(fem);
|
|
mean_field::mapping::GridFunctionMappingEvaluator mapping_evaluator(
|
|
*fem.domainMapperStateless, *fem.displacement, *fem.compactificationCoordinate
|
|
);
|
|
mfem::Vector physical_position(3);
|
|
mfem::Vector analytic_gradient(3);
|
|
mfem::Vector solution_reference_gradient(3);
|
|
mfem::Vector projection_reference_gradient(3);
|
|
mfem::Vector solution_physical_gradient(3);
|
|
mfem::Vector projection_physical_gradient(3);
|
|
|
|
for (int element_id = 0; element_id < fem.mesh->GetNE(); ++element_id) {
|
|
mfem::ElementTransformation *transformation = fem.mesh->GetElementTransformation(element_id);
|
|
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);
|
|
mean_field::mapping::MappingPointContext mapping_context;
|
|
MFEM_VERIFY(
|
|
mapping_evaluator.EvaluatePoint(*transformation, point, mapping_context) ==
|
|
mean_field::mapping::MappingStatus::valid,
|
|
"Invalid mapping in accuracy diagnostic."
|
|
);
|
|
physical_position = mapping_context.physical_position;
|
|
const mfem::DenseMatrix &mapping_jacobian = mapping_context.mapping_jacobian;
|
|
const double mapping_determinant = mapping_context.mapping_determinant;
|
|
MFEM_VERIFY(mapping_determinant > 0.0, "Non-positive mapping determinant in accuracy diagnostic.");
|
|
|
|
const double radius = physical_position.Norml2();
|
|
MFEM_VERIFY(std::isfinite(radius) && radius > 0.0, "Invalid radius in monopole diagnostic.");
|
|
|
|
analytic_gradient = physical_position;
|
|
double analytic_potential = 0.0;
|
|
if (transformation->Attribute == vacuum_attribute) {
|
|
analytic_gradient *= mean_field::utils::G * mass / (radius * radius * radius);
|
|
analytic_potential = -mean_field::utils::G * mass / radius;
|
|
} else {
|
|
analytic_gradient *= mean_field::utils::G * mass / (stellar_radius * stellar_radius * stellar_radius);
|
|
analytic_potential = -mean_field::utils::G * mass *
|
|
(3.0 * stellar_radius * stellar_radius - radius * radius) /
|
|
(2.0 * stellar_radius * stellar_radius * stellar_radius);
|
|
}
|
|
|
|
solution.gradPhi.GetVectorValue(element_id, point, solution_reference_gradient);
|
|
mapping_jacobian.Mult(solution_reference_gradient, solution_physical_gradient);
|
|
solution_physical_gradient /= mapping_determinant;
|
|
|
|
projected_gradient_grid_function.GetVectorValue(element_id, point, projection_reference_gradient);
|
|
mapping_jacobian.Mult(projection_reference_gradient, projection_physical_gradient);
|
|
projection_physical_gradient /= mapping_determinant;
|
|
|
|
const double solution_potential_value = solution.phi.GetValue(element_id, point);
|
|
const double projection_potential_value = projected_potential.GetValue(element_id, point);
|
|
const double weight = point.weight * transformation->Weight() * mapping_determinant;
|
|
|
|
solution_physical_gradient -= analytic_gradient;
|
|
projection_physical_gradient -= analytic_gradient;
|
|
local_solution_gradient_error += weight * (solution_physical_gradient * solution_physical_gradient);
|
|
local_projection_gradient_error += weight * (projection_physical_gradient * projection_physical_gradient);
|
|
local_gradient_norm += weight * (analytic_gradient * analytic_gradient);
|
|
local_solution_potential_error += weight * (solution_potential_value - analytic_potential) *
|
|
(solution_potential_value - analytic_potential);
|
|
local_projection_potential_error += weight * (projection_potential_value - analytic_potential) *
|
|
(projection_potential_value - analytic_potential);
|
|
local_potential_norm += weight * analytic_potential * analytic_potential;
|
|
}
|
|
}
|
|
|
|
const std::array<double, 6> local_values{
|
|
local_solution_gradient_error, local_projection_gradient_error, local_gradient_norm,
|
|
local_solution_potential_error, local_projection_potential_error, local_potential_norm
|
|
};
|
|
std::array<double, 6> global_values{};
|
|
MPI_Allreduce(
|
|
local_values.data(), global_values.data(), static_cast<int>(local_values.size()), MPI_DOUBLE, MPI_SUM,
|
|
fem.mesh->GetComm()
|
|
);
|
|
|
|
const AccuracyBudgetEnergies energies = measure_stellar_energies(fem, density, solution);
|
|
const double analytic_energy = -3.0 * mean_field::utils::G * mass * mass / (5.0 * stellar_radius);
|
|
|
|
REQUIRE(global_values[2] > 0.0);
|
|
REQUIRE(global_values[5] > 0.0);
|
|
|
|
AccuracyBudgetMetrics metrics;
|
|
metrics.direct_relative_residual = reduced_gravity_relative_residual(fem, density, displacement, solution);
|
|
metrics.gradient_relative_error = std::sqrt(global_values[0] / global_values[2]);
|
|
metrics.gradient_projection_relative_error = std::sqrt(global_values[1] / global_values[2]);
|
|
metrics.gradient_solution_projection_gap =
|
|
mapped_hdiv_relative_gap(fem, displacement, solution_gradient, projected_gradient);
|
|
metrics.potential_relative_error = std::sqrt(global_values[3] / global_values[5]);
|
|
metrics.potential_projection_relative_error = std::sqrt(global_values[4] / global_values[5]);
|
|
mfem::Vector potential_difference(solution_potential);
|
|
potential_difference -= projection_potential;
|
|
const double projection_potential_norm = global_norm(projection_potential, fem.gravityPotentialFes->GetComm());
|
|
REQUIRE(projection_potential_norm > 0.0);
|
|
metrics.potential_solution_projection_gap =
|
|
global_norm(potential_difference, fem.gravityPotentialFes->GetComm()) / projection_potential_norm;
|
|
metrics.binding_relative_error = std::abs(energies.binding - analytic_energy) / std::abs(analytic_energy);
|
|
metrics.virial_relative_error = std::abs(energies.virial - analytic_energy) / std::abs(analytic_energy);
|
|
metrics.virial_consistency_error = std::abs(energies.binding - energies.virial) /
|
|
std::max(std::abs(energies.binding), std::numeric_limits<double>::epsilon());
|
|
return metrics;
|
|
}
|
|
|
|
static void run_monopole_case(
|
|
const std::string &sweep_name,
|
|
const std::string &case_name,
|
|
mean_field::utils::Args args,
|
|
const double solver_tolerance,
|
|
const int quadrature_boost
|
|
) {
|
|
args.p.rtol = solver_tolerance;
|
|
args.p.atol = std::min(args.p.atol, solver_tolerance * 1.0e-2);
|
|
args.p.max_iters = std::max(args.p.max_iters, 2000);
|
|
args.quadrature.global_boost = quadrature_boost;
|
|
|
|
mean_field::fem::FEM fem = mean_field::fem::setup_fem(args.mesh_file, args, 0);
|
|
REQUIRE(fem.domainMapperStateless != nullptr);
|
|
|
|
const double stellar_radius = mean_field::utils::RADIUS;
|
|
const double mass = mean_field::utils::MASS;
|
|
const double density_value = mass / ((4.0 / 3.0) * M_PI * stellar_radius * stellar_radius * stellar_radius);
|
|
|
|
mfem::ParGridFunction displacement(fem.displacementFes.get());
|
|
displacement = 0.0;
|
|
*fem.displacement = 0.0;
|
|
mfem::GridFunction density(fem.densityFes.get());
|
|
density = density_value;
|
|
zero_vacuum_density(fem, density);
|
|
mean_field::analysis::conserve_mass(fem, density, mass);
|
|
fem.com = mean_field::analysis::get_com(fem, density);
|
|
fem.Q = mean_field::physics::compute_quadrupole_moment_tensor(fem, density, fem.com);
|
|
|
|
const mean_field::physics::GravitySolution solution =
|
|
mean_field::physics::solve_gravity_field(fem, args, density, displacement);
|
|
|
|
auto analytic_potential = [mass, stellar_radius](const mfem::Vector &position) {
|
|
const double radius = position.Norml2();
|
|
if (radius >= stellar_radius) {
|
|
return -mean_field::utils::G * mass / radius;
|
|
}
|
|
return -mean_field::utils::G * mass * (3.0 * stellar_radius * stellar_radius - radius * radius) /
|
|
(2.0 * stellar_radius * stellar_radius * stellar_radius);
|
|
};
|
|
mean_field::mapping::PhysicalPositionFunctionCoefficient potential_coefficient(
|
|
*fem.domainMapperStateless, *fem.displacement, *fem.compactificationCoordinate, analytic_potential
|
|
);
|
|
mfem::ParGridFunction projected_potential(fem.gravityPotentialFes.get());
|
|
projected_potential.ProjectCoefficient(potential_coefficient);
|
|
|
|
const mfem::Vector projected_gradient = project_monopole_gradient(fem, displacement, mass, stellar_radius);
|
|
|
|
const AccuracyBudgetMetrics metrics = measure_monopole_accuracy(
|
|
fem, density, displacement, solution, projected_potential, projected_gradient, mass, stellar_radius
|
|
);
|
|
|
|
REQUIRE(std::isfinite(metrics.direct_relative_residual));
|
|
REQUIRE(std::isfinite(metrics.gradient_relative_error));
|
|
REQUIRE(std::isfinite(metrics.potential_relative_error));
|
|
REQUIRE(std::isfinite(metrics.virial_consistency_error));
|
|
|
|
record_experiment_result(
|
|
sweep_name, case_name,
|
|
{{"solver_rtol", std::to_string(solver_tolerance)},
|
|
{"quadrature_global_boost", std::to_string(quadrature_boost)},
|
|
{"mesh_file", args.mesh_file}},
|
|
{{"direct_relative_residual", metrics.direct_relative_residual},
|
|
{"gradient_relative_error", metrics.gradient_relative_error},
|
|
{"gradient_projection_relative_error", metrics.gradient_projection_relative_error},
|
|
{"gradient_solution_projection_gap", metrics.gradient_solution_projection_gap},
|
|
{"potential_relative_error", metrics.potential_relative_error},
|
|
{"potential_projection_relative_error", metrics.potential_projection_relative_error},
|
|
{"potential_solution_projection_gap", metrics.potential_solution_projection_gap},
|
|
{"binding_relative_error", metrics.binding_relative_error},
|
|
{"virial_relative_error", metrics.virial_relative_error},
|
|
{"virial_consistency_error", metrics.virial_consistency_error}}
|
|
);
|
|
}
|
|
|
|
TEST_CASE(
|
|
"Uniform Monopole Accuracy Budget: Solver Tolerance",
|
|
tags::gravity_analytic_accuracy
|
|
) {
|
|
const mean_field::utils::Args args = test_utils::setup_args();
|
|
constexpr std::array<double, 4> solver_tolerances{1.0e-8, 1.0e-10, 1.0e-12, 1.0e-14};
|
|
|
|
for (const double solver_tolerance : solver_tolerances) {
|
|
run_monopole_case("solver_tolerance", "uniform_monopole", args, solver_tolerance, 0);
|
|
}
|
|
}
|
|
|
|
TEST_CASE(
|
|
"Uniform Monopole Accuracy Budget: Quadrature",
|
|
tags::gravity_analytic_accuracy
|
|
) {
|
|
const mean_field::utils::Args args = test_utils::setup_args();
|
|
constexpr std::array<int, 3> quadrature_boosts{0, 4, 8};
|
|
|
|
for (const int quadrature_boost : quadrature_boosts) {
|
|
run_monopole_case("quadrature", "uniform_monopole", args, 1.0e-13, quadrature_boost);
|
|
}
|
|
}
|
|
|
|
TEST_CASE(
|
|
"Uniform Monopole Accuracy Budget: Projection Decomposition",
|
|
tags::gravity_analytic_accuracy
|
|
) {
|
|
run_monopole_case("projection_decomposition", "uniform_monopole", test_utils::setup_args(), 1.0e-13, 0);
|
|
}
|