454 lines
26 KiB
C++
454 lines
26 KiB
C++
module;
|
|
|
|
#include <algorithm>
|
|
#include <array>
|
|
#include <cmath>
|
|
#include <concepts>
|
|
#include <cstdint>
|
|
#include <limits>
|
|
#include <memory>
|
|
#include <stdexcept>
|
|
#include <type_traits>
|
|
#include <utility>
|
|
|
|
#include <mfem.hpp>
|
|
#include <mpi.h>
|
|
|
|
export module mean_field:solver.stellar_newton_strategy;
|
|
|
|
export import :solver.newton_numerics;
|
|
export import :solver.gravity_reduction;
|
|
export import :solver.stellar_enthalpy_strategy;
|
|
|
|
namespace mean_field::solver::detail {
|
|
struct NewtonStrategyStorage {
|
|
virtual ~NewtonStrategyStorage() = default;
|
|
};
|
|
|
|
inline double StrategyDot(const mfem::Vector &a, const mfem::Vector &b, const MPI_Comm communicator) {
|
|
if (a.Size() != b.Size()) { throw std::invalid_argument("Strategy dot product has incompatible vectors."); }
|
|
double local = a * b;
|
|
double global = 0.0;
|
|
if (MPI_Allreduce(&local, &global, 1, MPI_DOUBLE, MPI_SUM, communicator) != MPI_SUCCESS) {
|
|
throw std::runtime_error("A numerical strategy could not reduce its vector product.");
|
|
}
|
|
return global;
|
|
}
|
|
|
|
inline double StrategyNorm(const mfem::Vector &vector, const MPI_Comm communicator) {
|
|
return std::sqrt(StrategyDot(vector, vector, communicator));
|
|
}
|
|
|
|
template <nonlinear::NewtonNumericalConfiguration Recipe>
|
|
void RequireCollectivelyIdenticalNumerics(const Recipe &recipe, const MPI_Comm communicator) {
|
|
using Enthalpy = typename Recipe::EnthalpyType;
|
|
using Gravity = typename Recipe::GravityType;
|
|
std::array<int, 8> integers{};
|
|
std::array<double, 2> reals{};
|
|
integers[0] = std::same_as<Enthalpy, nonlinear::PolynomialEnthalpy>;
|
|
integers[1] = std::same_as<Gravity, nonlinear::GravityCurvatureCorrection> ? 1 :
|
|
std::same_as<Gravity, nonlinear::EliminatedGravity> ? 2 : 0;
|
|
if constexpr (std::same_as<Enthalpy, nonlinear::PolynomialEnthalpy>) {
|
|
integers[2] = recipe.enthalpy.options.maximumSubdivisionDepth;
|
|
integers[3] = recipe.enthalpy.options.maximumSubcells;
|
|
}
|
|
if constexpr (!std::same_as<Gravity, nonlinear::MonolithicGravity>) {
|
|
const auto &options = recipe.gravity.gravity;
|
|
integers[4] = options.control.maximumIterations;
|
|
integers[5] = options.backend.restartLength;
|
|
integers[6] = options.backend.printLevel;
|
|
reals = {options.control.relativeTolerance, options.control.absoluteTolerance};
|
|
}
|
|
integers[7] = std::isfinite(reals[0]) && std::isfinite(reals[1]);
|
|
std::array<int, 8> minimumIntegers{}, maximumIntegers{};
|
|
std::array<double, 2> minimumReals{}, maximumReals{};
|
|
if (MPI_Allreduce(integers.data(), minimumIntegers.data(), 8, MPI_INT, MPI_MIN, communicator) != MPI_SUCCESS ||
|
|
MPI_Allreduce(integers.data(), maximumIntegers.data(), 8, MPI_INT, MPI_MAX, communicator) != MPI_SUCCESS ||
|
|
MPI_Allreduce(reals.data(), minimumReals.data(), 2, MPI_DOUBLE, MPI_MIN, communicator) != MPI_SUCCESS ||
|
|
MPI_Allreduce(reals.data(), maximumReals.data(), 2, MPI_DOUBLE, MPI_MAX, communicator) != MPI_SUCCESS) {
|
|
throw std::runtime_error("Newton could not compare numerical strategies across ranks.");
|
|
}
|
|
if (minimumIntegers != maximumIntegers || minimumReals != maximumReals ||
|
|
minimumIntegers[7] == 0) {
|
|
throw std::invalid_argument("Newton numerical strategies and controls must agree on every rank.");
|
|
}
|
|
if constexpr (!std::same_as<Gravity, nonlinear::MonolithicGravity>) { recipe.gravity.gravity.Validate(); }
|
|
}
|
|
|
|
template <typename State, typename Policy> class PreparedGravityStrategy;
|
|
|
|
template <typename State>
|
|
class PreparedGravityStrategy<State, nonlinear::MonolithicGravity> final {
|
|
public:
|
|
PreparedGravityStrategy(State &, const nonlinear::MonolithicGravity &) { }
|
|
bool Initialize(State &, nonlinear::NumericalStrategyDiagnostics &) { return true; }
|
|
LinearSolveReport Solve(State &state, const LinearSolveControl &control) {
|
|
return state.linearBackend->Solve(state.linearRightHandSide, state.normalizedCorrection, control);
|
|
}
|
|
bool CaptureDirection(State &) { return true; }
|
|
auto PrepareTrial(State &state, double) { return state.PrepareTrial(); }
|
|
bool TrialSolveSucceeded() const noexcept { return true; }
|
|
void CopyDiagnostics(nonlinear::NumericalStrategyDiagnostics &) const { }
|
|
};
|
|
|
|
template <typename State>
|
|
class StellarGravityWorkspace {
|
|
public:
|
|
using Operation = std::remove_cvref_t<decltype(*std::declval<State>().normalizedOperator)>;
|
|
using Inverse = std::remove_cvref_t<decltype(*std::declval<State>().normalizedInverse)>;
|
|
using Problem = std::remove_cvref_t<decltype(std::declval<State>().Problem())>;
|
|
using Subsystem = PreparedGravitySubsystem<Operation, Inverse>;
|
|
|
|
StellarGravityWorkspace(State &state, const nonlinear::GravitySolveOptions &options)
|
|
: m_options(options),
|
|
subsystem(*state.normalizedOperator, *state.normalizedInverse,
|
|
MakeGravityBlockPartition<typename Problem::FormType>(state.Problem().GetManifest().layout()),
|
|
state.Problem().GetCommunicator(), linear::FGMRES(options.backend)),
|
|
gravityResidual(subsystem.GetPartition().GravityResidualSize()),
|
|
gravityCorrection(subsystem.GetPartition().GravityValueSize()),
|
|
expectedResidual(subsystem.GetPartition().GravityResidualSize()),
|
|
gravityAction(subsystem.GetPartition().GravityResidualSize()),
|
|
desiredResidual(subsystem.GetPartition().GravityResidualSize()),
|
|
sourceResidual(subsystem.GetPartition().GravityResidualSize()),
|
|
scaleGravityValues(subsystem.GetPartition().GravityValueSize()),
|
|
scaleGravityAction(subsystem.GetPartition().GravityResidualSize()),
|
|
scaleRootValues(state.Problem().StateSize()),
|
|
fullCorrection(state.Problem().StateSize()), fullAction(state.Problem().EquationSize()) {
|
|
m_options.Validate();
|
|
}
|
|
|
|
bool TrialSolveSucceeded() const noexcept { return trialSolveSucceeded; }
|
|
|
|
void CopyDiagnostics(nonlinear::NumericalStrategyDiagnostics &report) const {
|
|
report.gravitySolves = subsystem.GetSolveCount();
|
|
report.gravityIterations = subsystem.GetIterationCount();
|
|
report.gravityFailures = subsystem.GetFailureCount();
|
|
report.gravityPreparedFailures = preparedFailures;
|
|
report.compactGravityOperator = subsystem.UsesCompactOperator();
|
|
report.compactGravityInverse = subsystem.UsesCompactPreconditioner();
|
|
report.lastLinearReportSpace = lastLinearReportSpace;
|
|
report.gravitySolveSeconds = subsystem.GetSolveSeconds();
|
|
if (const auto &last = subsystem.GetLastSolveReport(); last.has_value()) {
|
|
report.lastGravityTrueResidual = last->trueResidualNorm;
|
|
}
|
|
report.lastGravityCorrectionNorm = lastCorrectionNorm;
|
|
report.lastGravityPreparedDefect = lastPreparedDefect;
|
|
report.lastGravityAffineDisagreement = lastAffineDisagreement;
|
|
report.lastGravityArithmeticAllowance = lastArithmeticAllowance;
|
|
report.lastGravityPreparedTargetThreshold = lastPreparedTargetThreshold;
|
|
report.lastLiftedRelativeTrueResidual = liftedRelativeDefect;
|
|
}
|
|
|
|
protected:
|
|
void Refresh(State &state) {
|
|
// Preparing a trial invalidates both normalized and physical
|
|
// inverses. Refreshing only the restricted wrapper would hide this.
|
|
state.physicalInverse->Refresh();
|
|
subsystem.Refresh();
|
|
}
|
|
|
|
bool SolveCompletion(State &state) {
|
|
subsystem.GetPartition().GatherGravityResiduals(state.trialNormalizedResidual, gravityResidual);
|
|
gravityResidual *= -1.0;
|
|
desiredResidual = 0.0;
|
|
gravityCorrection = 0.0;
|
|
lastSolve = subsystem.Solve(gravityResidual, gravityCorrection, m_options.control);
|
|
return lastSolve.Converged();
|
|
}
|
|
|
|
auto ApplyCorrection(State &state) {
|
|
lastPreparedDefect = lastAffineDisagreement = lastArithmeticAllowance =
|
|
lastPreparedTargetThreshold = std::numeric_limits<double>::quiet_NaN();
|
|
const auto &partition = subsystem.GetPartition();
|
|
partition.InjectGravityValues(gravityCorrection, fullCorrection);
|
|
subsystem.GetOperator().Mult(gravityCorrection, gravityAction);
|
|
partition.GatherGravityResiduals(state.trialNormalizedResidual, expectedResidual);
|
|
sourceResidual = expectedResidual;
|
|
// Split the two gravity fields before taking norms: their large
|
|
// cancellation is precisely what the residual alone cannot expose.
|
|
using Form = typename Problem::FormType;
|
|
using namespace utils::blocks;
|
|
constexpr value_block<type_index_v<gravity::gradient::value, typename Form::value_blocks>> gradient{};
|
|
constexpr value_block<type_index_v<gravity::poisson::value, typename Form::value_blocks>> potential{};
|
|
const auto &layout = state.Problem().GetManifest().layout();
|
|
std::array<double, 8> contributionNorms{};
|
|
const auto measure = [&](const mfem::Vector &values, auto block, const bool source) {
|
|
scaleRootValues = 0.0;
|
|
for (int index = layout.offset(block); index < layout.offset(block) + layout.size(block); ++index)
|
|
scaleRootValues(index) = values(index);
|
|
partition.GatherGravityValues(scaleRootValues, scaleGravityValues);
|
|
subsystem.GetOperator().Mult(scaleGravityValues, scaleGravityAction);
|
|
if (source) sourceResidual -= scaleGravityAction;
|
|
return StrategyNorm(scaleGravityAction, subsystem.GetCommunicator());
|
|
};
|
|
contributionNorms[0] = measure(state.trialNormalizedState, gradient, true);
|
|
contributionNorms[1] = measure(state.trialNormalizedState, potential, true);
|
|
contributionNorms[2] = StrategyNorm(sourceResidual, subsystem.GetCommunicator());
|
|
contributionNorms[3] = measure(fullCorrection, gradient, false);
|
|
contributionNorms[4] = measure(fullCorrection, potential, false);
|
|
contributionNorms[5] = StrategyNorm(expectedResidual, subsystem.GetCommunicator());
|
|
contributionNorms[6] = StrategyNorm(desiredResidual, subsystem.GetCommunicator());
|
|
expectedResidual += gravityAction;
|
|
lastCorrectionNorm = StrategyNorm(gravityCorrection, subsystem.GetCommunicator());
|
|
partition.AddGravityValues(gravityCorrection, state.trialNormalizedState);
|
|
state.normalizedOperator->DenormalizeState(state.trialNormalizedState, state.candidatePhysicalState);
|
|
auto preparation = state.PrepareTrial();
|
|
if (preparation) {
|
|
partition.GatherGravityResiduals(state.trialNormalizedResidual, gravityResidual);
|
|
contributionNorms[7] = StrategyNorm(gravityResidual, subsystem.GetCommunicator());
|
|
expectedResidual -= gravityResidual;
|
|
const double disagreement = StrategyNorm(expectedResidual, subsystem.GetCommunicator());
|
|
gravityResidual -= desiredResidual;
|
|
const double targetDefect = StrategyNorm(gravityResidual, subsystem.GetCommunicator());
|
|
// G is affine in gravity for fixed retained variables. Verify
|
|
// that property against a fresh preparation, including its
|
|
// boundary rows, before trusting a completed candidate.
|
|
const auto assessment = AssessGravityPreparedResidual(
|
|
contributionNorms, disagreement, targetDefect, lastSolve.trueResidualNorm);
|
|
trialSolveSucceeded = assessment.admissible;
|
|
lastPreparedDefect = targetDefect;
|
|
lastAffineDisagreement = disagreement;
|
|
lastArithmeticAllowance = assessment.arithmeticAllowance;
|
|
lastPreparedTargetThreshold = assessment.targetThreshold;
|
|
}
|
|
if (!preparation || !trialSolveSucceeded) { ++preparedFailures; }
|
|
return preparation;
|
|
}
|
|
|
|
nonlinear::GravitySolveOptions m_options;
|
|
Subsystem subsystem;
|
|
mfem::Vector gravityResidual, gravityCorrection, expectedResidual, gravityAction, desiredResidual;
|
|
mfem::Vector sourceResidual, scaleGravityValues, scaleGravityAction, scaleRootValues;
|
|
mfem::Vector fullCorrection, fullAction;
|
|
LinearSolveReport lastSolve{};
|
|
bool trialSolveSucceeded{true};
|
|
double lastCorrectionNorm{0.0};
|
|
double lastPreparedDefect{0.0};
|
|
double lastAffineDisagreement{0.0};
|
|
double lastArithmeticAllowance{0.0};
|
|
double lastPreparedTargetThreshold{0.0};
|
|
double liftedRelativeDefect{0.0};
|
|
std::uint64_t preparedFailures{0};
|
|
nonlinear::LinearReportSpace lastLinearReportSpace{nonlinear::LinearReportSpace::full_system};
|
|
};
|
|
|
|
template <typename State>
|
|
class PreparedGravityStrategy<State, nonlinear::GravityCurvatureCorrection> final
|
|
: public StellarGravityWorkspace<State> {
|
|
using Base = StellarGravityWorkspace<State>;
|
|
public:
|
|
PreparedGravityStrategy(State &state, const nonlinear::GravityCurvatureCorrection &policy)
|
|
: Base(state, policy.gravity), m_curvature(this->subsystem) { }
|
|
bool Initialize(State &, nonlinear::NumericalStrategyDiagnostics &) { return true; }
|
|
LinearSolveReport Solve(State &state, const LinearSolveControl &control) {
|
|
return state.linearBackend->Solve(state.linearRightHandSide, state.normalizedCorrection, control);
|
|
}
|
|
bool CaptureDirection(State &state) {
|
|
this->subsystem.Refresh();
|
|
m_curvature.CaptureBase(state.acceptedNormalizedResidual, state.normalizedCorrection);
|
|
return true;
|
|
}
|
|
auto PrepareTrial(State &state, const double alpha) {
|
|
this->trialSolveSucceeded = true;
|
|
auto preparation = state.PrepareTrial();
|
|
if (!preparation) { return preparation; }
|
|
this->Refresh(state);
|
|
m_curvature.BuildLinearPrediction(alpha, this->desiredResidual);
|
|
this->lastSolve = m_curvature.SolveCorrection(alpha, state.trialNormalizedResidual,
|
|
this->gravityCorrection, this->m_options.control);
|
|
this->trialSolveSucceeded = this->lastSolve.Converged();
|
|
if (!this->trialSolveSucceeded) { return preparation; }
|
|
return this->ApplyCorrection(state);
|
|
}
|
|
private:
|
|
PreparedGravityCurvatureCorrection<typename Base::Subsystem> m_curvature;
|
|
};
|
|
|
|
template <typename State>
|
|
class PreparedGravityStrategy<State, nonlinear::EliminatedGravity> final
|
|
: public StellarGravityWorkspace<State> {
|
|
using Base = StellarGravityWorkspace<State>;
|
|
using ReducedOperator = GravityReducedOperator<typename Base::Subsystem>;
|
|
using ReducedInverse = GravityReducedPreconditioner<typename Base::Subsystem, typename Base::Inverse>;
|
|
using Configuration = std::remove_cvref_t<decltype(std::declval<State>().linearBackend->GetConfiguration())>;
|
|
using Backend = PreparedLinearBackendType<Configuration, ReducedOperator, ReducedInverse>;
|
|
public:
|
|
PreparedGravityStrategy(State &state, const nonlinear::EliminatedGravity &policy)
|
|
: Base(state, policy.gravity), m_operator(this->subsystem, policy.gravity.control),
|
|
m_inverse(m_operator, *state.normalizedInverse),
|
|
m_backend(prepareLinearBackend(state.linearBackend->GetConfiguration(), m_operator, m_inverse,
|
|
state.Problem().GetCommunicator())),
|
|
m_rhs(m_operator.Height()), m_correction(m_operator.Width()) { }
|
|
|
|
bool Initialize(State &state, nonlinear::NumericalStrategyDiagnostics &diagnostics) {
|
|
const auto &partition = this->subsystem.GetPartition();
|
|
partition.GatherGravityResiduals(state.acceptedNormalizedResidual, this->gravityResidual);
|
|
diagnostics.initialGravityResidual = StrategyNorm(this->gravityResidual, this->subsystem.GetCommunicator());
|
|
state.trialNormalizedState = state.acceptedNormalizedState;
|
|
state.trialNormalizedResidual = state.acceptedNormalizedResidual;
|
|
this->trialSolveSucceeded = this->SolveCompletion(state);
|
|
if (!this->trialSolveSucceeded) { return false; }
|
|
const auto preparation = this->ApplyCorrection(state);
|
|
if (!preparation || !this->trialSolveSucceeded ||
|
|
!std::isfinite(StrategyNorm(state.trialNormalizedResidual, this->subsystem.GetCommunicator()))) { return false; }
|
|
partition.GatherGravityResiduals(state.trialNormalizedResidual, this->gravityResidual);
|
|
diagnostics.completedInitialGravityResidual =
|
|
StrategyNorm(this->gravityResidual, this->subsystem.GetCommunicator());
|
|
state.PrepareTrialCommit();
|
|
state.CommitPreparedTrial(preparation->minimumJacobianDeterminant);
|
|
this->subsystem.Refresh();
|
|
return true;
|
|
}
|
|
|
|
LinearSolveReport Solve(State &state, const LinearSolveControl &control) {
|
|
this->subsystem.Refresh();
|
|
try {
|
|
this->lastLinearReportSpace = nonlinear::LinearReportSpace::gravity_subsystem;
|
|
auto inner = m_operator.TryCondenseRightHandSide(state.linearRightHandSide, m_rhs);
|
|
if (!inner.Converged()) { return inner; }
|
|
this->subsystem.GetPartition().GatherRetainedValues(state.normalizedCorrection, m_correction);
|
|
this->lastLinearReportSpace = nonlinear::LinearReportSpace::reduced_system;
|
|
auto report = m_backend.Solve(m_rhs, m_correction, control);
|
|
if (!report.Converged()) { return report; }
|
|
inner = m_operator.TryLiftCorrection(m_correction, state.linearRightHandSide, state.normalizedCorrection);
|
|
if (!inner.Converged()) {
|
|
this->lastLinearReportSpace = nonlinear::LinearReportSpace::gravity_subsystem;
|
|
return inner;
|
|
}
|
|
state.normalizedOperator->Mult(state.normalizedCorrection, this->fullAction);
|
|
this->fullAction -= state.linearRightHandSide;
|
|
const auto communicator = this->subsystem.GetCommunicator();
|
|
const double rhsNorm = StrategyNorm(state.linearRightHandSide, communicator);
|
|
const double defect = StrategyNorm(this->fullAction, communicator);
|
|
this->liftedRelativeDefect = rhsNorm > 0.0 ? defect / rhsNorm : (defect == 0.0 ? 0.0 : INFINITY);
|
|
// Inner Krylov truncation makes Schur actions approximate.
|
|
// The original equations decide whether the lifted result is
|
|
// actually an admissible inexact Newton correction.
|
|
// Keep backend recurrence/initial/true norms in their reduced
|
|
// coordinates. The separate lifted diagnostic reports the
|
|
// original-system check; these are different linear systems.
|
|
if (!std::isfinite(defect) || defect > control.ConvergenceThreshold(rhsNorm)) {
|
|
report.status = std::isfinite(defect) ? LinearSolveStatus::backend_failure : LinearSolveStatus::non_finite;
|
|
}
|
|
return report;
|
|
} catch (const GravityReductionFailure &failure) {
|
|
this->lastLinearReportSpace = nonlinear::LinearReportSpace::gravity_subsystem;
|
|
return failure.GetReport();
|
|
}
|
|
}
|
|
|
|
bool CaptureDirection(State &state) {
|
|
// A positivity repair changes retained coordinates. Lift that
|
|
// repaired tangent again so gravity follows the reduced manifold.
|
|
this->subsystem.GetPartition().GatherRetainedValues(state.normalizedCorrection, m_correction);
|
|
this->lastSolve = m_operator.TryLiftCorrection(m_correction, state.linearRightHandSide,
|
|
state.normalizedCorrection);
|
|
return this->lastSolve.Converged();
|
|
}
|
|
|
|
auto PrepareTrial(State &state, double) {
|
|
this->trialSolveSucceeded = true;
|
|
auto preparation = state.PrepareTrial();
|
|
if (!preparation) { return preparation; }
|
|
this->Refresh(state);
|
|
this->trialSolveSucceeded = this->SolveCompletion(state);
|
|
if (!this->trialSolveSucceeded) { return preparation; }
|
|
return this->ApplyCorrection(state);
|
|
}
|
|
private:
|
|
ReducedOperator m_operator;
|
|
ReducedInverse m_inverse;
|
|
Backend m_backend;
|
|
mfem::Vector m_rhs, m_correction;
|
|
};
|
|
|
|
template <typename State, typename Policy> class PreparedEnthalpyStrategy;
|
|
|
|
template <typename State>
|
|
class PreparedEnthalpyStrategy<State, nonlinear::SampledEnthalpy> final {
|
|
public:
|
|
PreparedEnthalpyStrategy(State &, const nonlinear::SampledEnthalpy &) { }
|
|
bool Check(const mfem::Vector &, nonlinear::NumericalStrategyDiagnostics &) { return true; }
|
|
bool Repair(State &, nonlinear::NumericalStrategyDiagnostics &) { return true; }
|
|
};
|
|
|
|
template <typename State>
|
|
class PreparedEnthalpyStrategy<State, nonlinear::PolynomialEnthalpy> final {
|
|
public:
|
|
PreparedEnthalpyStrategy(State &state, const nonlinear::PolynomialEnthalpy &policy)
|
|
: m_prepared(state, policy.options) { }
|
|
bool Check(const mfem::Vector &state, nonlinear::NumericalStrategyDiagnostics &diagnostics) {
|
|
++diagnostics.polynomialChecks;
|
|
return m_prepared.CheckState(state).status == PolynomialSignStatus::certified_nonnegative;
|
|
}
|
|
bool Repair(State &state, nonlinear::NumericalStrategyDiagnostics &diagnostics) {
|
|
const auto report = m_prepared.RepairDirection(state);
|
|
diagnostics.lastRepairFraction = report.chosenFraction;
|
|
diagnostics.lastMeritDirectionalDerivative = report.meritDirectionalDerivative;
|
|
if (report.chosenFraction > 0.0) { ++diagnostics.directionRepairs; }
|
|
return report.status == EnthalpyDirectionRepairStatus::feasible &&
|
|
std::isfinite(report.meritDirectionalDerivative) && report.meritDirectionalDerivative < 0.0;
|
|
}
|
|
private:
|
|
PreparedStellarEnthalpyStrategy<State> m_prepared;
|
|
};
|
|
|
|
template <typename State, nonlinear::NewtonNumericalConfiguration Recipe>
|
|
class PreparedNewtonStrategy final : public NewtonStrategyStorage {
|
|
public:
|
|
using Enthalpy = typename Recipe::EnthalpyType;
|
|
using Gravity = typename Recipe::GravityType;
|
|
static constexpr bool altersDirection = std::same_as<Enthalpy, nonlinear::PolynomialEnthalpy>;
|
|
|
|
PreparedNewtonStrategy(State &state, const Recipe &recipe)
|
|
: m_enthalpy(state, recipe.enthalpy), m_gravity(state, recipe.gravity) { }
|
|
|
|
bool CheckState(const mfem::Vector &state, nonlinear::NumericalStrategyDiagnostics &diagnostics) {
|
|
return m_enthalpy.Check(state, diagnostics);
|
|
}
|
|
bool Initialize(State &state, nonlinear::NumericalStrategyDiagnostics &diagnostics) {
|
|
const bool ready = m_gravity.Initialize(state, diagnostics);
|
|
m_gravity.CopyDiagnostics(diagnostics);
|
|
return ready;
|
|
}
|
|
LinearSolveReport Solve(State &state, const LinearSolveControl &control,
|
|
nonlinear::NumericalStrategyDiagnostics &diagnostics) {
|
|
auto report = m_gravity.Solve(state, control);
|
|
m_gravity.CopyDiagnostics(diagnostics);
|
|
return report;
|
|
}
|
|
bool CaptureDirection(State &state, nonlinear::NumericalStrategyDiagnostics &diagnostics) {
|
|
m_directionSolveFailed = false;
|
|
if (!m_enthalpy.Repair(state, diagnostics)) { return false; }
|
|
if constexpr (!std::same_as<Gravity, nonlinear::EliminatedGravity> || altersDirection) {
|
|
if (!m_gravity.CaptureDirection(state)) {
|
|
m_directionSolveFailed = true;
|
|
m_gravity.CopyDiagnostics(diagnostics);
|
|
return false;
|
|
}
|
|
}
|
|
if constexpr (altersDirection && std::same_as<Gravity, nonlinear::EliminatedGravity>) {
|
|
mfem::Vector action(state.Problem().EquationSize());
|
|
state.normalizedOperator->Mult(state.normalizedCorrection, action);
|
|
diagnostics.lastMeritDirectionalDerivative =
|
|
StrategyDot(state.acceptedNormalizedResidual, action, state.Problem().GetCommunicator());
|
|
if (!std::isfinite(diagnostics.lastMeritDirectionalDerivative) ||
|
|
diagnostics.lastMeritDirectionalDerivative >= 0.0) { return false; }
|
|
}
|
|
m_gravity.CopyDiagnostics(diagnostics);
|
|
return true;
|
|
}
|
|
auto PrepareTrial(State &state, const double alpha, nonlinear::NumericalStrategyDiagnostics &diagnostics) {
|
|
auto preparation = m_gravity.PrepareTrial(state, alpha);
|
|
m_gravity.CopyDiagnostics(diagnostics);
|
|
return preparation;
|
|
}
|
|
bool TrialSolveSucceeded() const noexcept { return m_gravity.TrialSolveSucceeded(); }
|
|
bool DirectionSolveFailed() const noexcept { return m_directionSolveFailed; }
|
|
private:
|
|
bool m_directionSolveFailed{false};
|
|
PreparedEnthalpyStrategy<State, Enthalpy> m_enthalpy;
|
|
PreparedGravityStrategy<State, Gravity> m_gravity;
|
|
};
|
|
}
|