Emily Boudreaux
75cc638739
This commit uses global pre allocated work space to dramatically reduce memory usage and allocation time
651 lines
30 KiB
C++
651 lines
30 KiB
C++
module;
|
|
|
|
#include <algorithm>
|
|
#include <array>
|
|
#include <bit>
|
|
#include <cmath>
|
|
#include <cstddef>
|
|
#include <cstdint>
|
|
#include <limits>
|
|
#include <mfem.hpp>
|
|
#include <mpi.h>
|
|
#include <numeric>
|
|
#include <stdexcept>
|
|
|
|
module mean_field;
|
|
|
|
namespace {
|
|
using Coefficients = std::array<double, 4>;
|
|
|
|
enum class InputFailure : int {
|
|
none,
|
|
invalid_options,
|
|
incompatible_space,
|
|
invalid_vector,
|
|
invalid_rule,
|
|
non_affine_exterior_map
|
|
};
|
|
|
|
enum class EvaluationFailure : int {
|
|
none,
|
|
invalid_accepted_mapping,
|
|
accepted_determinant_below_floor,
|
|
invalid_mapping_variation,
|
|
non_finite_polynomial
|
|
};
|
|
|
|
[[nodiscard]] bool vector_is_finite(const mfem::Vector &vector) noexcept {
|
|
for (int index = 0; index < vector.Size(); ++index) {
|
|
if (!std::isfinite(vector(index))) {
|
|
return false;
|
|
}
|
|
}
|
|
return true;
|
|
}
|
|
|
|
[[nodiscard]] bool matrix_is_finite(const mfem::DenseMatrix &matrix) noexcept {
|
|
for (int row = 0; row < matrix.Height(); ++row) {
|
|
for (int column = 0; column < matrix.Width(); ++column) {
|
|
if (!std::isfinite(matrix(row, column))) {
|
|
return false;
|
|
}
|
|
}
|
|
}
|
|
return true;
|
|
}
|
|
|
|
[[nodiscard]] bool
|
|
options_are_valid(const mean_field::deformation::LargestSafeNewtonStepSizeOptions &options) noexcept {
|
|
return std::isfinite(options.maximumStepSize) && options.maximumStepSize > 0.0 &&
|
|
std::isfinite(options.determinantFloor) && options.determinantFloor >= 0.0 &&
|
|
std::isfinite(options.fractionToBoundarySafety) && options.fractionToBoundarySafety > 0.0 &&
|
|
options.fractionToBoundarySafety < 1.0;
|
|
}
|
|
|
|
void require_mpi_success(
|
|
const int status,
|
|
const char *operation
|
|
) {
|
|
if (status != MPI_SUCCESS) {
|
|
throw std::runtime_error(operation);
|
|
}
|
|
}
|
|
|
|
[[nodiscard]] int collective_maximum(
|
|
const int localValue,
|
|
const MPI_Comm communicator,
|
|
const char *operation
|
|
) {
|
|
int globalValue = 0;
|
|
require_mpi_success(MPI_Allreduce(&localValue, &globalValue, 1, MPI_INT, MPI_MAX, communicator), operation);
|
|
return globalValue;
|
|
}
|
|
|
|
[[nodiscard]] double selected_entry(
|
|
const mfem::DenseMatrix &base,
|
|
const mfem::DenseMatrix &direction,
|
|
const unsigned int directionColumnMask,
|
|
const int row,
|
|
const int column,
|
|
const double maximumStepSize
|
|
) noexcept {
|
|
if ((directionColumnMask & (1U << static_cast<unsigned int>(column))) != 0U) {
|
|
return maximumStepSize * direction(row, column);
|
|
}
|
|
return base(row, column);
|
|
}
|
|
|
|
[[nodiscard]] double selected_column_determinant(
|
|
const mfem::DenseMatrix &base,
|
|
const mfem::DenseMatrix &direction,
|
|
const unsigned int directionColumnMask,
|
|
const int dimension,
|
|
const double maximumStepSize
|
|
) noexcept {
|
|
const auto entry = [&](const int row, const int column) {
|
|
return selected_entry(base, direction, directionColumnMask, row, column, maximumStepSize);
|
|
};
|
|
|
|
if (dimension == 1) {
|
|
return entry(0, 0);
|
|
}
|
|
if (dimension == 2) {
|
|
return entry(0, 0) * entry(1, 1) - entry(0, 1) * entry(1, 0);
|
|
}
|
|
|
|
return entry(0, 0) * (entry(1, 1) * entry(2, 2) - entry(1, 2) * entry(2, 1)) -
|
|
entry(0, 1) * (entry(1, 0) * entry(2, 2) - entry(1, 2) * entry(2, 0)) +
|
|
entry(0, 2) * (entry(1, 0) * entry(2, 1) - entry(1, 1) * entry(2, 0));
|
|
}
|
|
|
|
[[nodiscard]] Coefficients determinant_polynomial(
|
|
const mfem::DenseMatrix &base,
|
|
const mfem::DenseMatrix &direction,
|
|
const int dimension,
|
|
const double maximumStepSize,
|
|
const double determinantFloor
|
|
) noexcept {
|
|
Coefficients coefficients{};
|
|
const unsigned int termCount = 1U << static_cast<unsigned int>(dimension);
|
|
for (unsigned int mask = 0; mask < termCount; ++mask) {
|
|
const int degree = std::popcount(mask);
|
|
coefficients[static_cast<std::size_t>(degree)] +=
|
|
selected_column_determinant(base, direction, mask, dimension, maximumStepSize);
|
|
}
|
|
coefficients[0] -= determinantFloor;
|
|
return coefficients;
|
|
}
|
|
|
|
[[nodiscard]] double evaluate_polynomial(
|
|
const Coefficients &coefficients,
|
|
const double parameter
|
|
) noexcept {
|
|
return std::fma(
|
|
parameter, std::fma(parameter, std::fma(parameter, coefficients[3], coefficients[2]), coefficients[1]),
|
|
coefficients[0]
|
|
);
|
|
}
|
|
|
|
[[nodiscard]] int polynomial_degree(const Coefficients &coefficients) noexcept {
|
|
double scale = 0.0;
|
|
for (const double coefficient : coefficients) {
|
|
scale = std::max(scale, std::abs(coefficient));
|
|
}
|
|
const double tolerance = 64.0 * std::numeric_limits<double>::epsilon() * scale;
|
|
for (int degree = 3; degree > 0; --degree) {
|
|
if (std::abs(coefficients[static_cast<std::size_t>(degree)]) > tolerance) {
|
|
return degree;
|
|
}
|
|
}
|
|
return 0;
|
|
}
|
|
|
|
void append_unit_interval_root(
|
|
std::array<
|
|
double,
|
|
2> &roots,
|
|
int &rootCount,
|
|
const double root
|
|
) noexcept {
|
|
if (!std::isfinite(root) || root <= 0.0 || root >= 1.0) {
|
|
return;
|
|
}
|
|
if (rootCount > 0 && std::abs(root - roots[0]) <= 64.0 * std::numeric_limits<double>::epsilon()) {
|
|
return;
|
|
}
|
|
roots[static_cast<std::size_t>(rootCount)] = root;
|
|
++rootCount;
|
|
}
|
|
|
|
[[nodiscard]] int derivative_critical_points(
|
|
const Coefficients &coefficients,
|
|
const int degree,
|
|
std::array<
|
|
double,
|
|
2> &criticalPoints
|
|
) noexcept {
|
|
int count = 0;
|
|
if (degree == 2) {
|
|
append_unit_interval_root(criticalPoints, count, -coefficients[1] / (2.0 * coefficients[2]));
|
|
} else if (degree == 3) {
|
|
const double quadratic = 3.0 * coefficients[3];
|
|
const double linear = 2.0 * coefficients[2];
|
|
const double constant = coefficients[1];
|
|
const double discriminant = std::fma(linear, linear, -4.0 * quadratic * constant);
|
|
const double discriminantScale = linear * linear + std::abs(4.0 * quadratic * constant);
|
|
const double discriminantTolerance = 64.0 * std::numeric_limits<double>::epsilon() * discriminantScale;
|
|
|
|
if (discriminant >= -discriminantTolerance) {
|
|
const double squareRoot = std::sqrt(std::max(0.0, discriminant));
|
|
if (squareRoot == 0.0) {
|
|
append_unit_interval_root(criticalPoints, count, -linear / (2.0 * quadratic));
|
|
} else {
|
|
const double q = -0.5 * (linear + std::copysign(squareRoot, linear));
|
|
append_unit_interval_root(criticalPoints, count, q / quadratic);
|
|
append_unit_interval_root(criticalPoints, count, constant / q);
|
|
}
|
|
}
|
|
}
|
|
|
|
std::sort(criticalPoints.begin(), criticalPoints.begin() + count);
|
|
return count;
|
|
}
|
|
|
|
[[nodiscard]] double bisect_first_nonpositive_value(
|
|
const Coefficients &coefficients,
|
|
double lower,
|
|
double upper
|
|
) noexcept {
|
|
for (int iteration = 0; iteration < 80; ++iteration) {
|
|
const double middle = std::midpoint(lower, upper);
|
|
if (evaluate_polynomial(coefficients, middle) > 0.0) {
|
|
lower = middle;
|
|
} else {
|
|
upper = middle;
|
|
}
|
|
}
|
|
return upper;
|
|
}
|
|
|
|
[[nodiscard]] double first_boundary_parameter(const Coefficients &coefficients) noexcept {
|
|
const int degree = polynomial_degree(coefficients);
|
|
if (degree == 0) {
|
|
return std::numeric_limits<double>::infinity();
|
|
}
|
|
|
|
double coefficientScale = 0.0;
|
|
for (const double coefficient : coefficients) {
|
|
coefficientScale += std::abs(coefficient);
|
|
}
|
|
const double valueTolerance = 128.0 * std::numeric_limits<double>::epsilon() * coefficientScale;
|
|
|
|
std::array<double, 2> criticalPoints{};
|
|
const int criticalPointCount = derivative_critical_points(coefficients, degree, criticalPoints);
|
|
std::array<double, 4> intervalEnds{};
|
|
intervalEnds[0] = 0.0;
|
|
for (int index = 0; index < criticalPointCount; ++index) {
|
|
intervalEnds[static_cast<std::size_t>(index + 1)] = criticalPoints[static_cast<std::size_t>(index)];
|
|
}
|
|
intervalEnds[static_cast<std::size_t>(criticalPointCount + 1)] = 1.0;
|
|
|
|
for (int interval = 0; interval <= criticalPointCount; ++interval) {
|
|
const double lower = intervalEnds[static_cast<std::size_t>(interval)];
|
|
const double upper = intervalEnds[static_cast<std::size_t>(interval + 1)];
|
|
const double upperValue = evaluate_polynomial(coefficients, upper);
|
|
if (upperValue <= 0.0) {
|
|
return bisect_first_nonpositive_value(coefficients, lower, upper);
|
|
}
|
|
if (upperValue <= valueTolerance) {
|
|
// A repeated root only touches zero. Floating-point evaluation
|
|
// at the derivative root may land a few ulps above it.
|
|
return upper;
|
|
}
|
|
}
|
|
return std::numeric_limits<double>::infinity();
|
|
}
|
|
|
|
void true_to_local(
|
|
const mfem::ParFiniteElementSpace &finiteElementSpace,
|
|
const mfem::Vector &trueVector,
|
|
mfem::Vector &localVector
|
|
) {
|
|
localVector.SetSize(finiteElementSpace.GetVSize());
|
|
const mfem::Operator *prolongation = finiteElementSpace.GetProlongationMatrix();
|
|
if (prolongation != nullptr) {
|
|
prolongation->Mult(trueVector, localVector);
|
|
} else {
|
|
localVector = trueVector;
|
|
}
|
|
}
|
|
} // namespace
|
|
|
|
namespace mean_field::deformation {
|
|
void LargestSafeNewtonStepSizeOptions::Validate() const {
|
|
if (!std::isfinite(maximumStepSize) || maximumStepSize <= 0.0) {
|
|
throw std::invalid_argument("A geometry preflight requires a finite, positive maximum step size.");
|
|
}
|
|
if (!std::isfinite(determinantFloor) || determinantFloor < 0.0) {
|
|
throw std::invalid_argument("A geometry preflight requires a finite, non-negative determinant floor.");
|
|
}
|
|
if (!std::isfinite(fractionToBoundarySafety) || fractionToBoundarySafety <= 0.0 ||
|
|
fractionToBoundarySafety >= 1.0) {
|
|
throw std::invalid_argument(
|
|
"A geometry preflight requires a finite fraction-to-boundary safety factor strictly between zero "
|
|
"and one."
|
|
);
|
|
}
|
|
}
|
|
|
|
LargestSafeNewtonStepSizeEstimate estimate_largest_safe_newton_step_size(
|
|
const mapping::DomainMapper &domainMapper,
|
|
const mfem::ParFiniteElementSpace &displacementSpace,
|
|
const mfem::ParGridFunction &compactificationCoordinate,
|
|
const mfem::Vector &acceptedVolumeDisplacement,
|
|
const mfem::Vector &volumeNewtonDirection,
|
|
const std::span<const NewtonStepGeometryRule> geometryRules,
|
|
const LargestSafeNewtonStepSizeOptions &options
|
|
) {
|
|
const MPI_Comm communicator = displacementSpace.GetComm();
|
|
if (communicator == MPI_COMM_NULL) {
|
|
throw std::invalid_argument("A geometry preflight requires a valid displacement communicator.");
|
|
}
|
|
|
|
const mfem::FiniteElementSpace *compactificationSpace = compactificationCoordinate.FESpace();
|
|
mfem::Mesh *mesh = displacementSpace.GetMesh();
|
|
InputFailure localInputFailure = InputFailure::none;
|
|
const auto recordInputFailure = [&](const InputFailure failure) {
|
|
localInputFailure =
|
|
static_cast<InputFailure>(std::max(static_cast<int>(localInputFailure), static_cast<int>(failure)));
|
|
};
|
|
|
|
if (!options_are_valid(options)) {
|
|
recordInputFailure(InputFailure::invalid_options);
|
|
}
|
|
const int dimension = domainMapper.GetDimension();
|
|
const mfem::Ordering::Type ordering = displacementSpace.GetOrdering();
|
|
if (mesh == nullptr || compactificationSpace == nullptr || compactificationSpace->GetMesh() != mesh ||
|
|
dimension < 1 || dimension > 3 || (mesh != nullptr && mesh->SpaceDimension() != dimension) ||
|
|
displacementSpace.GetVDim() != dimension ||
|
|
(compactificationSpace != nullptr && compactificationSpace->GetVDim() != 1) ||
|
|
(compactificationSpace != nullptr &&
|
|
compactificationCoordinate.Size() != compactificationSpace->GetVSize()) ||
|
|
(ordering != mfem::Ordering::byNODES && ordering != mfem::Ordering::byVDIM)) {
|
|
recordInputFailure(InputFailure::incompatible_space);
|
|
}
|
|
if (acceptedVolumeDisplacement.Size() != displacementSpace.GetTrueVSize() ||
|
|
volumeNewtonDirection.Size() != displacementSpace.GetTrueVSize() ||
|
|
!vector_is_finite(acceptedVolumeDisplacement) || !vector_is_finite(volumeNewtonDirection)) {
|
|
recordInputFailure(InputFailure::invalid_vector);
|
|
}
|
|
if (geometryRules.size() > static_cast<std::size_t>(std::numeric_limits<int>::max())) {
|
|
recordInputFailure(InputFailure::invalid_rule);
|
|
}
|
|
|
|
std::uint64_t localPointCount = 0;
|
|
if (mesh != nullptr) {
|
|
for (const NewtonStepGeometryRule &entry : geometryRules) {
|
|
if (entry.element < 0 || entry.element >= mesh->GetNE() || entry.integrationRule == nullptr ||
|
|
entry.integrationRule->GetNPoints() <= 0) {
|
|
recordInputFailure(InputFailure::invalid_rule);
|
|
continue;
|
|
}
|
|
localPointCount += static_cast<std::uint64_t>(entry.integrationRule->GetNPoints());
|
|
mfem::ElementTransformation *transformation = mesh->GetElementTransformation(entry.element);
|
|
const mfem::FiniteElement *displacementElement = displacementSpace.GetFE(entry.element);
|
|
const mfem::FiniteElement *compactificationElement =
|
|
compactificationSpace != nullptr ? compactificationSpace->GetFE(entry.element) : nullptr;
|
|
if (transformation == nullptr || displacementElement == nullptr || compactificationElement == nullptr ||
|
|
transformation->GetSpaceDim() != dimension || displacementElement->GetDim() != dimension ||
|
|
compactificationElement->GetDim() != dimension ||
|
|
displacementElement->GetGeomType() != compactificationElement->GetGeomType() ||
|
|
displacementElement->GetRangeType() != mfem::FiniteElement::SCALAR ||
|
|
displacementElement->GetMapType() != mfem::FiniteElement::VALUE ||
|
|
displacementElement->GetDerivType() != mfem::FiniteElement::GRAD ||
|
|
compactificationElement->GetRangeType() != mfem::FiniteElement::SCALAR ||
|
|
compactificationElement->GetMapType() != mfem::FiniteElement::VALUE ||
|
|
compactificationElement->GetDerivType() != mfem::FiniteElement::GRAD) {
|
|
recordInputFailure(InputFailure::invalid_rule);
|
|
} else if (
|
|
domainMapper.IsCompactifiedElement(*transformation) &&
|
|
!domainMapper.GetExteriorMap().IsAffineInDisplacement()
|
|
) {
|
|
recordInputFailure(InputFailure::non_affine_exterior_map);
|
|
}
|
|
}
|
|
}
|
|
|
|
const int globalInputFailure = collective_maximum(
|
|
static_cast<int>(localInputFailure), communicator,
|
|
"The geometry preflight could not validate its distributed inputs."
|
|
);
|
|
if (globalInputFailure != static_cast<int>(InputFailure::none)) {
|
|
switch (static_cast<InputFailure>(globalInputFailure)) {
|
|
case InputFailure::invalid_options:
|
|
throw std::invalid_argument("The geometry preflight options are invalid on at least one rank.");
|
|
case InputFailure::incompatible_space:
|
|
throw std::invalid_argument(
|
|
"The geometry preflight requires compatible displacement and compactification spaces in one to "
|
|
"three dimensions."
|
|
);
|
|
case InputFailure::invalid_vector:
|
|
throw std::invalid_argument(
|
|
"The geometry preflight received an incompatible or non-finite true-DOF displacement vector."
|
|
);
|
|
case InputFailure::invalid_rule:
|
|
throw std::invalid_argument("The geometry preflight received an invalid local quadrature rule.");
|
|
case InputFailure::non_affine_exterior_map:
|
|
throw std::invalid_argument(
|
|
"The geometry preflight requires compactified mappings that are affine in displacement."
|
|
);
|
|
case InputFailure::none:
|
|
break;
|
|
}
|
|
}
|
|
|
|
const std::array<double, 3> localOptions{
|
|
options.maximumStepSize, options.determinantFloor, options.fractionToBoundarySafety
|
|
};
|
|
std::array<double, 3> minimumOptions{};
|
|
std::array<double, 3> maximumOptions{};
|
|
require_mpi_success(
|
|
MPI_Allreduce(
|
|
localOptions.data(), minimumOptions.data(), static_cast<int>(localOptions.size()), MPI_DOUBLE, MPI_MIN,
|
|
communicator
|
|
),
|
|
"The geometry preflight could not compare its distributed options."
|
|
);
|
|
require_mpi_success(
|
|
MPI_Allreduce(
|
|
localOptions.data(), maximumOptions.data(), static_cast<int>(localOptions.size()), MPI_DOUBLE, MPI_MAX,
|
|
communicator
|
|
),
|
|
"The geometry preflight could not compare its distributed options."
|
|
);
|
|
if (minimumOptions != maximumOptions) {
|
|
throw std::invalid_argument("The geometry preflight requires identical options on every rank.");
|
|
}
|
|
|
|
std::uint64_t globalPointCount = 0;
|
|
require_mpi_success(
|
|
MPI_Allreduce(&localPointCount, &globalPointCount, 1, MPI_UINT64_T, MPI_SUM, communicator),
|
|
"The geometry preflight could not count its distributed samples."
|
|
);
|
|
if (globalPointCount == 0) {
|
|
throw std::invalid_argument("The geometry preflight requires at least one quadrature point globally.");
|
|
}
|
|
|
|
mfem::Vector acceptedLocal;
|
|
mfem::Vector directionLocal;
|
|
true_to_local(displacementSpace, acceptedVolumeDisplacement, acceptedLocal);
|
|
true_to_local(displacementSpace, volumeNewtonDirection, directionLocal);
|
|
|
|
mapping::DomainMapper::Workspace workspace(dimension);
|
|
mapping::MappingPointContext mappingContext;
|
|
mapping::MappingPointVariation mappingVariation;
|
|
mfem::Array<int> displacementDofs;
|
|
mfem::Array<int> compactificationDofs;
|
|
mfem::Vector elementAcceptedDisplacement;
|
|
mfem::Vector elementDirection;
|
|
mfem::Vector elementCompactification;
|
|
|
|
double localBoundaryStep = std::numeric_limits<double>::infinity();
|
|
double localMinimumAtAccepted = std::numeric_limits<double>::infinity();
|
|
double localMinimumAtMaximum = std::numeric_limits<double>::infinity();
|
|
Coefficients localLimitingCoefficients{};
|
|
int localLimitingElement = -1;
|
|
int localLimitingRule = -1;
|
|
int localLimitingPoint = -1;
|
|
EvaluationFailure localEvaluationFailure = EvaluationFailure::none;
|
|
const auto recordEvaluationFailure = [&](const EvaluationFailure failure) {
|
|
localEvaluationFailure = static_cast<EvaluationFailure>(
|
|
std::max(static_cast<int>(localEvaluationFailure), static_cast<int>(failure))
|
|
);
|
|
};
|
|
|
|
for (std::size_t ruleIndex = 0; ruleIndex < geometryRules.size(); ++ruleIndex) {
|
|
const NewtonStepGeometryRule &entry = geometryRules[ruleIndex];
|
|
mfem::ElementTransformation *transformation = mesh->GetElementTransformation(entry.element);
|
|
mfem::DofTransformation *displacementDofTransformation =
|
|
displacementSpace.GetElementVDofs(entry.element, displacementDofs);
|
|
mfem::DofTransformation *compactificationDofTransformation =
|
|
compactificationSpace->GetElementDofs(entry.element, compactificationDofs);
|
|
|
|
acceptedLocal.GetSubVector(displacementDofs, elementAcceptedDisplacement);
|
|
directionLocal.GetSubVector(displacementDofs, elementDirection);
|
|
compactificationCoordinate.GetSubVector(compactificationDofs, elementCompactification);
|
|
if (displacementDofTransformation != nullptr) {
|
|
displacementDofTransformation->InvTransformPrimal(elementAcceptedDisplacement);
|
|
displacementDofTransformation->InvTransformPrimal(elementDirection);
|
|
}
|
|
if (compactificationDofTransformation != nullptr) {
|
|
compactificationDofTransformation->InvTransformPrimal(elementCompactification);
|
|
}
|
|
|
|
const mfem::FiniteElement &displacementElement = *displacementSpace.GetFE(entry.element);
|
|
const mfem::FiniteElement &compactificationElement = *compactificationSpace->GetFE(entry.element);
|
|
const mapping::ElementDisplacementData acceptedData(
|
|
displacementElement, elementAcceptedDisplacement, displacementSpace.GetOrdering()
|
|
);
|
|
const mapping::ElementDisplacementData directionData(
|
|
displacementElement, elementDirection, displacementSpace.GetOrdering()
|
|
);
|
|
const mapping::ElementCompactificationData compactificationData(
|
|
compactificationElement, elementCompactification
|
|
);
|
|
const mapping::ElementMappingData elementData{
|
|
.displacement = acceptedData, .compactification = compactificationData
|
|
};
|
|
|
|
for (int point = 0; point < entry.integrationRule->GetNPoints(); ++point) {
|
|
const mfem::IntegrationPoint &integrationPoint = entry.integrationRule->IntPoint(point);
|
|
const mapping::MappingStatus mappingStatus = domainMapper.EvaluatePoint(
|
|
elementData, *transformation, integrationPoint, workspace, mappingContext
|
|
);
|
|
if (mappingStatus != mapping::MappingStatus::valid) {
|
|
recordEvaluationFailure(EvaluationFailure::invalid_accepted_mapping);
|
|
continue;
|
|
}
|
|
if (mappingContext.mapping_determinant <= options.determinantFloor) {
|
|
recordEvaluationFailure(EvaluationFailure::accepted_determinant_below_floor);
|
|
continue;
|
|
}
|
|
|
|
const mapping::MappingStatus variationStatus = domainMapper.EvaluatePointVariation(
|
|
elementData, directionData, *transformation, integrationPoint, mappingContext, workspace,
|
|
mappingVariation
|
|
);
|
|
if (variationStatus != mapping::MappingStatus::valid) {
|
|
recordEvaluationFailure(EvaluationFailure::invalid_mapping_variation);
|
|
continue;
|
|
}
|
|
if (!matrix_is_finite(mappingContext.mapping_jacobian) ||
|
|
!matrix_is_finite(mappingVariation.mapping_jacobian_variation)) {
|
|
recordEvaluationFailure(EvaluationFailure::invalid_mapping_variation);
|
|
continue;
|
|
}
|
|
|
|
Coefficients coefficients = determinant_polynomial(
|
|
mappingContext.mapping_jacobian, mappingVariation.mapping_jacobian_variation, dimension,
|
|
options.maximumStepSize, options.determinantFloor
|
|
);
|
|
// Use the mapper's own determinant at the accepted state to
|
|
// avoid a second, slightly different round-off path.
|
|
coefficients[0] = mappingContext.mapping_determinant - options.determinantFloor;
|
|
const double determinantAtMaximum = evaluate_polynomial(coefficients, 1.0) + options.determinantFloor;
|
|
if (!std::isfinite(determinantAtMaximum)) {
|
|
recordEvaluationFailure(EvaluationFailure::non_finite_polynomial);
|
|
continue;
|
|
}
|
|
|
|
localMinimumAtAccepted = std::min(localMinimumAtAccepted, mappingContext.mapping_determinant);
|
|
localMinimumAtMaximum = std::min(localMinimumAtMaximum, determinantAtMaximum);
|
|
|
|
const double boundaryParameter = first_boundary_parameter(coefficients);
|
|
if (std::isfinite(boundaryParameter)) {
|
|
const double boundaryStep = options.maximumStepSize * boundaryParameter;
|
|
if (boundaryStep < localBoundaryStep) {
|
|
localBoundaryStep = boundaryStep;
|
|
localLimitingCoefficients = coefficients;
|
|
localLimitingElement = entry.element;
|
|
localLimitingRule = static_cast<int>(ruleIndex);
|
|
localLimitingPoint = point;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
const int globalEvaluationFailure = collective_maximum(
|
|
static_cast<int>(localEvaluationFailure), communicator,
|
|
"The geometry preflight could not combine its distributed mapping status."
|
|
);
|
|
if (globalEvaluationFailure != static_cast<int>(EvaluationFailure::none)) {
|
|
switch (static_cast<EvaluationFailure>(globalEvaluationFailure)) {
|
|
case EvaluationFailure::invalid_accepted_mapping:
|
|
throw std::domain_error(
|
|
"The geometry preflight received an accepted displacement with an invalid mapped geometry."
|
|
);
|
|
case EvaluationFailure::accepted_determinant_below_floor:
|
|
throw std::domain_error(
|
|
"The accepted displacement does not lie strictly above the requested determinant floor."
|
|
);
|
|
case EvaluationFailure::invalid_mapping_variation:
|
|
throw std::domain_error("The geometry preflight could not evaluate the mapping direction.");
|
|
case EvaluationFailure::non_finite_polynomial:
|
|
throw std::domain_error("The geometry preflight produced a non-finite determinant polynomial.");
|
|
case EvaluationFailure::none:
|
|
break;
|
|
}
|
|
}
|
|
|
|
double globalMinimumAtAccepted = 0.0;
|
|
double globalMinimumAtMaximum = 0.0;
|
|
require_mpi_success(
|
|
MPI_Allreduce(&localMinimumAtAccepted, &globalMinimumAtAccepted, 1, MPI_DOUBLE, MPI_MIN, communicator),
|
|
"The geometry preflight could not reduce its accepted-state determinant."
|
|
);
|
|
require_mpi_success(
|
|
MPI_Allreduce(&localMinimumAtMaximum, &globalMinimumAtMaximum, 1, MPI_DOUBLE, MPI_MIN, communicator),
|
|
"The geometry preflight could not reduce its maximum-step determinant."
|
|
);
|
|
|
|
int rank = 0;
|
|
require_mpi_success(MPI_Comm_rank(communicator, &rank), "The geometry preflight could not identify its rank.");
|
|
struct BoundaryLocation {
|
|
double step;
|
|
int rank;
|
|
};
|
|
const BoundaryLocation localLocation{.step = localBoundaryStep, .rank = rank};
|
|
BoundaryLocation globalLocation{};
|
|
require_mpi_success(
|
|
MPI_Allreduce(&localLocation, &globalLocation, 1, MPI_DOUBLE_INT, MPI_MINLOC, communicator),
|
|
"The geometry preflight could not select its limiting point."
|
|
);
|
|
|
|
const bool limitedByGeometry = std::isfinite(globalLocation.step);
|
|
std::array<int, 3> limitingLocation{-1, -1, -1};
|
|
Coefficients limitingCoefficients{};
|
|
if (limitedByGeometry) {
|
|
if (rank == globalLocation.rank) {
|
|
limitingLocation = {localLimitingElement, localLimitingRule, localLimitingPoint};
|
|
limitingCoefficients = localLimitingCoefficients;
|
|
}
|
|
require_mpi_success(
|
|
MPI_Bcast(
|
|
limitingLocation.data(), static_cast<int>(limitingLocation.size()), MPI_INT, globalLocation.rank,
|
|
communicator
|
|
),
|
|
"The geometry preflight could not broadcast its limiting location."
|
|
);
|
|
require_mpi_success(
|
|
MPI_Bcast(
|
|
limitingCoefficients.data(), static_cast<int>(limitingCoefficients.size()), MPI_DOUBLE,
|
|
globalLocation.rank, communicator
|
|
),
|
|
"The geometry preflight could not broadcast its limiting polynomial."
|
|
);
|
|
}
|
|
|
|
const double boundaryStepSize = limitedByGeometry ? globalLocation.step : options.maximumStepSize;
|
|
const double stepSize =
|
|
limitedByGeometry ? options.fractionToBoundarySafety * boundaryStepSize : options.maximumStepSize;
|
|
const double limitingPointDeterminantAtStepSize =
|
|
limitedByGeometry ? evaluate_polynomial(limitingCoefficients, stepSize / options.maximumStepSize) +
|
|
options.determinantFloor
|
|
: globalMinimumAtMaximum;
|
|
|
|
return {
|
|
.stepSize = stepSize,
|
|
.boundaryStepSize = boundaryStepSize,
|
|
.minimumDeterminantAtAcceptedState = globalMinimumAtAccepted,
|
|
.minimumDeterminantAtMaximumStepSize = globalMinimumAtMaximum,
|
|
.limitingPointDeterminantAtStepSize = limitingPointDeterminantAtStepSize,
|
|
.sampledQuadraturePointCount = globalPointCount,
|
|
.limitedByGeometry = limitedByGeometry,
|
|
.limitingRank = limitedByGeometry ? globalLocation.rank : -1,
|
|
.limitingElement = limitingLocation[0],
|
|
.limitingRule = limitingLocation[1],
|
|
.limitingQuadraturePoint = limitingLocation[2]
|
|
};
|
|
}
|
|
} // namespace mean_field::deformation
|