From 61fe4d8ae80ed5a0f813065e5f28ec66e598ae6a Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Fri, 27 Mar 2026 11:36:40 -0400 Subject: [PATCH] feat(mapping): major work on parallel version and performance --- default.toml | 13 +- mapping.cpp | 1923 +++++++++++++++++++++++++++++++++++++++----------- meson.build | 14 + 3 files changed, 1519 insertions(+), 431 deletions(-) create mode 100644 meson.build diff --git a/default.toml b/default.toml index 467b617..f79dd93 100644 --- a/default.toml +++ b/default.toml @@ -1,10 +1,15 @@ [main] +core_id = 1 core_steepness = 1.0 +envelope_id = 2 flattening = 0.0 -include_external_domain = false -order = 2 -r_core = 0.1 -r_infinity = 6.0 +include_external_domain = true +inf_bdr_id = 2 +order = 3 +r_core = 0.2 +r_infinity = 3.0 r_instability = 1e-14 r_star = 1.0 refinement_levels = 2 +surface_bdr_id = 1 +vacuum_id = 3 \ No newline at end of file diff --git a/mapping.cpp b/mapping.cpp index 84da56d..8747672 100644 --- a/mapping.cpp +++ b/mapping.cpp @@ -13,15 +13,68 @@ #include #include +#include +#include //endregion -//region Test Utilities +static size_t s_total_num_tests = 0; + +#define LOG_POINT std::println("Log Point: {}", __COUNTER__) + +#define CONCAT_INNER(a, b) a##b +#define CONCAT(a, b) CONCAT_INNER(a, b) + constexpr std::string_view ANSI_GREEN = "\033[32m"; constexpr std::string_view ANSI_RED = "\033[31m"; constexpr std::string_view ANSI_YELLOW = "\033[33m"; constexpr std::string_view ANSI_BLUE = "\033[34m"; +constexpr std::string_view ANSI_MAGENTA = "\033[35m"; +constexpr std::string_view ANSI_CYAN = "\033[36m"; constexpr std::string_view ANSI_RESET = "\033[0m"; +#define MAKE_UNIQUE_VAR_NAME(prefix) CONCAT(prefix, __COUNTER__) + +#define START_TIMER(timer_name) \ +auto start_##timer_name = std::chrono::steady_clock::now() + +#define REPORT_TIMER(timer_name) \ +std::println("Timer '{}': {}", #timer_name, elapsed_##timer_name) + +#define END_TIMER(timer_name) \ +auto end_##timer_name = std::chrono::steady_clock::now(); \ +std::chrono::duration elapsed_##timer_name = end_##timer_name - start_##timer_name; \ +REPORT_TIMER(timer_name) + + +#define RUN_TEST_IMPL(test_name, test_func_w_args, id) \ +const int CONCAT(mpi_world_rank_, id) = mfem::Mpi::WorldRank(); \ +auto CONCAT(start_var_, id) = std::chrono::steady_clock::now(); \ +if (CONCAT(mpi_world_rank_, id) == 0) { \ +std::println("{}===== TEST: {} ====={}", ANSI_MAGENTA, test_name, ANSI_RESET); \ +} \ +test_func_w_args; \ +if (CONCAT(mpi_world_rank_, id) == 0) { \ +auto CONCAT(end_var_, id) = std::chrono::steady_clock::now(); \ +std::chrono::duration CONCAT(elapsed_var_, id) = \ +CONCAT(end_var_, id) - CONCAT(start_var_, id); \ +std::println("{}===== END TEST: {} ({}runtime: {:+0.2f}ms{}) ====={}", ANSI_MAGENTA, \ +test_name, \ +ANSI_YELLOW, \ +CONCAT(elapsed_var_, id).count(), \ +ANSI_MAGENTA, \ +ANSI_RESET); \ +} + +#define RANK_GUARD(proc) \ + if (mfem::Mpi::WorldRank() == 0) { \ + proc \ + } + +#define RUN_TEST(test_name, test_func_w_args) \ +RUN_TEST_IMPL(test_name, test_func_w_args, __COUNTER__) + +//region Test Utilities + enum class TEST_RESULT_TYPE : uint8_t { SUCCESS, FAILURE, @@ -86,7 +139,8 @@ using EOS_P = std::function; * User Args *********************/ struct potential { - double tol; + double rtol; + double atol; int max_iters; }; @@ -104,6 +158,8 @@ struct Args { double mass{}; double c{}; + int quad_boost{0}; + int max_iters{}; double tol{}; }; @@ -129,11 +185,11 @@ enum BoundsError : uint8_t { //region Domain Enums enum class Domains : uint8_t { - NONE = 0, CORE = 1 << 0, ENVELOPE = 1 << 1, VACUUM = 1 << 2, - ALL = 0x7 + STELLAR = CORE | ENVELOPE, + ALL = CORE | ENVELOPE | VACUUM }; inline Domains operator|(Domains lhs, Domains rhs) { @@ -144,7 +200,15 @@ inline Domains operator&(Domains lhs, Domains rhs) { return static_cast(static_cast(lhs) & static_cast(rhs)); } -const Domains STELLAR = Domains::CORE | Domains::ENVELOPE; +enum class Boundaries : uint8_t { + STELLAR_SURFACE = 1, + INF_SURFACE = 2 +}; + +inline int operator-(Boundaries b, const int a) { + return static_cast(static_cast(b) - static_cast(a)); +} + //endregion //region Domain Mapper @@ -159,7 +223,9 @@ public: ) : m_d(nullptr), m_r_star_ref(r_star_ref), - m_r_inf_ref(r_inf_ref) {} + m_r_inf_ref(r_inf_ref) { + InitAllScratchSpaces(); + } explicit DomainMapper( const mfem::GridFunction &d, @@ -169,7 +235,10 @@ public: m_d(&d), m_dim(d.FESpace()->GetMesh()->Dimension()), m_r_star_ref(r_star_ref), - m_r_inf_ref(r_inf_ref) {}; + m_r_inf_ref(r_inf_ref) { + InitAllScratchSpaces(); + }; + [[nodiscard]] bool is_vacuum(const mfem::ElementTransformation &T) const { if (T.ElementType == mfem::ElementTransformation::ELEMENT) { @@ -186,6 +255,7 @@ public: throw std::invalid_argument(err_msg); } m_d = &d; + InvalidateCache(); } [[nodiscard]] bool IsIdentity() const { @@ -194,90 +264,124 @@ public: void ResetDisplacement() { m_d = nullptr; + InvalidateCache(); } void ComputeJacobian(mfem::ElementTransformation &T, mfem::DenseMatrix &J) const { - mfem::DenseMatrix J_D(m_dim, m_dim); J.SetSize(m_dim, m_dim); J = 0.0; - J_D = 0.0; + m_J_D = 0.0; if (IsIdentity()) { for (int i = 0; i < m_dim; ++i) { - J_D(i, i) = 1.0; // Identity mapping + m_J_D(i, i) = 1.0; // Identity mapping } } else { - m_d->GetVectorGradient(T, J_D); + UpdateElementCache(T); + m_dshape.SetSize(m_fe->GetDof(), m_dim); + m_fe->CalcPhysDShape(T, m_dshape); + mfem::MultAtB(m_dof_mat, m_dshape, m_J_D); + for (int i = 0; i < m_dim; ++i) { - J_D(i, i) += 1.0; // Add identity to get the total Jacobian of the mapping + m_J_D(i, i) += 1.0; } + } if (is_vacuum(T)) { - mfem::Vector x_ref(m_dim); - T.Transform(T.GetIntPoint(), x_ref); + T.Transform(T.GetIntPoint(), m_x_ref); - mfem::Vector d_val(m_dim), x_disp(m_dim); if (IsIdentity()) { - x_disp = x_ref; + m_x_disp = m_x_ref; } else { - m_d->GetVectorValue(T, T.GetIntPoint(), d_val); - add(x_ref, d_val, x_disp); + m_shape.SetSize(m_fe->GetDof()); + m_fe->CalcShape(T.GetIntPoint(), m_shape); + m_dof_mat.MultTranspose(m_shape, m_d_val); + add(m_x_ref, m_d_val, m_x_disp); } - mfem::DenseMatrix J_K(m_dim, m_dim); - ComputeKelvinJacobian(x_ref, x_disp, J_K); - mfem::Mult(J_K, J_D, J); + ComputeKelvinJacobian(m_x_ref, m_x_disp, m_J_D, J); } else { - J = J_D; + J = m_J_D; } } double ComputeDetJ(mfem::ElementTransformation& T, const mfem::IntegrationPoint& ip) const { if (IsIdentity() && !is_vacuum(T)) return 1.0; // If no mapping, the determinant of the Jacobian is 1 + T.SetIntPoint(&ip); mfem::DenseMatrix J; ComputeJacobian(T, J); return J.Det(); } void ComputeMappedDiffusionTensor(mfem::ElementTransformation &T, mfem::DenseMatrix &D) const { - mfem::DenseMatrix J(m_dim, m_dim), JInv(m_dim, m_dim); - - ComputeJacobian(T, J); - const double detJ = J.Det(); - mfem::CalcInverse(J, JInv); + ComputeJacobian(T, m_J_temp); + const double detJ = m_J_temp.Det(); + mfem::CalcInverse(m_J_temp, m_JInv_temp); D.SetSize(m_dim, m_dim); - mfem::MultABt(JInv, JInv, D); + mfem::MultABt(m_JInv_temp, m_JInv_temp, D); D *= fabs(detJ); } void ComputeInverseJacobian(mfem::ElementTransformation &T, mfem::DenseMatrix &JInv) const { - mfem::DenseMatrix J(m_dim, m_dim); - ComputeJacobian(T, J); + ComputeJacobian(T, m_J_temp); JInv.SetSize(m_dim, m_dim); - mfem::CalcInverse(J, JInv); + mfem::CalcInverse(m_J_temp, JInv); } void GetPhysicalPoint(mfem::ElementTransformation& T, const mfem::IntegrationPoint& ip, mfem::Vector& x_phys) const { x_phys.SetSize(m_dim); - mfem::Vector x_ref(m_dim); - T.Transform(ip, x_ref); + T.Transform(ip, m_x_ref); if (IsIdentity()) { - x_phys = x_ref; + x_phys = m_x_ref; } else { - mfem::Vector d_val(m_dim); - m_d->GetVectorValue(T, ip, d_val); - add(x_ref, d_val, x_phys); + UpdateElementCache(T); + + m_shape.SetSize(m_fe->GetDof()); + m_fe->CalcShape(ip, m_shape); + + m_dof_mat.MultTranspose(m_shape, m_d_val); + add(m_x_ref, m_d_val, x_phys); } if (is_vacuum(T)) { - ApplyKelvinMapping(x_ref, x_phys); + ApplyKelvinMapping(m_x_ref, x_phys); } } [[nodiscard]] const mfem::GridFunction* GetDisplacement() const { return m_d; } + [[nodiscard]] double GetPhysInfRadius() const { + return 1.0 / (1.0 - m_xi_clamp); + } + + [[nodiscard]] size_t GetCacheHits() const { + return m_cache_hits; + } + + [[nodiscard]] size_t GetCacheMisses() const { + return m_cache_misses; + } + + [[nodiscard]] double GetCacheHitRate() const { + return (static_cast(m_cache_hits)) / static_cast(m_cache_misses + m_cache_hits); + } + + void ResetCacheStats() const { + m_cache_hits = 0; + m_cache_misses = 0; + } + private: + void InitAllScratchSpaces() const { + m_J_D.SetSize(m_dim, m_dim); + m_J_temp.SetSize(m_dim, m_dim); + m_JInv_temp.SetSize(m_dim, m_dim); + m_x_ref.SetSize(m_dim); + m_x_disp.SetSize(m_dim); + m_d_val.SetSize(m_dim); + } + void ApplyKelvinMapping(const mfem::Vector& x_ref, mfem::Vector& x_phys) const { const double r_ref = x_ref.Norml2(); double xi = (r_ref - m_r_star_ref) / (m_r_inf_ref - m_r_star_ref); @@ -286,7 +390,7 @@ private: x_phys *= factor; } - void ComputeKelvinJacobian(const mfem::Vector& x_ref, const mfem::Vector &x_disp, mfem::DenseMatrix& J_K) const { + void ComputeKelvinJacobian(const mfem::Vector& x_ref, const mfem::Vector &x_disp, const mfem::DenseMatrix &J_D, mfem::DenseMatrix& J) const { const double r_ref = x_ref.Norml2(); const double delta_R = m_r_inf_ref - m_r_star_ref; @@ -299,18 +403,51 @@ private: const double dk_dr = m_r_star_ref * (( 1.0 / (delta_R* r_ref * denom * denom)) - ( 1.0 / (r_ref * r_ref * denom))); - J_K.SetSize(m_dim, m_dim); + J.SetSize(m_dim, m_dim); const double outer_factor = dk_dr / r_ref; for (int i = 0; i < m_dim; ++i) { for (int j = 0; j < m_dim; ++j) { - J_K(i, j) = outer_factor * x_disp(i) * x_ref(j); - if (i == j) { - J_K(i, j) += k; - } + J(i, j) = outer_factor * x_disp(i) * x_ref(j) + k * J_D(i, j); + } } } + + void InvalidateCache() const { + m_cached_elem_id = -1; + } + + void UpdateElementCache(const mfem::ElementTransformation& T) const { + if (IsIdentity()) return; + + if (T.ElementNo != m_cached_elem_id || T.ElementType != m_cached_elem_type) { + m_cache_misses++; + m_cached_elem_id = T.ElementNo; + m_cached_elem_type = T.ElementType; + + const mfem::FiniteElementSpace *fes = m_d->FESpace(); + mfem::Array vdofs; + + if (T.ElementType == mfem::ElementTransformation::ELEMENT) { + m_fe = fes->GetFE(m_cached_elem_id); + fes->GetElementVDofs(m_cached_elem_id, vdofs); + } else { + m_fe = fes->GetBE(m_cached_elem_id); + fes->GetBdrElementVDofs(m_cached_elem_id, vdofs); + } + + m_d->GetSubVector(vdofs, m_elem_dofs); + + const int nd = m_fe->GetDof(); + const int vd = fes->GetVDim(); + + m_dof_mat.UseExternalData(m_elem_dofs.GetData(), nd, vd); + } else { + m_cache_hits++; + } + } + private: const mfem::GridFunction *m_d; std::unique_ptr m_internal_d; @@ -318,33 +455,111 @@ private: const int m_vacuum_attr{3}; const double m_r_star_ref{1.0}; const double m_r_inf_ref{2.0}; - const double m_xi_clamp{0.999}; + const double m_xi_clamp{0.99999999}; + + mutable int m_cached_elem_id{-1}; + mutable int m_cached_elem_type{mfem::ElementTransformation::ELEMENT}; + mutable const mfem::FiniteElement* m_fe{nullptr}; + + mutable mfem::Vector m_elem_dofs; + mutable mfem::DenseMatrix m_dof_mat; + mutable mfem::DenseMatrix m_dshape; + mutable mfem::Vector m_shape; + + mutable size_t m_cache_hits{0}; + mutable size_t m_cache_misses{0}; + + mutable mfem::DenseMatrix m_J_D; + mutable mfem::DenseMatrix m_J_temp; + mutable mfem::DenseMatrix m_JInv_temp; + mutable mfem::Vector m_x_ref; + mutable mfem::Vector m_x_disp; + mutable mfem::Vector m_d_val; }; //endregion +/******************** + * Cache Types + *********************/ + //region State Types +class MappedScalarCoefficient; + /******************** * State Types *********************/ +struct LORPrecWrapper : public mfem::Solver { + mfem::Solver& m_amg; + explicit LORPrecWrapper(mfem::Solver& amg) : mfem::Solver(amg.Height(), amg.Width()), m_amg(amg) {} + + void SetOperator(const Operator &op) override {}; + + void Mult(const mfem::Vector &x, mfem::Vector &y) const override { m_amg.Mult(x, y); } +}; + +struct GravityContext { + std::unique_ptr ho_laplacian; + std::unique_ptr lor_laplacian; + + std::unique_ptr phi; + + + + mfem::Array ho_ess_tdof_list; + mfem::Array lor_ess_tdof_list; + + mfem::Array stellar_mask; + + std::unique_ptr diff_coeff; + + std::unique_ptr b; + std::unique_ptr rho_coeff; + std::unique_ptr four_pi_G_coeff; + std::unique_ptr rhs_coeff; + std::unique_ptr mapped_rhs_coeff; + std::unique_ptr unit_coeff; + + std::unique_ptr amg_prec; + std::unique_ptr amg_wrapper; + std::unique_ptr solver; + + mfem::OperatorPtr A_ho; + mfem::OperatorPtr A_lor; + + mfem::Vector B_true; + mfem::Vector X_true; + +}; + +struct BoundaryContext { + mfem::Array inf_bounds; + mfem::Array stellar_bounds; +}; + struct FEM { - std::unique_ptr mesh; + std::unique_ptr mesh; std::unique_ptr H1_fec; - std::unique_ptr H1_fes; - std::unique_ptr Vec_H1_fes; + std::unique_ptr H1_fes; + std::unique_ptr Vec_H1_fes; + std::unique_ptr H1_lor_disc; + const mfem::ParFiniteElementSpace* H1_lor_fes{nullptr}; std::unique_ptr mapping; - mfem::Array block_offsets; - std::unique_ptr reference_x; + mfem::Array block_true_offsets; + std::unique_ptr reference_x; mfem::Vector com; mfem::DenseMatrix Q; mfem::Array ess_tdof_x; - std::unique_ptr> star_marker; - std::unique_ptr> vacuum_marker; - std::unique_ptr> surface_marker; - std::unique_ptr> inf_marker; + + int int_order{3}; + std::unique_ptr int_rule; + + GravityContext gravity_context; + BoundaryContext boundary_context; + [[nodiscard]] bool okay() const { return (mesh != nullptr) && (H1_fec != nullptr) && (H1_fes != nullptr) && (Vec_H1_fes != nullptr); } @@ -357,7 +572,7 @@ struct CoupledState { mfem::GridFunction d; // Stability depends on solving for the displacement vector not the nodal position vector, those live on a reference domain. explicit CoupledState(const FEM& fem) { - U = std::make_unique(fem.block_offsets); + U = std::make_unique(fem.block_true_offsets); rho.MakeRef(fem.H1_fes.get(), U->GetBlock(0), 0); d.MakeRef(fem.Vec_H1_fes.get(), U->GetBlock(1), 0); *U = 0.0; @@ -365,13 +580,14 @@ struct CoupledState { } }; + //endregion //region Function Definitions /******************** * Core Setup Functions *********************/ -FEM setup_fem(const std::string& filename, bool verbose=true); +FEM setup_fem(const std::string& filename, const Args &args); /******************** * Utility Functions @@ -383,21 +599,22 @@ void get_physical_coordinates(const mfem::GridFunction& reference_pos, const mfe void populate_element_mask(const FEM& fem, Domains domain, mfem::Array& mask); std::expected DiscoverBounds(const mfem::Mesh *mesh, int vacuum_attr); +int get_mesh_order(const mfem::Mesh &mesh); +void conserve_mass(const FEM& fem, mfem::GridFunction& rho, double target_mass); + /******************** * Physics Functions *********************/ double centrifugal_potential(const mfem::Vector& phys_x, double omega); double get_moment_of_inertia(const FEM& fem, const mfem::GridFunction& rho); double oblate_spheroid_surface_potential(const mfem::Vector& x, double a, double c, double total_mass); -std::unique_ptr grav_potential(const FEM& fem, const Args &args, const mfem::GridFunction& rho); -std::unique_ptr get_potential(const FEM& fem, const Args& args, const mfem::GridFunction& rho); + +const mfem::GridFunction &grav_potential(FEM &fem, const Args &args, const mfem::GridFunction &rho, + bool phi_warm = false); + +mfem::GridFunction get_potential(FEM &fem, const Args &args, const mfem::GridFunction &rho); mfem::DenseMatrix compute_quadrupole_moment_tensor(const FEM& fem, const mfem::GridFunction& rho, const mfem::Vector& com); double l2_multipole_potential(const FEM& fem, double total_mass, const mfem::Vector& phys_x); - -/******************** - * Conservation Functions - *********************/ -void conserve_mass(const FEM& fem, mfem::GridFunction& rho, double target_mass); //endregion //region Mapping Coefficients @@ -560,9 +777,11 @@ class FluidIntegrator : public mfem::NonlinearFormIntegrator { using Scalar = EOS_T::value_type; public: explicit FluidIntegrator( + const FEM& fem, EOS_P eos, const DomainMapper* map = nullptr ) : + m_fem(fem), m_eos(std::move(eos)), m_map(map) {}; @@ -577,7 +796,7 @@ public: elvect.SetSize(dof); elvect = 0.0; - const mfem::IntegrationRule *ir = IntRule ? IntRule : &mfem::IntRules.Get(el.GetGeomType(), 2 * el.GetOrder() + 1); + const mfem::IntegrationRule *ir = m_fem.int_rule.get(); mfem::Vector shape(dof); for (int i = 0; i < ir->GetNPoints(); i++) { @@ -606,7 +825,7 @@ public: elmat.SetSize(dof); elmat = 0.0; - const mfem::IntegrationRule *ir = IntRule ? IntRule : &mfem::IntRules.Get(el.GetGeomType(), 2 * el.GetOrder() + 1); + const mfem::IntegrationRule *ir = m_fem.int_rule.get(); mfem::Vector shape(dof); for (int i = 0; i < ir->GetNPoints(); i++) { @@ -643,6 +862,7 @@ public: void clear_mapping() { m_map = nullptr; } private: + const FEM& m_fem; EOS_P m_eos; const DomainMapper* m_map{nullptr}; }; @@ -961,7 +1181,7 @@ public: const EOS_P& eos_pressure, const double alpha ) : - Operator(fem.block_offsets.Last()), + Operator(fem.block_true_offsets.Last()), m_fem(fem), m_args(args), m_eos_enthalpy(eos_enthalpy), @@ -970,21 +1190,34 @@ public: m_fluid_nlf(m_fem.H1_fes.get()), m_reference_stiffness(m_fem.Vec_H1_fes.get()) { - m_fluid_nlf.AddDomainIntegrator(new FluidIntegrator(m_eos_enthalpy, m_fem.mapping.get())); + auto* fluid_integrator = new FluidIntegrator(m_fem, m_eos_enthalpy, m_fem.has_mapping() ? m_fem.mapping.get() : nullptr); + fluid_integrator->SetIntRule(fem.int_rule.get()); - m_reference_stiffness.AddDomainIntegrator(new mfem::VectorMassIntegrator(*m_alpha)); - m_reference_stiffness.AddDomainIntegrator(new mfem::VectorDiffusionIntegrator()); + populate_element_mask(m_fem, Domains::STELLAR, m_stellar_mask); + m_bdr_mask.SetSize(m_fem.mesh->bdr_attributes.Max()); + m_bdr_mask = 0; + m_bdr_mask[0] = 1; + + m_fluid_nlf.AddDomainIntegrator(fluid_integrator, m_stellar_mask); + + auto* alpha_integrator = new mfem::VectorMassIntegrator(*m_alpha); + alpha_integrator->SetIntRule(fem.int_rule.get()); + + auto* diff_integrator = new mfem::VectorDiffusionIntegrator(); + diff_integrator->SetIntRule(fem.int_rule.get()); + + m_reference_stiffness.AddDomainIntegrator(alpha_integrator, m_stellar_mask); + m_reference_stiffness.AddDomainIntegrator(diff_integrator, m_stellar_mask); m_reference_stiffness.Assemble(); m_reference_stiffness.Finalize(); - std::println("ResidualOperator initialized with XAD-enabled FluidIntegrator and reference stiffness."); }; void Mult(const mfem::Vector &u, mfem::Vector &r) const override { mfem::GridFunction rho, d; - rho.MakeRef(m_fem.H1_fes.get(), u.GetData() + m_fem.block_offsets[0]); - d.MakeRef(m_fem.Vec_H1_fes.get(), u.GetData() + m_fem.block_offsets[1]); - double lambda = u(m_fem.block_offsets[2]); + rho.MakeRef(m_fem.H1_fes.get(), u.GetData() + m_fem.block_true_offsets[0]); + d.MakeRef(m_fem.Vec_H1_fes.get(), u.GetData() + m_fem.block_true_offsets[1]); + double lambda = u(m_fem.block_true_offsets[2]); m_fem.mapping->SetDisplacement(d); @@ -992,29 +1225,36 @@ public: m_fem.Q = compute_quadrupole_moment_tensor(m_fem, rho, m_fem.com); mfem::GridFunction r_rho, r_d; - r_rho.MakeRef(m_fem.H1_fes.get(), r.GetData() + m_fem.block_offsets[0]); - r_d.MakeRef(m_fem.Vec_H1_fes.get(), r.GetData() + m_fem.block_offsets[1]); + r_rho.MakeRef(m_fem.H1_fes.get(), r.GetData() + m_fem.block_true_offsets[0]); + r_d.MakeRef(m_fem.Vec_H1_fes.get(), r.GetData() + m_fem.block_true_offsets[1]); - double &r_lambda = r(m_fem.block_offsets[2]); + double &r_lambda = r(m_fem.block_true_offsets[2]); m_fluid_nlf.Mult(rho, r_rho); auto phi = get_potential(m_fem, m_args, rho); - mfem::GridFunctionCoefficient phi_c(phi.get()); + mfem::GridFunctionCoefficient phi_c(&phi); MappedScalarCoefficient mapped_phi_c(*m_fem.mapping, phi_c); mfem::LinearForm phi_lf(m_fem.H1_fes.get()); - phi_lf.AddDomainIntegrator(new mfem::DomainLFIntegrator(mapped_phi_c)); + + auto* potential_lf_integrator = new mfem::DomainLFIntegrator(mapped_phi_c); + potential_lf_integrator -> SetIntRule(m_fem.int_rule.get()); + + phi_lf.AddDomainIntegrator(potential_lf_integrator, m_stellar_mask); phi_lf.Assemble(); r_rho += phi_lf; mfem::ConstantCoefficient lambda_c(lambda); MappedScalarCoefficient mapped_lambda_c(*m_fem.mapping, lambda_c); mfem::LinearForm mass_grad_lf(m_fem.H1_fes.get()); - mass_grad_lf.AddDomainIntegrator(new mfem::DomainLFIntegrator(mapped_lambda_c)); + auto* lambda_lf_integrator = new mfem::DomainLFIntegrator(mapped_lambda_c); + lambda_lf_integrator->SetIntRule(m_fem.int_rule.get()); + + mass_grad_lf.AddDomainIntegrator(lambda_lf_integrator, m_stellar_mask); mass_grad_lf.Assemble(); // ReSharper disable once CppDFAUnusedValue - r_rho += mass_grad_lf; + r_rho -= mass_grad_lf; m_reference_stiffness.Mult(d, r_d); @@ -1026,7 +1266,11 @@ public: m_args.c ); mfem::LinearForm pbf_lf(m_fem.Vec_H1_fes.get()); - pbf_lf.AddBoundaryIntegrator(new mfem::VectorBoundaryLFIntegrator(pbf)); + + auto* vec_bdr_integrator = new mfem::VectorBoundaryLFIntegrator(pbf); + vec_bdr_integrator->SetIntRule(m_fem.int_rule.get()); + + pbf_lf.AddBoundaryIntegrator(vec_bdr_integrator, m_bdr_mask); pbf_lf.Assemble(); r_d -= pbf_lf; @@ -1035,21 +1279,26 @@ public: r_d(m_fem.ess_tdof_x[i]) = 0.0; } - double current_mass = domain_integrate_grid_function(m_fem, rho); + double current_mass = domain_integrate_grid_function(m_fem, rho, Domains::STELLAR); r_lambda = current_mass - m_args.mass; + + auto com = get_com(m_fem, rho); + std::cout << std::flush; } [[nodiscard]] Operator& GetGradient(const mfem::Vector &u) const override { + mfem::Array stellar_mask; + populate_element_mask(m_fem, Domains::STELLAR, stellar_mask); + mfem::GridFunction rho, d; - rho.MakeRef(m_fem.H1_fes.get(), u.GetData() + m_fem.block_offsets[0]); - d.MakeRef(m_fem.Vec_H1_fes.get(), u.GetData() + m_fem.block_offsets[1]); + rho.MakeRef(m_fem.H1_fes.get(), u.GetData() + m_fem.block_true_offsets[0]); + d.MakeRef(m_fem.Vec_H1_fes.get(), u.GetData() + m_fem.block_true_offsets[1]); m_fem.mapping->SetDisplacement(d); - m_approx_jacobian = std::make_unique(m_fem.block_offsets); + m_approx_jacobian = std::make_unique(m_fem.block_true_offsets); - const mfem::GridFunction grad(m_fem.Vec_H1_fes.get(), u.GetData() + m_fem.block_offsets[1]); - // m_fem.mesh->SetNodes(grad); + const mfem::GridFunction grad(m_fem.Vec_H1_fes.get(), u.GetData() + m_fem.block_true_offsets[1]); m_approx_jacobian->SetBlock(0, 0, &m_fluid_nlf.GetGradient(rho)); m_approx_jacobian->SetBlock(1, 1, &m_reference_stiffness); @@ -1060,13 +1309,18 @@ public: mfem::ConstantCoefficient one(1.0); MappedScalarCoefficient mapped_b_vec(*m_fem.mapping, one); mfem::LinearForm b_lf(m_fem.H1_fes.get()); - b_lf.AddDomainIntegrator(new mfem::DomainLFIntegrator(mapped_b_vec)); + + auto* lf_integrator = new mfem::DomainLFIntegrator(mapped_b_vec); + lf_integrator->SetIntRule(m_fem.int_rule.get()); + b_lf.AddDomainIntegrator(lf_integrator, stellar_mask); b_lf.Assemble(); B_vec += b_lf; - m_B_vec_op_col.SetVector(B_vec, true); m_B_vec_op_row.SetVector(B_vec, false); + B_vec *= -1.0; + m_B_vec_op_col.SetVector(B_vec, true); + m_approx_jacobian->SetBlock(0, 2, &m_B_vec_op_col); m_approx_jacobian->SetBlock(2, 0, &m_B_vec_op_row); @@ -1086,6 +1340,9 @@ private: const EOS_P& m_eos_pressure; const std::unique_ptr m_alpha; + mutable mfem::Array m_stellar_mask; + mutable mfem::Array m_bdr_mask; + mutable mfem::NonlinearForm m_fluid_nlf; mutable mfem::BilinearForm m_reference_stiffness; mutable std::unique_ptr m_approx_jacobian = nullptr; @@ -1095,18 +1352,54 @@ private: mutable std::unique_ptr> m_C; mutable std::unique_ptr> m_D; }; + +class JFNKOperator : public mfem::Operator { +public: + JFNKOperator(const Operator& F, const mfem::Vector& u0, const mfem::Vector& F_u0) : + Operator(F.Width()), m_F(F), m_u0(u0), m_F_u0(F_u0), m_w(u0.Size()), m_eps(1e-8) {} + + void Mult(const mfem::Vector &dx, mfem::Vector &Jdx) const override { + add(m_u0, m_eps, dx, m_w); + m_F.Mult(m_w, Jdx); + + Jdx.Add(-1.0, m_F_u0); + Jdx /= m_eps; + } +private: + const Operator& m_F; + const mfem::Vector& m_u0; + const mfem::Vector& m_F_u0; + mutable mfem::Vector m_w; + double m_eps; +}; + //endregion //region Utility Functions -FEM setup_fem(const std::string& filename, const bool verbose) { +FEM setup_fem(const std::string& filename, const Args &args) { FEM fem; - fem.mesh = std::make_unique(filename, 0, 0); - fem.mesh->EnsureNodes(); - int dim = fem.mesh->Dimension(); - fem.H1_fec = std::make_unique(2, dim); - fem.H1_fes = std::make_unique(fem.mesh.get(), fem.H1_fec.get()); - fem.Vec_H1_fes = std::make_unique(fem.mesh.get(), fem.H1_fec.get(), dim, mfem::Ordering::byNODES); + //================================================================== + // Section 1: Mesh and FE Space Setup + //================================================================== + mfem::Mesh serial_mesh(filename, 0, 0); + fem.mesh = std::make_unique(MPI_COMM_WORLD, serial_mesh); + fem.mesh->EnsureNodes(); + const int geom_order = get_mesh_order(*fem.mesh); + int dim = fem.mesh->Dimension(); + + // Assume solution space order == geometrical order + fem.H1_fec = std::make_unique(geom_order, dim); + fem.H1_fes = std::make_unique(fem.mesh.get(), fem.H1_fec.get()); + fem.Vec_H1_fes = std::make_unique(fem.mesh.get(), fem.H1_fec.get(), dim, mfem::Ordering::byNODES); + + // LOR discretization for fast preconditioning + fem.H1_lor_disc = std::make_unique(*fem.H1_fes); + fem.H1_lor_fes = &fem.H1_lor_disc->GetParFESpace(); + + //================================================================== + // Section 2: Domain Mapping + //================================================================== auto [r_star_ref, r_inf_ref] = DiscoverBounds(fem.mesh.get(), 3) .or_else([](const BoundsError& err)->std::expected { throw std::runtime_error("Unable to determine vacuum domain reference boundary..."); @@ -1114,77 +1407,250 @@ FEM setup_fem(const std::string& filename, const bool verbose) { fem.mapping = std::make_unique(r_star_ref, r_inf_ref); - fem.block_offsets.SetSize(4); - fem.block_offsets[0] = 0; - fem.block_offsets[1] = fem.H1_fes->GetTrueVSize(); - fem.block_offsets[2] = fem.block_offsets[1] + fem.Vec_H1_fes->GetTrueVSize(); - fem.block_offsets[3] = fem.block_offsets[2] + 1; + //================================================================== + // Section 3: Multi-physics Block-offsets + // Layout: [rho (scalar H1) | d (vector H1) | λ (scalar)] + //================================================================== + fem.block_true_offsets.SetSize(4); + fem.block_true_offsets[0] = 0; + fem.block_true_offsets[1] = fem.H1_fes->GetTrueVSize(); + fem.block_true_offsets[2] = fem.block_true_offsets[1] + fem.Vec_H1_fes->GetTrueVSize(); + fem.block_true_offsets[3] = fem.block_true_offsets[2] + 1; - // Setup Domain Markers to constrain integration - fem.star_marker = std::make_unique>(fem.mesh->attributes.Max()); - fem.star_marker->operator=(0); - fem.star_marker.operator->()[0] = 1; // core - fem.star_marker.operator->()[1] = 1; // envelope - - fem.vacuum_marker = std::make_unique>(fem.mesh->attributes.Max()); - fem.vacuum_marker->operator=(0); - fem.vacuum_marker.operator->()[2] = 1; // vacuum - - fem.surface_marker = std::make_unique>(fem.mesh->bdr_attributes.Max()); - fem.surface_marker->operator=(0); - fem.surface_marker.operator->()[0] = 1; // surface - - fem.vacuum_marker = std::make_unique>(fem.mesh->bdr_attributes.Max()); - fem.vacuum_marker->operator=(0); - fem.vacuum_marker.operator->()[1] = 1; // Infinity - - // Initial COM should be 0, initial mass distribution should be uniform + //================================================================== + // Section 4: Multipole BC setup. + // Assumes: + // - COM starts at origin + // - Geometry starts out as a spherically symmetric + //================================================================== fem.com.SetSize(dim); fem.com = 0.0; fem.Q.SetSize(dim, dim); fem.Q = 0.0; - fem.reference_x = std::make_unique(fem.Vec_H1_fes.get()); + fem.reference_x = std::make_unique(fem.Vec_H1_fes.get()); fem.mesh->GetNodes(*fem.reference_x); + //================================================================== + // Section 5: Integration Rules + // Assumes: + // - Mesh is composed of hex elements (CUBE) + //================================================================== + MFEM_ASSERT(fem.mesh->GetElementGeometry(0) == mfem::Geometry::CUBE, "Currently only hexahedral meshes are supported"); + const int element_order = fem.H1_fes->GetMaxElementOrder(); + fem.int_order = 2 * element_order + geom_order - 2 + args.quad_boost; + + fem.int_rule = std::make_unique(mfem::IntRules.Get(mfem::Geometry::CUBE, fem.int_order)); + + //================================================================== + // Section 6: Essential Degrees of Freedom and Boundary Markers + //================================================================== fem.ess_tdof_x.SetSize(0); // No essential boundary conditions for the displacement, the null space here is handled with a weak penalty term + + populate_element_mask(fem, Domains::STELLAR, fem.gravity_context.stellar_mask); + + const int n_bdr_attrs = fem.mesh->bdr_attributes.Max(); + fem.boundary_context.inf_bounds.SetSize(n_bdr_attrs); + fem.boundary_context.stellar_bounds.SetSize(n_bdr_attrs); + + fem.boundary_context.inf_bounds = 0; + fem.boundary_context.stellar_bounds = 0; + + fem.boundary_context.inf_bounds[Boundaries::INF_SURFACE - 1] = 1; // Vacuum boundary + fem.boundary_context.stellar_bounds[Boundaries::STELLAR_SURFACE - 1] = 1; // Stellar boundary + + fem.H1_fes->GetEssentialTrueDofs(fem.boundary_context.inf_bounds, fem.gravity_context.ho_ess_tdof_list); + fem.H1_lor_fes->GetEssentialTrueDofs(fem.boundary_context.inf_bounds, fem.gravity_context.lor_ess_tdof_list); + + //================================================================== + // Section 7: Gravity Solution Field + //================================================================== + fem.gravity_context.phi = std::make_unique(fem.H1_fes.get()); + *fem.gravity_context.phi = 0.0; + + //================================================================== + // Section 8: Laplacian Coefficients + //================================================================== + fem.gravity_context.unit_coeff = std::make_unique(1.0); + if (fem.has_mapping()) { + fem.gravity_context.diff_coeff = std::make_unique(*fem.mapping, *fem.gravity_context.unit_coeff, fem.mesh->Dimension()); + } else { + throw std::runtime_error("Unable to determine domain reference boundary or mapping..."); + } + + + //================================================================== + // Section 9: High order Laplacian (Partial Assembly) + //================================================================== + fem.gravity_context.ho_laplacian = std::make_unique(fem.H1_fes.get()); + { + auto* ho_di = new mfem::DiffusionIntegrator(*fem.gravity_context.diff_coeff); + ho_di->SetIntRule(fem.int_rule.get()); + fem.gravity_context.ho_laplacian->AddDomainIntegrator(ho_di); + } + fem.gravity_context.ho_laplacian->SetAssemblyLevel(mfem::AssemblyLevel::PARTIAL); + fem.gravity_context.ho_laplacian->Assemble(); + + //================================================================== + // Section 10: Low order Laplacian (preconditioning) + //================================================================== + fem.gravity_context.lor_laplacian = std::make_unique(&fem.H1_lor_disc->GetParFESpace()); + { + auto* lo_di = new mfem::DiffusionIntegrator(*fem.gravity_context.diff_coeff); + fem.gravity_context.lor_laplacian->AddDomainIntegrator(lo_di); + } + fem.gravity_context.lor_laplacian->Assemble(); + fem.gravity_context.lor_laplacian->Finalize(); + + //================================================================== + // Section 11: Constrained System Operators + //================================================================== + fem.gravity_context.ho_laplacian->FormSystemMatrix( + fem.gravity_context.ho_ess_tdof_list, + fem.gravity_context.A_ho + ); + + fem.gravity_context.lor_laplacian->FormSystemMatrix( + fem.gravity_context.lor_ess_tdof_list, + fem.gravity_context.A_lor + ); + + //================================================================== + // Section 12: AMG Preconditioners and CG Solver + //================================================================== + fem.gravity_context.amg_prec = std::make_unique(*fem.gravity_context.A_lor.As()); + fem.gravity_context.amg_prec->SetPrintLevel(0); + + fem.gravity_context.amg_wrapper = std::make_unique(*fem.gravity_context.amg_prec); + + fem.gravity_context.solver = std::make_unique(MPI_COMM_WORLD); + fem.gravity_context.solver->SetOperator(*fem.gravity_context.A_ho.Ptr()); + fem.gravity_context.solver->SetPreconditioner(*fem.gravity_context.amg_wrapper); + fem.gravity_context.solver->SetRelTol(args.p.rtol); + fem.gravity_context.solver->SetAbsTol(args.p.atol); + fem.gravity_context.solver->SetMaxIter(1000); + fem.gravity_context.solver->SetKDim(50); + fem.gravity_context.solver->SetPrintLevel(0); + + //================================================================== + // Section 13: RHS Coefficient Chain + // Assumes: + // - fem.gravity_context.rho_coeff->SetGridFunction(&rho) + // must be called before each assembly of the RHS + //================================================================== + fem.gravity_context.four_pi_G_coeff = std::make_unique(-4.0 * M_PI * G); + fem.gravity_context.rho_coeff = std::make_unique(); + fem.gravity_context.rhs_coeff = std::make_unique(*fem.gravity_context.four_pi_G_coeff, *fem.gravity_context.rho_coeff); + fem.gravity_context.mapped_rhs_coeff = std::make_unique(*fem.mapping, *fem.gravity_context.rhs_coeff, MappedScalarCoefficient::EVAL_POINTS::REFERENCE); + + //================================================================== + // Section 14: RHS Linear Form + //================================================================== + fem.gravity_context.b = std::make_unique(fem.H1_fes.get()); + + { + auto* grav_rhs_integrator = new mfem::DomainLFIntegrator(*fem.gravity_context.mapped_rhs_coeff); + grav_rhs_integrator->SetIntRule(fem.int_rule.get()); + fem.gravity_context.b->AddDomainIntegrator(grav_rhs_integrator, fem.gravity_context.stellar_mask); + } + + + //========================================================= + // Section 15: Diagnostic Output + //========================================================= + std::println( + "{}Setup (rank: {}): element_order={}, geom_order={}, int_order={}, dim={}, " + "r_star={:.3f}, r_inf={:.3f}, HO_ndofs={}, LOR_ndofs={}{}", + ANSI_BLUE, + mfem::Mpi::WorldRank(), + element_order, geom_order, fem.int_order, dim, + r_star_ref, r_inf_ref, + fem.H1_fes->GlobalTrueVSize(), // Global DOF count for HO space + fem.H1_lor_fes->GlobalTrueVSize(), // Global DOF count for LOR space + ANSI_RESET); return fem; } +void update_stiffness_matrix(FEM& fem) { + auto& ctx = fem.gravity_context; + + const mfem::GridFunction* saved_d = fem.mapping->GetDisplacement(); + fem.mapping->ResetDisplacement(); + + + ctx.lor_laplacian = std::make_unique(&fem.H1_lor_disc->GetParFESpace()); + ctx.lor_laplacian->AddDomainIntegrator(new mfem::DiffusionIntegrator(*ctx.diff_coeff)); + ctx.lor_laplacian->Assemble(); + ctx.lor_laplacian->Finalize(); + + ctx.A_lor.Clear(); + ctx.lor_laplacian->FormSystemMatrix(ctx.lor_ess_tdof_list, ctx.A_lor); + + ctx.amg_prec = std::make_unique(*ctx.A_lor.As()); + ctx.amg_prec->SetPrintLevel(0); + + ctx.amg_wrapper = std::make_unique(*ctx.amg_prec); + ctx.solver->SetPreconditioner(*ctx.amg_wrapper); + + if (saved_d) { + fem.mapping->SetDisplacement(*saved_d); + } + + ctx.ho_laplacian = std::make_unique(fem.H1_fes.get()); + auto* ho_di = new mfem::DiffusionIntegrator(*ctx.diff_coeff); + ho_di->SetIntRule(fem.int_rule.get()); + ctx.ho_laplacian->AddDomainIntegrator(ho_di); + ctx.ho_laplacian->SetAssemblyLevel(mfem::AssemblyLevel::PARTIAL); + ctx.ho_laplacian->Assemble(); + + ctx.A_ho.Clear(); + ctx.ho_laplacian->FormSystemMatrix(ctx.ho_ess_tdof_list, ctx.A_ho); + ctx.solver->SetOperator(*ctx.A_ho.Ptr());} + void view_mesh(const std::string& host, int port, const mfem::Mesh& mesh, const mfem::GridFunction& gf, const std::string& title) { mfem::socketstream sol_sock(host.c_str(), port); if (!sol_sock.is_open()) return; sol_sock << "solution\n" << mesh << gf; - sol_sock << "window_title '" << title << "'\n" << std::flush; + sol_sock << "window_title '" << title << "\n" << std::flush; } double domain_integrate_grid_function(const FEM& fem, const mfem::GridFunction& gf, Domains domain) { mfem::LinearForm lf(fem.H1_fes.get()); mfem::GridFunctionCoefficient gf_c(&gf); - double integral; + double local_integral; mfem::Array elem_markers; populate_element_mask(fem, domain, elem_markers); if (fem.has_mapping()) { MappedScalarCoefficient mapped_gf_c(*fem.mapping, gf_c); - lf.AddDomainIntegrator(new mfem::DomainLFIntegrator(mapped_gf_c), elem_markers); + + // ReSharper disable once CppDFAMemoryLeak // Disabled because MFEM takes ownership so memory is not leaked + auto* lf_integrator = new mfem::DomainLFIntegrator(mapped_gf_c); + lf_integrator->SetIntRule(fem.int_rule.get()); + lf.AddDomainIntegrator(lf_integrator, elem_markers); lf.Assemble(); - integral = lf.Sum(); + + + local_integral = lf.Sum(); } else { lf.AddDomainIntegrator(new mfem::DomainLFIntegrator(gf_c), elem_markers); lf.Assemble(); - integral = lf.Sum(); + local_integral = lf.Sum(); } - return integral; + + double global_integral = 0.0; + MPI_Allreduce(&local_integral, &global_integral, 1, MPI_DOUBLE, MPI_SUM, fem.H1_fes->GetComm()); + return global_integral; } mfem::Vector get_com(const FEM& fem, const mfem::GridFunction &rho) { const int dim = fem.mesh->Dimension(); - mfem::Vector com(dim); - com = 0.0; - double total_mass = 0.0; + mfem::Vector local_com(dim); + local_com = 0.0; + double local_mass = 0.0; for (int i = 0; i < fem.H1_fes->GetNE(); ++i) { + if (fem.mesh->GetAttribute(i) == 3) continue; mfem::ElementTransformation *trans = fem.H1_fes->GetElementTransformation(i); - const mfem::IntegrationRule &ir = mfem::IntRules.Get(trans->GetGeometryType(), fem.H1_fes->GetOrder(0) + trans->OrderW()); + const mfem::IntegrationRule &ir = *fem.int_rule; for (int j = 0; j < ir.GetNPoints(); ++j) { const mfem::IntegrationPoint &ip = ir.IntPoint(j); @@ -1204,15 +1670,29 @@ mfem::Vector get_com(const FEM& fem, const mfem::GridFunction &rho) { } const double mass_term = rho_val * weight; - total_mass += mass_term; + local_mass += mass_term; for (int d = 0; d < dim; ++d) { - com(d) += phys_point(d) * mass_term; + local_com(d) += phys_point(d) * mass_term; } } } - com /= total_mass; - return com; + + double global_mass = 0.0; + mfem::Vector global_com(dim); + MPI_Comm comm = fem.H1_fes->GetComm(); + + MPI_Allreduce(&local_mass, &global_mass, 1, MPI_DOUBLE, MPI_SUM, comm); + + MPI_Allreduce(local_com.GetData(), global_com.GetData(), dim, MPI_DOUBLE, MPI_SUM, comm); + + if (global_mass > 1e-18) { + global_com /= global_mass; + } else { + global_com = 0.0; + } + + return global_com; } void get_physical_coordinates(const mfem::GridFunction& reference_pos, const mfem::GridFunction& displacement, mfem::GridFunction& physical_pos) { @@ -1237,9 +1717,8 @@ void populate_element_mask(const FEM &fem, Domains domain, mfem::Array &mas } std::expected DiscoverBounds(const mfem::Mesh *mesh, const int vacuum_attr) { - double min_r = std::numeric_limits::max(); - double max_r = -std::numeric_limits::max(); - + double local_min_r = std::numeric_limits::max(); + double local_max_r = -std::numeric_limits::max(); bool found_vacuum = false; for (int i = 0; i < mesh->GetNE(); ++i) { @@ -1247,27 +1726,48 @@ std::expected DiscoverBounds(const mfem::Mesh *mesh, const found_vacuum = true; mfem::Array vertices; mesh->GetElementVertices(i, vertices); - for (const int v: vertices) { + for (const int v : vertices) { const double* coords = mesh->GetVertex(v); - double r = std::sqrt(coords[0] * coords[0] + coords[1] * coords[1] + coords[2] * coords[2]); - min_r = std::min(min_r, r); - max_r = std::max(max_r, r); + double r = std::sqrt(coords[0]*coords[0] + coords[1]*coords[1] + coords[2]*coords[2]); + local_min_r = std::min(local_min_r, r); + local_max_r = std::max(local_max_r, r); } } } - if (found_vacuum) { - return Bounds(min_r, max_r); + + double global_min_r, global_max_r; + int global_found_vacuum; + int l_found = found_vacuum ? 1 : 0; + + MPI_Comm comm = MPI_COMM_WORLD; + if (const auto* pmesh = dynamic_cast(mesh)) { + comm = pmesh->GetComm(); } - return std::unexpected(BoundsError::CANNOT_FIND_VACUUM); + + MPI_Allreduce(&local_min_r, &global_min_r, 1, MPI_DOUBLE, MPI_MIN, comm); + MPI_Allreduce(&local_max_r, &global_max_r, 1, MPI_DOUBLE, MPI_MAX, comm); + MPI_Allreduce(&l_found, &global_found_vacuum, 1, MPI_INT, MPI_MAX, comm); + + if (global_found_vacuum) { + return Bounds(global_min_r, global_max_r); + } + return std::unexpected(CANNOT_FIND_VACUUM); +} + +int get_mesh_order(const mfem::Mesh &mesh) { + if (mesh.GetNodes() != nullptr) { + return mesh.GetNodes() -> FESpace() -> GetMaxElementOrder(); + } + return 1; } void conserve_mass(const FEM& fem, mfem::GridFunction& rho, const double target_mass) { - if (const double current_mass = domain_integrate_grid_function(fem, rho); current_mass > 1e-15) rho *= (target_mass / current_mass); + if (const double current_mass = domain_integrate_grid_function(fem, rho, Domains::STELLAR); current_mass > 1e-15) rho *= (target_mass / current_mass); } //endregion //region Physics Functions -double centrifugal_potential(const mfem::Vector& phys_x, double omega) { +double centrifugal_potential(const mfem::Vector& phys_x, const double omega) { const double s2 = std::pow(phys_x(0), 2) + std::pow(phys_x(1), 2); return -0.5 * s2 * std::pow(omega, 2); } @@ -1290,6 +1790,7 @@ double get_moment_of_inertia(const FEM& fem, const mfem::GridFunction& rho) { mfem::LinearForm I_lf(fem.H1_fes.get()); double I = 0.0; + // TODO: Need to filter here to just the stellar domain and also update the IntRule if (fem.has_mapping()) { MappedScalarCoefficient mapped_integrand(*fem.mapping, I_integrand); I_lf.AddDomainIntegrator(new mfem::DomainLFIntegrator(mapped_integrand)); @@ -1305,215 +1806,44 @@ double get_moment_of_inertia(const FEM& fem, const mfem::GridFunction& rho) { } //endregion -//region Analytic External Potential for an oblate spheroid (doi.10.16.j.pss.2018.01.005) - -inline double sq(const double x) { return x * x; } - -double eccentricity(const double a, const double c) { - return std::sqrt(1.0 - sq(c)/sq(a)); -} - -double equatorial_radius(const mfem::Vector& x) { - return std::sqrt(sq(x(0)) + sq(x(1))); -} - -double beta_ext(const mfem::Vector& x, const double a, const double c) { - const double e = eccentricity(a, c); - const double z = x(2); - const double r = equatorial_radius(x); - const double ae_sq = sq(a*e); - - if (std::abs(r) < 1e-12) { - return std::atan2(a*e, std::abs(x(2))); +//region potentials +const mfem::GridFunction &grav_potential(FEM &fem, const Args &args, const mfem::GridFunction &rho, const bool phi_warm) { + auto& ctx = fem.gravity_context; + if (!phi_warm) { + *ctx.phi = 0.0; } - const double fraction = (sq(z) + ae_sq) / sq(r); - const double first_term = -fraction; + double total_mass = domain_integrate_grid_function(fem, rho, Domains::STELLAR); - const double radical_one = sq(1.0 + fraction); - const double radical_two = 4.0 * ae_sq / sq(r); - const double second_term = std::sqrt(std::max(0.0, radical_one - radical_two)); - const double cos_term = first_term + second_term; - return 0.5 * std::acos(std::clamp(cos_term, -1.0, 1.0)); -} - -double oblate_spheroid_surface_potential(const mfem::Vector& x, const double a, const double c, const double total_mass) { - const double beta = beta_ext(x, a, c); - const double r = equatorial_radius(x); - const double e = eccentricity(a, c); - const double z = x(2); - - const double G_M = G * total_mass; - const double e3a3 = std::pow(e, 3) * std::pow(a, 3); - - const double t1 = (3.0 * G_M * beta) / (2.0 * e * a); - - const double t2aa = (3.0 * G_M * sq(r)) / (2.0 * e3a3); - const double t2ab = beta - std::sin(beta) * std::cos(beta); - - const double t2ba = (3.0 * G_M * sq(z)) / e3a3; - const double t2bb = std::tan(beta) - beta; - - return -t1 + 0.5 * (t2aa * t2ab + t2ba * t2bb); -} -//endregion End Analytic Potential - -//region potentials -std::unique_ptr grav_potential(const FEM& fem, const Args& args, const mfem::GridFunction& rho) { - auto phi = std::make_unique(fem.H1_fes.get()); - - mfem::Array ess_bdr(fem.mesh->bdr_attributes.Max()); - ess_bdr = 0; - assert(ess_bdr.Size() == 2); - ess_bdr[1] = 1; // The only boundary condition is that the potential goes to 0 at infinity - - // mfem::ConstantCoefficient zero(0.0); - // phi->ProjectBdrCoefficient(zero, ess_bdr); // Set the potential to 0 at infinity; - - double total_mass = domain_integrate_grid_function(fem, rho, STELLAR); auto bdr_func = [&fem, total_mass](const mfem::Vector& x) { - double r = x.Norml2(); - double theta = std::atan2(x(1), x(0)); - double phi = std::acos(x(2) / r); - double val = l2_multipole_potential(fem, total_mass, x); - // std::println("{},{},{},{},{}", r, theta, phi, val, -G*total_mass / r); return l2_multipole_potential(fem, total_mass, x); }; - std::unique_ptr phi_bdr_coeff; - if (fem.has_mapping()) { - phi_bdr_coeff = std::make_unique(*fem.mapping, bdr_func); - } else { - phi_bdr_coeff = std::make_unique(bdr_func); - } - phi->ProjectBdrCoefficient(*phi_bdr_coeff, ess_bdr); + PhysicalPositionFunctionCoefficient phi_bdr_coeff(*fem.mapping, bdr_func); + ctx.phi->ProjectBdrCoefficient(phi_bdr_coeff, fem.boundary_context.inf_bounds); + ctx.rho_coeff->SetGridFunction(&rho); + ctx.b->Assemble(); - mfem::Array ess_tdof_list; - fem.H1_fes->GetEssentialTrueDofs(ess_bdr, ess_tdof_list); + ctx.ho_laplacian->FormLinearSystem( + ctx.ho_ess_tdof_list, + *ctx.phi, + *ctx.b, + ctx.A_ho, + ctx.X_true, + ctx.B_true + ); - mfem::Array stellar_marker; - populate_element_mask(fem, STELLAR, stellar_marker); + ctx.solver->SetOperator(*ctx.A_ho.Ptr()); - auto laplacian = std::make_unique(fem.H1_fes.get()); + ctx.solver->Mult(ctx.B_true, ctx.X_true); + ctx.ho_laplacian->RecoverFEMSolution(ctx.X_true, *ctx.b, *ctx.phi); - mfem::ConstantCoefficient one_coeff(1.0); - std::unique_ptr mapped_diff_coeff; - if (fem.has_mapping()) { - mapped_diff_coeff = std::make_unique(*fem.mapping, one_coeff, fem.mesh->Dimension()); - laplacian->AddDomainIntegrator(new mfem::DiffusionIntegrator(*mapped_diff_coeff)); // The laplacian is global - } else { - laplacian->AddDomainIntegrator(new mfem::DiffusionIntegrator()); - } - - laplacian->Assemble(); - laplacian->Finalize(); - - mfem::LinearForm b(fem.H1_fes.get()); - mfem::ConstantCoefficient four_pi_G(-4.0 * M_PI * G); - mfem::GridFunctionCoefficient rho_coeff(&rho); - mfem::ProductCoefficient rhs_coeff(rho_coeff, four_pi_G); - - auto mapped_rhs = std::make_unique(*fem.mapping, rhs_coeff, MappedScalarCoefficient::EVAL_POINTS::REFERENCE); - if (fem.has_mapping()) { - b.AddDomainIntegrator(new mfem::DomainLFIntegrator(*mapped_rhs), stellar_marker); // The mass contribution is local to the stellar domain - } else { - b.AddDomainIntegrator(new mfem::DomainLFIntegrator(rhs_coeff), stellar_marker); - } - - b.Assemble(); - - mfem::OperatorPtr A; - mfem::Vector B, X; - laplacian->FormLinearSystem(ess_tdof_list, *phi, b, A, X, B); - - mfem::GSSmoother prec; - mfem::CGSolver cg; - cg.SetPreconditioner(prec); - cg.SetOperator(*A); - cg.SetRelTol(args.p.tol); - cg.SetMaxIter(args.p.max_iters); - cg.Mult(B, X); - - laplacian->RecoverFEMSolution(X, b, *phi); - return phi; + return *ctx.phi; } -// std::unique_ptr grav_potential(const FEM& fem, const Args &args, const mfem::GridFunction& rho, std::optional oblate) { -// auto phi = std::make_unique(fem.H1_fes.get()); -// -// mfem::Array ess_bdr(fem.mesh->bdr_attributes.Max()); -// ess_bdr = 1; -// -// mfem::GridFunctionCoefficient rho_coeff(&rho); -// double total_mass = domain_integrate_grid_function(fem, rho); -// -// auto grav_potential = [&fem, &total_mass, &oblate](const mfem::Vector& x) { -// if (!oblate.has_value() or oblate->use == false) { -// return l2_multipole_potential(fem, total_mass, x); -// } -// return oblate_spheroid_surface_potential(x, oblate->a, oblate->c, total_mass); -// }; -// -// std::unique_ptr phi_bdr_coeff; -// if (fem.has_mapping()) { -// phi_bdr_coeff = std::make_unique(*fem.mapping, grav_potential); -// } else { -// phi_bdr_coeff = std::make_unique(grav_potential); -// } -// phi->ProjectBdrCoefficient(*phi_bdr_coeff, ess_bdr); -// -// mfem::Array ess_tdof_list; -// fem.H1_fes->GetEssentialTrueDofs(ess_bdr, ess_tdof_list); -// -// auto laplacian = std::make_unique(fem.H1_fes.get()); -// // ReSharper disable once CppTooWideScope -// std::unique_ptr one_coeff; -// // ReSharper disable once CppTooWideScope -// std::unique_ptr mapped_diff_coeff; -// if (fem.has_mapping()) { -// one_coeff = std::make_unique(1.0); -// mapped_diff_coeff = std::make_unique(*fem.mapping, *one_coeff, fem.mesh->Dimension()); -// laplacian->AddDomainIntegrator(new mfem::DiffusionIntegrator(*mapped_diff_coeff)); -// } else { -// laplacian->AddDomainIntegrator(new mfem::DiffusionIntegrator()); -// } -// laplacian->Assemble(); -// laplacian->Finalize(); -// -// mfem::ConstantCoefficient four_pi_G(-4.0 * M_PI * G); -// mfem::ProductCoefficient rhs_coeff(rho_coeff, four_pi_G); -// mfem::LinearForm b(fem.H1_fes.get()); -// -// // ReSharper disable once CppTooWideScope -// std::unique_ptr mapped_rhs; -// if (fem.has_mapping()) { -// mapped_rhs = std::make_unique(*fem.mapping, rhs_coeff, MappedScalarCoefficient::EVAL_POINTS::REFERENCE); -// b.AddDomainIntegrator(new mfem::DomainLFIntegrator(*mapped_rhs)); -// } else { -// b.AddDomainIntegrator(new mfem::DomainLFIntegrator(rhs_coeff)); -// } -// b.Assemble(); -// -// mfem::OperatorPtr A; -// mfem::Vector B, X; -// laplacian->FormLinearSystem(ess_tdof_list, *phi, b, A, X, B); -// -// mfem::GSSmoother prec; -// mfem::CGSolver cg; -// cg.SetPreconditioner(prec); -// cg.SetOperator(*A); -// cg.SetRelTol(args.p.tol); -// cg.SetMaxIter(args.p.max_iters); -// cg.SetPrintLevel(0); -// -// cg.Mult(B, X); -// laplacian->RecoverFEMSolution(X, b, *phi); -// return phi; -// } - -std::unique_ptr get_potential(const FEM &fem, const Args &args, const mfem::GridFunction &rho) { +mfem::GridFunction get_potential(FEM &fem, const Args &args, const mfem::GridFunction &rho) { auto phi = grav_potential(fem, args, rho); if (args.r.enabled) { @@ -1529,7 +1859,7 @@ std::unique_ptr get_potential(const FEM &fem, const Args &ar } mfem::GridFunction centrifugal_gf(fem.H1_fes.get()); centrifugal_gf.ProjectCoefficient(*centrifugal_coeff); - (*phi) += centrifugal_gf; + phi += centrifugal_gf; } return phi; @@ -1537,12 +1867,15 @@ std::unique_ptr get_potential(const FEM &fem, const Args &ar mfem::DenseMatrix compute_quadrupole_moment_tensor(const FEM& fem, const mfem::GridFunction& rho, const mfem::Vector& com) { const int dim = fem.mesh->Dimension(); - mfem::DenseMatrix Q(dim, dim); - Q = 0.0; + mfem::DenseMatrix local_Q(dim, dim); + local_Q = 0.0; for (int i = 0; i < fem.H1_fes->GetNE(); ++i) { + if (fem.mesh->GetAttribute(i) == 3) continue; + mfem::ElementTransformation *trans = fem.mesh->GetElementTransformation(i); - const mfem::IntegrationRule &ir = mfem::IntRules.Get(trans->GetGeometryType(), 2 * fem.H1_fes->GetOrder(0) + trans->OrderW()); + const mfem::IntegrationRule &ir = *fem.int_rule; + for (int j = 0; j < ir.GetNPoints(); ++j) { const mfem::IntegrationPoint &ip = ir.IntPoint(j); trans->SetIntPoint(&ip); @@ -1574,14 +1907,17 @@ mfem::DenseMatrix compute_quadrupole_moment_tensor(const FEM& fem, const mfem::G for (int n = 0; n < dim; ++n) { const double delta = (m == n) ? 1.0 : 0.0; const double contrib = 3.0 * x_prime(m) * x_prime(n) - delta * r_sq; - Q(m, n) += rho_val * contrib * weight; + local_Q(m, n) += rho_val * contrib * weight; } } } } - return Q; -} + mfem::DenseMatrix global_Q(dim, dim); + MPI_Allreduce(local_Q.GetData(), global_Q.GetData(), dim * dim, MPI_DOUBLE, MPI_SUM, fem.H1_fes->GetComm()); + + return global_Q; +} double l2_multipole_potential(const FEM &fem, const double total_mass, const mfem::Vector &phys_x) { const double r = phys_x.Norml2(); if (r < 1e-12) return 0.0; @@ -1598,7 +1934,7 @@ double l2_multipole_potential(const FEM &fem, const double total_mass, const mfe } } - const double l2_contrib = (G / (2.0 * std::pow(r, 3))) * l2_mult_factor; + const double l2_contrib = - (G / (2.0 * std::pow(r, 3))) * l2_mult_factor; const double l0_contrib = -G * total_mass / r; @@ -1618,10 +1954,10 @@ void test_mesh_load(const FEM& fem) { const int n_vector = fem.Vec_H1_fes->GetTrueVSize(); if (n_vector != dim * n_scalar) ++failed; - if (fem.block_offsets[0] != 0) ++failed; - if (fem.block_offsets[1] != n_scalar) ++failed; - if (fem.block_offsets[2] != n_scalar + n_vector) ++failed; - if (fem.block_offsets[3] != n_scalar + n_vector + 1) ++failed; + if (fem.block_true_offsets[0] != 0) ++failed; + if (fem.block_true_offsets[1] != n_scalar) ++failed; + if (fem.block_true_offsets[2] != n_scalar + n_vector) ++failed; + if (fem.block_true_offsets[3] != n_scalar + n_vector + 1) ++failed; constexpr size_t num_tests = 6; auto result_type = TEST_RESULT_TYPE::FAILURE; @@ -1630,15 +1966,18 @@ void test_mesh_load(const FEM& fem) { } else if (failed < num_tests) { result_type = TEST_RESULT_TYPE::PARTIAL; } - std::string test_msg = fmt_test_msg("Mesh Load Test", result_type, failed, num_tests); - std::println("{}", test_msg); + + RANK_GUARD( + std::string test_msg = fmt_test_msg("Mesh Load Test", result_type, failed, num_tests); + std::println("{}", test_msg); + ) assert(dim == 3); assert(n_vector == (n_scalar * dim)); - assert (fem.block_offsets[0] == 0); - assert (fem.block_offsets[1] == n_scalar); - assert (fem.block_offsets[2] == n_scalar + n_vector); - assert (fem.block_offsets[3] == n_scalar + n_vector + 1); + assert (fem.block_true_offsets[0] == 0); + assert (fem.block_true_offsets[1] == n_scalar); + assert (fem.block_true_offsets[2] == n_scalar + n_vector); + assert (fem.block_true_offsets[3] == n_scalar + n_vector + 1); } void test_ref_coord_storage(const FEM& fem) { @@ -1668,8 +2007,10 @@ void test_ref_coord_storage(const FEM& fem) { result_type = TEST_RESULT_TYPE::PARTIAL; } - std::string test_msg = fmt_test_msg("Mesh Ref Coord", result_type, failed, num_tests); - std::println("{}", test_msg); + RANK_GUARD( + std::string test_msg = fmt_test_msg("Mesh Ref Coord", result_type, failed, num_tests); + std::println("{}", test_msg); + ) } void test_reference_volume_integral(const FEM& fem) { @@ -1678,7 +2019,7 @@ void test_reference_volume_integral(const FEM& fem) { mfem::GridFunction ones(fem.H1_fes.get()); ones = 1.0; - double vol = domain_integrate_grid_function(fem, ones, STELLAR); + double vol = domain_integrate_grid_function(fem, ones, Domains::STELLAR); double expected = (4.0/3.0) * M_PI * std::pow(RADIUS, 3.0); double rel_err = std::abs(vol - expected) / expected; @@ -1690,11 +2031,13 @@ void test_reference_volume_integral(const FEM& fem) { result_type = TEST_RESULT_TYPE::SUCCESS; } - std::println("{}", fmt_test_msg("Reference Volume Integral", result_type, failed, num_tests)); + RANK_GUARD( + std::println("{}", fmt_test_msg("Reference Volume Integral", result_type, failed, num_tests)); - if (result_type == TEST_RESULT_TYPE::FAILURE) { - std::println("\tFAILURE: Volume: {}, Expected: {}, Error (rel): {}", vol, expected, rel_err); - } + if (result_type == TEST_RESULT_TYPE::FAILURE) { + std::println("\tFAILURE: Volume: {}, Expected: {}, Error (rel): {}", vol, expected, rel_err); + } + ) } void test_spherically_symmetric_com(const FEM& fem) { @@ -1718,11 +2061,13 @@ void test_spherically_symmetric_com(const FEM& fem) { result_type = TEST_RESULT_TYPE::PARTIAL; } - std::println("{}", fmt_test_msg("Uniform COM", result_type, failed, num_tests)); + RANK_GUARD( + std::println("{}", fmt_test_msg("Uniform COM", result_type, failed, num_tests)); - if (result_type == TEST_RESULT_TYPE::FAILURE) { - std::println("\t COM=<{:+0.3E}, {:+0.3E}, {:+0.3E}>", com(0), com(1), com(2)); - } + if (result_type == TEST_RESULT_TYPE::FAILURE) { + std::println("\t COM=<{:+0.3E}, {:+0.3E}, {:+0.3E}>", com(0), com(1), com(2)); + } + ) } void test_com_variance_to_displacement(const FEM& fem) { @@ -1750,11 +2095,13 @@ void test_com_variance_to_displacement(const FEM& fem) { result_type = TEST_RESULT_TYPE::PARTIAL; } - std::println("{}", fmt_test_msg("COM variance to displacement", result_type, failed, num_tests)); + RANK_GUARD( + std::println("{}", fmt_test_msg("COM variance to displacement", result_type, failed, num_tests)); - if (result_type == TEST_RESULT_TYPE::FAILURE) { - std::println("\tFAILURE COM=<{:+0.2E}, {:+0.2E}, {:+0.2E}>", mapped_com(0), mapped_com(1), mapped_com(2)); - } + if (result_type == TEST_RESULT_TYPE::FAILURE) { + std::println("\tFAILURE COM=<{:+0.2E}, {:+0.2E}, {:+0.2E}>", mapped_com(0), mapped_com(1), mapped_com(2)); + } + ) fem.mapping->ResetDisplacement(); } @@ -1768,7 +2115,7 @@ void test_volume_invariance_to_displacement(const FEM& fem) { mfem::GridFunction ones(fem.H1_fes.get()); ones = 1.0; - double mapped_vol = domain_integrate_grid_function(fem, ones, STELLAR); + double mapped_vol = domain_integrate_grid_function(fem, ones, Domains::STELLAR); double expected = (4.0/3.0) * M_PI * std::pow(RADIUS, 3.0); double rel_err = std::abs(mapped_vol - expected) / expected; @@ -1779,11 +2126,13 @@ void test_volume_invariance_to_displacement(const FEM& fem) { result_type = TEST_RESULT_TYPE::SUCCESS; } - std::println("{}", fmt_test_msg("Invariance of volume against translation", result_type, failed, num_tests)); + RANK_GUARD( + std::println("{}", fmt_test_msg("Invariance of volume against translation", result_type, failed, num_tests)); - if (result_type == TEST_RESULT_TYPE::FAILURE) { - std::println("\tFAILURE: Volume: {}, Expected: {}", mapped_vol, expected); - } + if (result_type == TEST_RESULT_TYPE::FAILURE) { + std::println("\tFAILURE: Volume: {}, Expected: {}", mapped_vol, expected); + } + ) fem.mapping->ResetDisplacement(); } @@ -1813,12 +2162,14 @@ void test_volume_ellipsoid_deformation(const FEM& fem) { ++num_tests; mfem::GridFunction ones(fem.H1_fes.get()); ones = 1.0; - const double mapped_vol = domain_integrate_grid_function(fem, ones, STELLAR); + const double mapped_vol = domain_integrate_grid_function(fem, ones, Domains::STELLAR); const double rel_err = std::abs(mapped_vol - expected_vol) / expected_vol; if (rel_err > 1e-3) { ++failed; - std::println("\tFAILURE Test 1: Mapped volume = {:.6f}, expected = {:.6f}, rel_err = {:.2e}", - mapped_vol, expected_vol, rel_err); + RANK_GUARD( + std::println("\tFAILURE Test 1: Mapped volume = {:.6f}, expected = {:.6f}, rel_err = {:.2e}", + mapped_vol, expected_vol, rel_err); + ) } } @@ -1840,12 +2191,11 @@ void test_volume_ellipsoid_deformation(const FEM& fem) { ); x_phys_sq.ProjectCoefficient(x_phys_sq_coeff); - const double mapped_x2_integral = domain_integrate_grid_function(fem, x_phys_sq, STELLAR); - const double rel_err = std::abs(mapped_x2_integral - expected_x2_integral) / expected_x2_integral; - if (rel_err > 1e-3) { + const double mapped_x2_integral = domain_integrate_grid_function(fem, x_phys_sq, Domains::STELLAR); + if (const double rel_err = std::abs(mapped_x2_integral - expected_x2_integral) / expected_x2_integral; rel_err > 1e-3) { ++failed; - std::println("\tFAILURE Test 2: integral x_phys^2 = {:.6f}, expected = {:.6f}, rel_err = {:.2e}", - mapped_x2_integral, expected_x2_integral, rel_err); + RANK_GUARD(std::println("\tFAILURE Test 2: integral x_phys^2 = {:.6f}, expected = {:.6f}, rel_err = {:.2e}", + mapped_x2_integral, expected_x2_integral, rel_err);) } } @@ -1854,8 +2204,12 @@ void test_volume_ellipsoid_deformation(const FEM& fem) { constexpr double expected_detJ = a * b * c; double max_detJ_err = 0.0; for (int e = 0; e < std::min(5, fem.mesh->GetNE()); ++e) { + if (fem.mesh->GetAttribute(e) == 3) { // We want to ignore vacuum elements for this test + e--; + continue; + } auto* trans = fem.mesh->GetElementTransformation(e); - const auto& ir = mfem::IntRules.Get(trans->GetGeometryType(), 2); + const auto& ir = *fem.int_rule; for (int q = 0; q < ir.GetNPoints(); ++q) { const auto& ip = ir.IntPoint(q); trans->SetIntPoint(&ip); @@ -1865,7 +2219,7 @@ void test_volume_ellipsoid_deformation(const FEM& fem) { } if (max_detJ_err > 1e-10) { ++failed; - std::println("\tFAILURE Test 3: max pointwise |det(J) - a*b*c| = {:.2e}", max_detJ_err); + RANK_GUARD(std::println("\tFAILURE Test 3: max pointwise |det(J) - a*b*c| = {:.2e}", max_detJ_err);) } } @@ -1876,12 +2230,15 @@ void test_volume_ellipsoid_deformation(const FEM& fem) { result_type = TEST_RESULT_TYPE::PARTIAL; } - std::println("{}", fmt_test_msg("Volume under ellipsoidal deformation", result_type, failed, num_tests)); + RANK_GUARD(std::println("{}", fmt_test_msg("Volume under ellipsoidal deformation", result_type, failed, num_tests));) fem.mapping->ResetDisplacement(); } void test_uniform_potential(FEM& fem, const Args& args) { + fem.mapping->ResetDisplacement(); + update_stiffness_matrix(fem); + const double analytic_vol = (4.0/3.0) * M_PI * std::pow(RADIUS, 3); const double rho0 = MASS / analytic_vol; @@ -1893,53 +2250,69 @@ void test_uniform_potential(FEM& fem, const Args& args) { const auto phi = grav_potential(fem, args, rho_uniform); - double max_abs_err = 0.0; - double max_rel_err = 0.0; + double local_max_abs_err = 0.0; + double local_max_rel_err = 0.0; constexpr double tol = 1e-3; - size_t failed = 0; - size_t num_tests = 0; + size_t local_failed = 0; + size_t local_num_tests = 0; + const size_t num_elemIDs = std::min(30, fem.mesh->GetNE()); for (int elemID = 0; elemID < num_elemIDs; ++elemID) { - num_tests++; + local_num_tests++; auto* trans = fem.mesh->GetElementTransformation(elemID); const auto& ir = mfem::IntRules.Get(trans->GetGeometryType(), 2); const auto& ip = ir.IntPoint(0); trans->SetIntPoint(&ip); - mfem::Vector x_ref, x_phys; - trans->Transform(ip, x_ref); + mfem::Vector x_phys; fem.mapping->GetPhysicalPoint(*trans, ip, x_phys); const double r = x_phys.Norml2(); - if (r < 1e-9) continue; const double phi_analytic = (-G * MASS / (2.0 * std::pow(RADIUS, 3.0))) * (3.0*RADIUS*RADIUS - r*r); - - const double phi_fem = phi->GetValue(elemID, ip); + const double phi_fem = phi.GetValue(elemID, ip); // Evaluates local ParGridFunction part const double abs_err = std::abs(phi_fem - phi_analytic); const double rel_err = abs_err / std::abs(phi_analytic); - max_abs_err = std::max(max_abs_err, abs_err); - max_rel_err = std::max(max_rel_err, rel_err); + local_max_abs_err = std::max(local_max_abs_err, abs_err); + local_max_rel_err = std::max(local_max_rel_err, rel_err); - if (rel_err > tol) ++failed; + if (rel_err > tol) ++local_failed; } + double global_max_abs_err = 0.0; + double global_max_rel_err = 0.0; + long global_failed = 0; + long global_num_tests = 0; + + MPI_Comm comm = fem.H1_fes->GetComm(); + + MPI_Allreduce(&local_max_abs_err, &global_max_abs_err, 1, MPI_DOUBLE, MPI_MAX, comm); + MPI_Allreduce(&local_max_rel_err, &global_max_rel_err, 1, MPI_DOUBLE, MPI_MAX, comm); + + long l_failed = static_cast(local_failed); + long l_tests = static_cast(local_num_tests); + MPI_Allreduce(&l_failed, &global_failed, 1, MPI_LONG, MPI_SUM, comm); + MPI_Allreduce(&l_tests, &global_num_tests, 1, MPI_LONG, MPI_SUM, comm); + auto result_type = TEST_RESULT_TYPE::FAILURE; - if (failed == 0) { + if (global_failed == 0) { result_type = TEST_RESULT_TYPE::SUCCESS; - } else if (failed < num_tests) { + } else if (global_failed < global_num_tests) { result_type = TEST_RESULT_TYPE::PARTIAL; } - std::println("{}", fmt_test_msg("Test Uniform Potential", result_type, failed, num_tests)); + RANK_GUARD( + std::println("{}", fmt_test_msg("Test Uniform Potential", result_type, global_failed, global_num_tests)); - if (result_type == TEST_RESULT_TYPE::FAILURE) { - std::println("\tFAILURE: max abs error: {:+0.2E}, max rel error: {:+0.2E}", max_abs_err, max_rel_err); - } + if (result_type == TEST_RESULT_TYPE::FAILURE || result_type == TEST_RESULT_TYPE::PARTIAL) { + std::println("\tFAILURE: global max abs error: {:+0.2E}, global max rel error: {:+0.2E}", + global_max_abs_err, global_max_rel_err); + } + ) } void test_ellipsoidal_potential(FEM& fem, const Args& args) { @@ -1959,6 +2332,7 @@ void test_ellipsoidal_potential(FEM& fem, const Args& args) { }); ellipsoidal_disp.ProjectCoefficient(disp_coeff); fem.mapping->SetDisplacement(ellipsoidal_disp); + update_stiffness_matrix(fem); mfem::GridFunction rho(fem.H1_fes.get()); rho = rho0; @@ -1994,7 +2368,7 @@ void test_ellipsoidal_potential(FEM& fem, const Args& args) { const double z2 = x_phys(2)*x_phys(2); const double phi_analytic = -M_PI * G * rho0 * (a*a*I_const - A_R * R2 - A_z * z2); - const double phi_fem = phi->GetValue(elemID, ip); + const double phi_fem = phi.GetValue(elemID, ip); const double rel_err = std::abs(phi_fem - phi_analytic) / std::abs(phi_analytic); max_rel_err = std::max(max_rel_err, rel_err); total_err += rel_err; @@ -2009,15 +2383,698 @@ void test_ellipsoidal_potential(FEM& fem, const Args& args) { result_type = TEST_RESULT_TYPE::PARTIAL; } - std::println("{}", fmt_test_msg("Test Ellipsoidal Potential", result_type, failed, num_tests)); - if (result_type == TEST_RESULT_TYPE::FAILURE) { - std::println("\tFAILURE: max rel error: {:+0.2E}, mean rel error: {:+0.2E}", max_rel_err, total_err/static_cast(num_tests)); + RANK_GUARD( + std::println("{}", fmt_test_msg("Test Ellipsoidal Potential", result_type, failed, num_tests)); + if (result_type == TEST_RESULT_TYPE::FAILURE) { + std::println("\tFAILURE: max rel error: {:+0.2E}, mean rel error: {:+0.2E}", max_rel_err, total_err/static_cast(num_tests)); + } + ) +} + +void test_ferrers_sphere_potential(FEM& fem, const Args& args) { + constexpr double R = RADIUS; + constexpr double rho0 = 1.0; + + [[maybe_unused]] const double expected_mass = (8.0 / 15.0) * M_PI * rho0 * std::pow(R, 3.0); + + fem.mapping->ResetDisplacement(); + update_stiffness_matrix(fem); + + mfem::FunctionCoefficient rho_coeff([&](const mfem::Vector& x) { + double r = x.Norml2(); + if (r > R) return 0.0; + return rho0 * (1.0 - std::pow(r / R, 2.0)); + }); + + mfem::GridFunction rho(fem.H1_fes.get()); + rho.ProjectCoefficient(rho_coeff); + + fem.com = get_com(fem, rho); + fem.Q = compute_quadrupole_moment_tensor(fem, rho, fem.com); + + const auto phi = grav_potential(fem, args, rho); + + size_t failed = 0; + size_t num_tests = 0; + double max_rel_err = 0.0; + + // Phi(r) = 4*pi*G*rho0 * (r^2/6 - r^4/(20*R^2)) - pi*G*rho0*R^2 + const size_t check_count = std::min(50, fem.mesh->GetNE()); + for (int elemID = 0; elemID < check_count; ++elemID) { + auto* trans = fem.mesh->GetElementTransformation(elemID); + const auto& ip = mfem::IntRules.Get(trans->GetGeometryType(), 2).IntPoint(0); + trans->SetIntPoint(&ip); + + mfem::Vector x_phys; + fem.mapping->GetPhysicalPoint(*trans, ip, x_phys); + const double r = x_phys.Norml2(); + + if (r < R) { + num_tests++; + const double term1 = 4.0 * M_PI * G * rho0 * ((r*r / 6.0) - (std::pow(r, 4.0) / (20.0 * R*R))); + constexpr double term2 = M_PI * G * rho0 * R*R; + const double phi_analytic = term1 - term2; + + const double phi_fem = phi.GetValue(elemID, ip); + const double rel_err = std::abs(phi_fem - phi_analytic) / std::abs(phi_analytic); + + max_rel_err = std::max(max_rel_err, rel_err); + if (rel_err > 1e-4) ++failed; + } } + RANK_GUARD( + auto result_type = (failed == 0) ? TEST_RESULT_TYPE::SUCCESS : TEST_RESULT_TYPE::FAILURE; + std::println("{}", fmt_test_msg("Test Ferrers Inhomogeneous Potential", result_type, failed, num_tests)); + + if (result_type == TEST_RESULT_TYPE::FAILURE) { + std::println("\tFAILURE: max rel error: {:+0.2E}", max_rel_err); + } + ) } + +void test_force_continuity(FEM& fem, const Args& args) { + constexpr double rho0 = 1.0; + constexpr double R = RADIUS; + + mfem::GridFunction rho(fem.H1_fes.get()); + rho = rho0; + + fem.mapping->ResetDisplacement(); + update_stiffness_matrix(fem); + + const auto phi = grav_potential(fem, args, rho); + + size_t failed = 0; + size_t num_tests = 0; + double max_jump = 0.0; + + for (int i = 0; i < fem.mesh->GetNE(); i++) { + mfem::ElementTransformation *T = fem.mesh->GetElementTransformation(i); + const mfem::IntegrationRule &ir = mfem::IntRules.Get(T->GetGeometryType(), 2); + + for (int j = 0; j < ir.GetNPoints(); j++) { + const mfem::IntegrationPoint &ip = ir.IntPoint(j); + T->SetIntPoint(&ip); + + mfem::Vector x; + fem.mapping->GetPhysicalPoint(*T, ip, x); + const double r = x.Norml2(); + + // Check very close to the surface R + if (std::abs(r - R) < 0.05) { + num_tests++; + mfem::Vector grad_phi; + phi.GetGradient(*T, grad_phi); + + const double g_mag_fem = grad_phi.Norml2(); + const double total_mass = (4.0/3.0) * M_PI * std::pow(R, 3.0) * rho0; + const double g_mag_analytic = (r <= R) ? (4.0/3.0)*M_PI*G*rho0*r : (G*total_mass)/(r*r); + + const double rel_err = std::abs(g_mag_fem - g_mag_analytic) / g_mag_analytic; + max_jump = std::max(max_jump, rel_err); + if (rel_err > 5e-3) ++failed; // Gradients are usually 1 order less accurate than the solution + } + } + } + + RANK_GUARD( + std::println("{}", fmt_test_msg("Test Force Continuity", + (failed == 0) ? TEST_RESULT_TYPE::SUCCESS : TEST_RESULT_TYPE::FAILURE, failed, num_tests)); + ) +} + +void test_ferrers_ellipsoid_potential(FEM& fem, const Args& args) { + constexpr double a = 1.1 * RADIUS; + constexpr double b = 1.0 * RADIUS; + constexpr double c = 0.9 * RADIUS; + constexpr double rho0 = 1.0; + + mfem::GridFunction ferrers_disp(fem.Vec_H1_fes.get()); + mfem::VectorFunctionCoefficient disp_coeff(3, [&](const mfem::Vector& x, mfem::Vector& d) { + d.SetSize(3); + d(0) = (a/RADIUS - 1.0) * x(0); + d(1) = (b/RADIUS - 1.0) * x(1); + d(2) = (c/RADIUS - 1.0) * x(2); + }); + ferrers_disp.ProjectCoefficient(disp_coeff); + fem.mapping->SetDisplacement(ferrers_disp); + update_stiffness_matrix(fem); + + auto rho_func = [&](const mfem::Vector& x_phys) { + double m2 = std::pow(x_phys(0)/a, 2) + std::pow(x_phys(1)/b, 2) + std::pow(x_phys(2)/c, 2); + return (m2 < 1.0) ? rho0 * (1.0 - m2) : 0.0; + }; + PhysicalPositionFunctionCoefficient rho_coeff(*fem.mapping, rho_func); + + mfem::GridFunction rho(fem.H1_fes.get()); + rho.ProjectCoefficient(rho_coeff); + + fem.com = get_com(fem, rho); + fem.Q = compute_quadrupole_moment_tensor(fem, rho, fem.com); + const auto phi = grav_potential(fem, args, rho); + + auto calc_analytic_phi = [&](const mfem::Vector& x) { + auto integrand = [&](double theta) { + if (theta == 0.0) return 0.0; + if (theta >= M_PI / 2.0) return 1.0 / a; + + const double tan_t = std::tan(theta); + const double sec_t = 1.0 / std::cos(theta); + + const double u = a * a * tan_t * tan_t; + const double du_dtheta = 2.0 * a * a * tan_t * sec_t * sec_t; + + const double a2_u = a*a + u; + const double b2_u = b*b + u; + const double c2_u = c*c + u; + + const double delta = std::sqrt(a2_u * b2_u * c2_u); + const double m_u2 = (x(0)*x(0))/a2_u + (x(1)*x(1))/b2_u + (x(2)*x(2))/c2_u; + + const double val = 0.5 * std::pow(1.0 - m_u2, 2) / delta; + + return val * du_dtheta; + }; + + int n_steps = 1000; + double dtheta = (M_PI / 2.0) / n_steps; + double sum = integrand(0.0) + integrand(M_PI / 2.0); + for (int i = 1; i < n_steps; ++i) { + const double theta = i * dtheta; + sum += (i % 2 == 0 ? 2.0 : 4.0) * integrand(theta); + } + sum *= dtheta / 3.0; + + return -M_PI * G * a * b * c * rho0 * sum; + }; + + size_t failed = 0; + size_t num_tests = 0; + double max_rel_err = 0.0; + double total_rel_err = 0.0; + + const size_t check_count = std::min(50, fem.mesh->GetNE()); + for (int elemID = 0; elemID < check_count; ++elemID) { + if (fem.mesh->GetAttribute(elemID) != 1) continue; + + auto* trans = fem.mesh->GetElementTransformation(elemID); + const auto& ip = mfem::IntRules.Get(trans->GetGeometryType(), 2).IntPoint(0); + trans->SetIntPoint(&ip); + + mfem::Vector x_phys; + fem.mapping->GetPhysicalPoint(*trans, ip, x_phys); + + double phi_fem = phi.GetValue(elemID, ip); + double phi_analytic = calc_analytic_phi(x_phys); + + double rel_err = std::abs(phi_fem - phi_analytic) / std::abs(phi_analytic); + max_rel_err = std::max(max_rel_err, rel_err); + total_rel_err += rel_err; + num_tests++; + + if (rel_err > 1e-4) ++failed; + } + + RANK_GUARD( + auto result_type = (failed == 0) ? TEST_RESULT_TYPE::SUCCESS : TEST_RESULT_TYPE::FAILURE; + std::println("{}", fmt_test_msg("Test Ferrers Ellipsoid (n=1)", result_type, failed, num_tests)); + if (result_type == TEST_RESULT_TYPE::FAILURE) { + std::println("\tFAILURE: Max Rel Error: {:+0.2E}, Mean Rel Error: {:+0.2E}", max_rel_err, total_rel_err/static_cast(num_tests)); + } + ) +} + +void test_mass_conservation_constraint(const FEM& fem, const Args& args) { + constexpr double target_mass = MASS; + constexpr double R = RADIUS; + + fem.mapping->ResetDisplacement(); + + mfem::FunctionCoefficient rho_coeff([&](const mfem::Vector& x) { + double r = x.Norml2(); + return (r < R) ? (1.0 - (r/R)*(r/R)) : 0.0; + }); + + mfem::GridFunction rho(fem.H1_fes.get()); + rho.ProjectCoefficient(rho_coeff); + + auto enforce_mass = [&](mfem::GridFunction& gf) { + double current_m = domain_integrate_grid_function(fem, gf, Domains::STELLAR); + if (current_m > 0.0) { + gf *= (target_mass / current_m); + } + }; + + size_t failed = 0; + size_t num_tests = 0; + + enforce_mass(rho); + double mass_undeformed = domain_integrate_grid_function(fem, rho, Domains::STELLAR); + num_tests++; + if (std::abs(mass_undeformed - target_mass) / target_mass > 1e-12) { + failed++; + RANK_GUARD( + std::println("\tFAILURE (Undeformed Conservation):"); + std::println("\t\tExpected Mass: {:.6e}", target_mass); + std::println("\t\tActual Mass: {:.6e}", mass_undeformed); + std::println("\t\tRel Error: {:+0.2E}", std::abs(mass_undeformed - target_mass) / target_mass); + ) + } + + constexpr double a = 1.5 * RADIUS; + constexpr double b = 0.8 * RADIUS; + constexpr double c = 0.5 * RADIUS; + + mfem::GridFunction disp(fem.Vec_H1_fes.get()); + mfem::VectorFunctionCoefficient disp_coeff(3, [&](const mfem::Vector& x, mfem::Vector& d) { + d.SetSize(3); + d(0) = (a/R - 1.0) * x(0); + d(1) = (b/R - 1.0) * x(1); + d(2) = (c/R - 1.0) * x(2); + }); + disp.ProjectCoefficient(disp_coeff); + fem.mapping->SetDisplacement(disp); + + double expected_scaled_mass = target_mass * (a * b * c) / (R * R * R); + double mass_deformed_unscaled = domain_integrate_grid_function(fem, rho, Domains::STELLAR); + num_tests++; + if (std::abs(mass_deformed_unscaled - expected_scaled_mass) / expected_scaled_mass > 1e-10) { + failed++; + RANK_GUARD( + std::println("\tFAILURE (Geometric Scaling):"); + std::println("\t\tDisplacement mapping did not correctly scale the physical mass integral."); + std::println("\t\tExpected Scaled Mass: {:.6e}", expected_scaled_mass); + std::println("\t\tActual Scaled Mass: {:.6e}", mass_deformed_unscaled); + std::println("\t\tRel Error: {:+0.2E}", std::abs(mass_deformed_unscaled - expected_scaled_mass) / expected_scaled_mass); + ) + } + + enforce_mass(rho); + double mass_deformed_conserved = domain_integrate_grid_function(fem, rho, Domains::STELLAR); + num_tests++; + if (std::abs(mass_deformed_conserved - target_mass) / target_mass > 1e-12) { + failed++; + RANK_GUARD( + std::println("\tFAILURE (Deformed Conservation):"); + std::println("\t\tFailed to enforce mass constraint on the deformed geometry."); + std::println("\t\tExpected Mass: {:.6e}", target_mass); + std::println("\t\tActual Mass: {:.6e}", mass_deformed_conserved); + std::println("\t\tRel Error: {:+0.2E}", std::abs(mass_deformed_conserved - target_mass) / target_mass); + ) + } + + RANK_GUARD( + auto result_type = (failed == 0) ? TEST_RESULT_TYPE::SUCCESS : (failed < num_tests ? TEST_RESULT_TYPE::PARTIAL : TEST_RESULT_TYPE::FAILURE); + std::println("{}", fmt_test_msg("Test Mass Conservation Constraint", result_type, failed, num_tests)); + ) +} + +void test_xad_eos_derivative(const FEM& fem, const Args& args) { + constexpr double K = 2.5; + constexpr double gamma = 5.0 / 3.0; + + auto polytropic_eos = [&](const auto& rho) { + using std::pow; // Allow ADL to find xad::pow if necessary + return K * pow(rho, gamma); + }; + + fem.mapping->ResetDisplacement(); + mfem::FunctionCoefficient rho_coeff([&](const mfem::Vector& x) { + double r = x.Norml2(); + return std::max(0.1, 1.0 - (r / RADIUS)); + }); + mfem::GridFunction rho(fem.H1_fes.get()); + rho.ProjectCoefficient(rho_coeff); + + size_t failed = 0; + size_t num_tests = 0; + double max_rel_err_p = 0.0; + double max_rel_err_dp = 0.0; + + const size_t check_count = std::min(100, fem.mesh->GetNE()); + for (int elemID = 0; elemID < check_count; ++elemID) { + if (fem.mesh->GetAttribute(elemID) != 1) continue; // Only Core elements + + auto* trans = fem.mesh->GetElementTransformation(elemID); + const auto& ir = *fem.int_rule; + + for (int q = 0; q < ir.GetNPoints(); ++q) { + const auto& ip = ir.IntPoint(q); + trans->SetIntPoint(&ip); + + double rho_val = rho.GetValue(elemID, ip); + + double p_analytic = K * std::pow(rho_val, gamma); + double dp_drho_analytic = K * gamma * std::pow(rho_val, gamma - 1.0); + + xad::FReal rho_ad(rho_val, 1.0); + xad::FReal p_ad = polytropic_eos(rho_ad); + + const double p_xad = p_ad.getValue(); + const double dp_drho_xad = p_ad.getDerivative(); + + double rel_err_p = std::abs(p_xad - p_analytic) / std::abs(p_analytic); + double rel_err_dp = std::abs(dp_drho_xad - dp_drho_analytic) / std::abs(dp_drho_analytic); + + max_rel_err_p = std::max(max_rel_err_p, rel_err_p); + max_rel_err_dp = std::max(max_rel_err_dp, rel_err_dp); + + num_tests++; + if (rel_err_p > 1e-12 || rel_err_dp > 1e-12) { + failed++; + } + } + } + + RANK_GUARD( + const auto result_type = (failed == 0) ? TEST_RESULT_TYPE::SUCCESS : TEST_RESULT_TYPE::FAILURE; + std::println("{}", fmt_test_msg("Test XAD EOS Derivative", result_type, failed, num_tests)); + + if (result_type == TEST_RESULT_TYPE::FAILURE) { + std::println("\tFAILURE: Max Rel Error (P): {:+0.2E}, Max Rel Error (dP/drho): {:+0.2E}", + max_rel_err_p, max_rel_err_dp); + } + ) +} + +void test_domain_mapper_state_isolation(const FEM& fem) { + size_t failed = 0; + + mfem::GridFunction messy_disp(fem.Vec_H1_fes.get()); + mfem::VectorFunctionCoefficient disp_coeff(3, [](const mfem::Vector& x, mfem::Vector& d) { + d.SetSize(3); + d(0) = std::sin(x(0) * M_PI); + d(1) = std::cos(x(1) * M_PI); + d(2) = x(0) * x(1); + }); + messy_disp.ProjectCoefficient(disp_coeff); + fem.mapping->SetDisplacement(messy_disp); + + int elem_A = 0; + int elem_B = fem.mesh->GetNE() / 2; + + auto* trans_A = fem.mesh->GetElementTransformation(elem_A); + auto* trans_B = fem.mesh->GetElementTransformation(elem_B); + const auto& ip = mfem::IntRules.Get(trans_A->GetGeometryType(), 2).IntPoint(0); + + mfem::DenseMatrix J_A_baseline, J_B_baseline, J_A_test, J_B_test; + + trans_A->SetIntPoint(&ip); + fem.mapping->ComputeJacobian(*trans_A, J_A_baseline); + + trans_B->SetIntPoint(&ip); + fem.mapping->ComputeJacobian(*trans_B, J_B_baseline); + + for (int i = 0; i < 5; ++i) { + trans_A->SetIntPoint(&ip); + fem.mapping->ComputeJacobian(*trans_A, J_A_test); + + trans_B->SetIntPoint(&ip); + fem.mapping->ComputeJacobian(*trans_B, J_B_test); + + J_A_test.Add(-1.0, J_A_baseline); + if (J_A_test.MaxMaxNorm() > 1e-12) ++failed; + + J_B_test.Add(-1.0, J_B_baseline); + if (J_B_test.MaxMaxNorm() > 1e-12) ++failed; + } + + RANK_GUARD( + std::println("{}", fmt_test_msg("Test DomainMapper State Isolation", + (failed == 0) ? TEST_RESULT_TYPE::SUCCESS : TEST_RESULT_TYPE::FAILURE, + failed, 10)); + ) + + fem.mapping->ResetDisplacement(); +} + +void test_hydrostatic_zero_residual(FEM& fem, const Args& i_args) { + fem.mapping->ResetCacheStats(); + Args args = i_args; + + constexpr double R = RADIUS; + constexpr double k = M_PI / R; + constexpr double rho_c = 1.0; + constexpr double K_poly = (2.0 * G * R * R) / M_PI; + + constexpr double exact_mass = (4.0 * R * R * R * rho_c) / M_PI; + constexpr double exact_lambda = -(4.0 * G * R * R * rho_c) / M_PI; + + args.mass = exact_mass; + + auto eos_pressure = [&](const auto& rho_val, const auto& temp_val) { + using std::pow; + return K_poly * pow(rho_val, 2.0); + }; + + auto eos_enthalpy = [&](const auto& rho_val, const auto& temp_val) { + return 2.0 * K_poly * rho_val; + }; + + fem.mapping->ResetDisplacement(); + update_stiffness_matrix(fem); + + mfem::FunctionCoefficient rho_coeff([&](const mfem::Vector& x) -> double { + double r = x.Norml2(); + if (r < 1e-8) return rho_c; + if (r >= R) return 0.0; + return rho_c * std::sin(k * r) / (k * r); + }); + + mfem::GridFunction rho(fem.H1_fes.get()); + rho.ProjectCoefficient(rho_coeff); + + mfem::BlockVector state(fem.block_true_offsets); + state = 0.0; + + mfem::Vector& state_rho = state.GetBlock(0); + state_rho = rho; // TODO: How to this in a safer way which does not involve slicing a GridFunction into a vector + + state(fem.block_true_offsets[2]) = exact_lambda; + + mfem::BlockVector residual(fem.block_true_offsets); + residual = 0.0; + + ResidualOperator> coupled_operator(fem, args, eos_enthalpy, eos_pressure, 1.0); + coupled_operator.Mult(state, residual); + + mfem::Vector& r_rho = residual.GetBlock(0); + mfem::Vector& r_d = residual.GetBlock(1); + double r_lambda = residual(fem.block_true_offsets[2]); + + double force_scale = K_poly * rho_c * rho_c * R * R; + + double norm_rho = r_rho.Norml2() / force_scale; + double norm_d = r_d.Norml2(); + double abs_r_lambda = std::abs(r_lambda) / exact_mass; + + size_t failed = 0; + size_t num_tests = 3; + + if (norm_rho > 1e-4) ++failed; + if (norm_d > 1e-4) ++failed; + if (abs_r_lambda > 1e-5) ++failed; + + RANK_GUARD( + auto result_type = (failed == 0) ? TEST_RESULT_TYPE::SUCCESS : + (failed < num_tests ? TEST_RESULT_TYPE::PARTIAL : TEST_RESULT_TYPE::FAILURE); + + std::println("{}", fmt_test_msg("Test Hydrostatic Zero Residual (n=1)", result_type, failed, num_tests)); + + if (result_type != TEST_RESULT_TYPE::SUCCESS) { + std::println("\tFAILURE: The Residual Operator did not achieve equilibrium."); + std::println("\t\tDensity Rel Residual (H + Phi - Lambda): {:+0.2E}", norm_rho); + std::println("\t\tDisplacement Abs Residual (Stiffness - PBF): {:+0.2E}", norm_d); + std::println("\t\tLambda Rel Residual (Mass Conservation): {:+0.2E}", abs_r_lambda); + } + ) + + mfem::BlockVector perturbed_state(state); + mfem::Vector& perturbed_rho = perturbed_state.GetBlock(0); + mfem::Vector& perturbed_d = perturbed_state.GetBlock(1); + perturbed_rho *= 1.05; + + mfem::GridFunction p_d_gf(fem.Vec_H1_fes.get()); + mfem::VectorFunctionCoefficient pd_coeff(3, [&](const mfem::Vector& x, mfem::Vector& v) { + v.SetSize(3); + v(0) = 0.01 + x(0); + v(1) = 0.01 + x(1); + v(2) = 0.01 + x(2); + }); + p_d_gf.ProjectCoefficient(pd_coeff); + perturbed_d = p_d_gf; + + mfem::BlockVector perturbed_residual(fem.block_true_offsets); + perturbed_residual = 0.0; + coupled_operator.Mult(perturbed_state, perturbed_residual); + + mfem::Vector& p_r_rho = perturbed_residual.GetBlock(0); + mfem::Vector& p_r_d = perturbed_residual.GetBlock(1); + double p_r_lambda = perturbed_residual(fem.block_true_offsets[2]); + + double p_norm_rho = p_r_rho.Norml2() / force_scale; + double p_norm_d = p_r_d.Norml2(); + double p_abs_r_lambda = std::abs(p_r_lambda) / exact_mass; + + size_t perturb_failed = 0; + size_t perturb_tests = 3; + + if (p_norm_rho <= norm_rho) ++perturb_failed; + + if (p_norm_d <= norm_d) ++perturb_failed; + + if (p_abs_r_lambda <= abs_r_lambda) ++perturb_failed; + + RANK_GUARD( + auto perturb_result = (perturb_failed == 0) ? TEST_RESULT_TYPE::SUCCESS : + (perturb_failed < perturb_tests ? TEST_RESULT_TYPE::PARTIAL : TEST_RESULT_TYPE::FAILURE); + + std::println("{}", fmt_test_msg("Test Hydrostatic Perturbation Spike", perturb_result, perturb_failed, perturb_tests)); + + if (perturb_result != TEST_RESULT_TYPE::SUCCESS) { + std::println("\tFAILURE: Perturbing the state did not worsen the residual."); + std::println("\t\tDensity Res jump: {:.2e} -> {:.2e}", norm_rho, p_norm_rho); + std::println("\t\tDisp Res jump: {:.2e} -> {:.2e}", norm_d, p_norm_d); + std::println("\t\tLambda Res jump: {:.2e} -> {:.2e}", abs_r_lambda, p_abs_r_lambda); + }; + ) +} + +void test_single_newton_step(FEM& fem, const Args& i_args) { + Args args = i_args; + + constexpr double R = RADIUS; + constexpr double k = M_PI / R; + constexpr double rho_c = 1.0; + constexpr double K_poly = (2.0 * G * R * R) / M_PI; + + constexpr double exact_mass = (4.0 * R * R * R * rho_c) / M_PI; + constexpr double exact_lambda = -(4.0 * G * R * R * rho_c) / M_PI; + + args.mass = exact_mass; + + auto eos_pressure = [&](const auto& rho_val, const auto& temp_val) { + using std::pow; + return K_poly * pow(rho_val, 2.0); + }; + + auto eos_enthalpy = [&](const auto& rho_val, const auto& temp_val) { + return 2.0 * K_poly * rho_val; + }; + + fem.mapping->ResetDisplacement(); + mfem::FunctionCoefficient rho_coeff([&](const mfem::Vector& x) -> double { + double r = x.Norml2(); + if (r < 1e-8) return rho_c; + if (r >= R) return 0.0; + return rho_c * std::sin(k * r) / (k * r); + }); + + mfem::GridFunction rho(fem.H1_fes.get()); + rho.ProjectCoefficient(rho_coeff); + + mfem::BlockVector state(fem.block_true_offsets); + state = 0.0; + + mfem::Vector& state_rho = state.GetBlock(0); + state_rho = rho; // TODO: How to this in a safer way which does not involve slicing a GridFunction into a vector + + state(fem.block_true_offsets[2]) = exact_lambda; + + mfem::BlockVector residual(fem.block_true_offsets); + residual = 0.0; + + ResidualOperator> coupled_operator(fem, args, eos_enthalpy, eos_pressure, 1.0); + + mfem::BlockVector u_perturbed(state); + mfem::Vector& perturbed_rho = u_perturbed.GetBlock(0); + mfem::Vector& perturbed_d = u_perturbed.GetBlock(1); + perturbed_rho *= 1.05; + + mfem::GridFunction p_d_gf(fem.Vec_H1_fes.get()); + mfem::VectorFunctionCoefficient pd_coeff(3, [&](const mfem::Vector& x, mfem::Vector& v) { + v.SetSize(3); + v(0) = 0.01 * x(0); + v(1) = 0.01 * x(1); + v(2) = 0.01 * x(2); + }); + p_d_gf.ProjectCoefficient(pd_coeff); + perturbed_d = p_d_gf; + + mfem::BlockVector initial_res(fem.block_true_offsets); + coupled_operator.Mult(u_perturbed, initial_res); + double initial_norm = initial_res.Norml2(); + + // mfem::Operator& J = coupled_operator.GetGradient(u_perturbed); + JFNKOperator jfnk_j(coupled_operator, u_perturbed, initial_res); + + mfem::BlockVector du(fem.block_true_offsets); + du = 0.0; + + mfem::GMRESSolver solver; + solver.SetOperator(jfnk_j); + solver.SetAbsTol(1e-12); + solver.SetRelTol(1e-10); + solver.SetMaxIter(500); + solver.SetPrintLevel(0); + + mfem::Vector neg_res = initial_res; + neg_res *= -1.0; + + solver.Mult(neg_res, du); + + mfem::BlockVector u_new(u_perturbed); + u_new += du; + + // 5. Evaluate new residual + mfem::BlockVector final_res(fem.block_true_offsets); + coupled_operator.Mult(u_new, final_res); + double final_norm = final_res.Norml2(); + + bool success = (final_norm < initial_norm * 0.01); + + RANK_GUARD( + std::println("{}", fmt_test_msg("Test Single Newton Step", + success ? TEST_RESULT_TYPE::SUCCESS : TEST_RESULT_TYPE::FAILURE, + success ? 1 : 0, 1)); + + std::println("\tInitial Residual Norm: {:.6e}", initial_norm); + std::println("\tFinal Residual Norm: {:.6e}", final_norm); + std::println("\tConvergence Factor: {:.2e}", final_norm / initial_norm); + ) +} + //endregion +void run_grav_potential_timer(FEM& fem, const Args& args) { + mfem::GridFunction rho(fem.H1_fes.get()); + rho = 1.0; + fem.com = get_com(fem, rho); + fem.Q = compute_quadrupole_moment_tensor(fem, rho, fem.com); + + START_TIMER(grav_potential_cold_start); + auto phi_1 = grav_potential(fem, args, rho); + END_TIMER(grav_potential_cold_start); + rho *= 1.01; + START_TIMER(grav_potential_warm_start); + auto phi_2 = grav_potential(fem, args, rho, true); + END_TIMER(grav_potential_warm_start); +} + int main(int argc, char** argv) { + mfem::Mpi::Init(argc, argv); + + std::string device_config = "cpu"; + mfem::Device device(device_config); + + const int myid = mfem::Mpi::WorldRank(); + const int num_procs = mfem::Mpi::WorldSize(); + + if (myid == 0) { + std::println("Starting MFEM run on {} processors.", num_procs); + } + Args args; CLI::App app{"Mapped Coordinates XAD-Enabled Non-Linear Solver"}; @@ -2026,29 +3083,41 @@ int main(int argc, char** argv) { app.add_option("--index", args.index)->default_val(1); app.add_option("--mass", args.mass)->default_val(MASS); app.add_option("--surface-pressure-scalar", args.c)->default_val(1e-4); + app.add_option("-q,--quad-boost", args.quad_boost)->default_val(0); args.r.enabled = false; - args.p.tol = 1e-6; + args.p.rtol = 1e-8; + args.p.atol = 1e-7; args.p.max_iters = 1000; CLI11_PARSE(app, argc, argv); - FEM fem = setup_fem(args.mesh_file, true); - // fem.mesh->UniformRefinement(); + FEM fem = setup_fem(args.mesh_file, args); + // run_grav_potential_timer(fem, args); + RUN_TEST("Mesh Loading", test_mesh_load(fem)); + RUN_TEST("Test Reference Coordinates", test_ref_coord_storage(fem)); + RUN_TEST("Test Reference Volume Integral", test_reference_volume_integral(fem)); + RUN_TEST("Test Spherically Symmetric Center of Mass", test_spherically_symmetric_com(fem)); + RUN_TEST("Test COM variance to displacement", test_com_variance_to_displacement(fem)); + RUN_TEST("Test Volume Invariance to Displacement", test_volume_invariance_to_displacement(fem)) + RUN_TEST("Test Volume of Ellipsoid Deformation", test_volume_ellipsoid_deformation(fem)); - test_mesh_load(fem); - test_ref_coord_storage(fem); - test_reference_volume_integral(fem); - test_spherically_symmetric_com(fem); + RUN_TEST("Test Uniform Potential", test_uniform_potential(fem, args)); + RUN_TEST("Test Ellipsoidal Potential", test_ellipsoidal_potential(fem, args)); + RUN_TEST("Test Ferrers Sphere Potential", test_ferrers_sphere_potential(fem, args)); + RUN_TEST("Test Ferrers Ellipsoid Potential", test_ferrers_ellipsoid_potential(fem, args)); - test_com_variance_to_displacement(fem); - test_volume_invariance_to_displacement(fem); - test_volume_ellipsoid_deformation(fem); - // test_kelvin_jacobian(fem); - test_uniform_potential(fem, args); - test_ellipsoidal_potential(fem, args); // Note that this test is currently predicated on an analytic solution for the surface boundary potential + RUN_TEST("Test Mass Conservation Constraint", test_mass_conservation_constraint(fem, args)); + RUN_TEST("Test XAD EOS Derivative", test_xad_eos_derivative(fem, args)); + RUN_TEST("Test Force Continuity", test_force_continuity(fem, args)); - CoupledState state(fem); + RUN_TEST("Test Domain Mapper State Isolation", test_domain_mapper_state_isolation(fem)); + + RUN_TEST("Test Hydrostatic Zero Residuals", test_hydrostatic_zero_residual(fem, args)); + + // test_single_newton_step(fem, args); + + // CoupledState state(fem); // typedef xad::AReal adouble; // diff --git a/meson.build b/meson.build new file mode 100644 index 0000000..fd9e6c2 --- /dev/null +++ b/meson.build @@ -0,0 +1,14 @@ +project('RBPoly', ['c', 'cpp'], + version : '0.1.0', + default_options : ['cpp_std=c++23'] +) + +mfem_dep = dependency('mfem', required : true) +xad_dep = dependency('XAD', method : 'cmake', modules : ['XAD::xad']) +suitesparse_dep = dependency('UMFPACK', required: true) +mpi_dep = dependency('MPI', required: true) +hypre_dep = dependency('hypre', required: true) + +#executable('fp', 'main.cpp', dependencies : mfem_dep) +#executable('free_energy', 'free_energy.cpp', dependencies : mfem_dep) +executable('mapping', 'mapping.cpp', dependencies : [mfem_dep, xad_dep, suitesparse_dep, mpi_dep, hypre_dep]) \ No newline at end of file