Emily Boudreaux
cb20faeb3d
The test suite is extensive. Note that LLMs were used to aid in the development of the test suite. All LLM generated code was manually reviewed.
48 lines
2.4 KiB
C++
48 lines
2.4 KiB
C++
#include <print>
|
|
|
|
#include "serif/runtime.hpp"
|
|
#include "serif/model/finite_elements.hpp"
|
|
#include "serif/solver/equilibrium.hpp"
|
|
#include "serif/solver/seed/state.hpp"
|
|
#include "polytropic_model.hpp"
|
|
|
|
/** @brief Solve one nonrotating polytrope with the default Newton, normalization, and inverse choices. */
|
|
int main() {
|
|
serif::initialize();
|
|
namespace solver = serif::solver;
|
|
namespace normalization = solver::normalization;
|
|
|
|
// Physical prescriptions stay fixed while Newton changes the numerical state.
|
|
const auto model = serif::examples::make_polytropic_model();
|
|
const auto radial = solver::seed::make_polytropic_seed(model);
|
|
if (!radial) {
|
|
std::println(stderr, "Cannot prepare the radial star: {}", radial.error());
|
|
return 1;
|
|
}
|
|
const auto finite_elements = serif::model::make_stellar_discretization(model,
|
|
{.radius = radial->scales.radius},
|
|
serif::geometry::deformation::DomainDeformation<serif::geometry::deformation::SphericalHarmonicSurface>{
|
|
.surface = {.maximum_degree = 2}, .stellar_interior = {}, .vacuum = {}});
|
|
auto initial = solver::seed::make_polytropic_state(model, finite_elements, *radial);
|
|
if (!initial) {
|
|
std::println(stderr, "Cannot project the initial star: {}", initial.error());
|
|
return 1;
|
|
}
|
|
|
|
// The reference is numerical bookkeeping; it does not constrain central density.
|
|
auto configuration = serif::model::make_discretization(serif::model::DiscretizationOptions{
|
|
.model = model,
|
|
.finite_elements = finite_elements,
|
|
.normalization = normalization::make_normalization(normalization::NormalizationOptions{
|
|
.scheme = normalization::LaneEmdenApproximation{.n = 1.0, .central_density = radial->scales.central_density}})
|
|
});
|
|
const auto result = solver::solve_equilibrium(std::move(configuration), std::move(*initial));
|
|
std::println("{}; accepted Newton steps: {}; final normalized residual: {}", result.message(),
|
|
result.diagnostics().accepted_iterations, result.diagnostics().final_residual_norm.value_or(-1.0));
|
|
if (!result.converged() && result.diagnostics().last_linear_solve) {
|
|
const auto& linear = *result.diagnostics().last_linear_solve;
|
|
std::println(stderr, "Last linear solve: {} iterations, true relative residual {}", linear.iterations, linear.relative_true_residual_norm);
|
|
}
|
|
return result.converged() ? 0 : 1;
|
|
}
|