diff --git a/src/smith/differentiable_numerics/nonlinear_block_solver.cpp b/src/smith/differentiable_numerics/nonlinear_block_solver.cpp index c5ac92d3b..d62f554fb 100644 --- a/src/smith/differentiable_numerics/nonlinear_block_solver.cpp +++ b/src/smith/differentiable_numerics/nonlinear_block_solver.cpp @@ -9,6 +9,8 @@ #include "smith/physics/state/finite_element_dual.hpp" #include "smith/numerics/stdfunction_operator.hpp" #include "smith/numerics/equation_solver.hpp" +#include "smith/numerics/block_preconditioner.hpp" +#include "smith/numerics/solver_with_preconditioner.hpp" #include "smith/physics/boundary_conditions/boundary_condition_manager.hpp" #include "smith/physics/mesh.hpp" #include "mfem.hpp" @@ -95,9 +97,57 @@ void NonlinearBlockSolver::completeSetup(const std::vector& us) const auto* mfem_solver = &nonlinear_solver_->preconditioner(); - auto* amg_prec = dynamic_cast(mfem_solver); - if (amg_prec) { - amg_prec->SetSystemsOptions(best->space().GetVDim(), smith::ordering == mfem::Ordering::byNODES); + auto configure_amg = [](mfem::Solver* s, int vdim) { + if (!s) { + return; + } + if (auto* amg = dynamic_cast(s)) { + amg->SetSystemsOptions(vdim, smith::ordering == mfem::Ordering::byNODES); + return; + } + if (auto* wrapped = dynamic_cast(s)) { + if (auto* amg = dynamic_cast(wrapped->preconditioner())) { + amg->SetSystemsOptions(vdim, smith::ordering == mfem::Ordering::byNODES); + } + } + }; + + // Top-level AMG preconditioner + configure_amg(mfem_solver, best->space().GetVDim()); + + // Block preconditioners: configure per-block AMG with each block's vdim. + struct BlockSolverAccessor { + int n = 0; + mfem::Solver* (*get)(void*, int) = nullptr; + void* self = nullptr; + }; + + const BlockSolverAccessor acc = [&]() { + if (auto* p = dynamic_cast(mfem_solver)) { + return BlockSolverAccessor{p->numSubSolvers(), static_cast([](void* self, int i) { + return static_cast(self)->subSolver(i); + }), + p}; + } + if (auto* p = dynamic_cast(mfem_solver)) { + return BlockSolverAccessor{p->numSubSolvers(), static_cast([](void* self, int i) { + return static_cast(self)->subSolver(i); + }), + p}; + } + if (auto* p = dynamic_cast(mfem_solver)) { + return BlockSolverAccessor{p->numSubSolvers(), static_cast([](void* self, int i) { + return static_cast(self)->subSolver(i); + }), + p}; + } + return BlockSolverAccessor{}; + }(); + + if (acc.get) { + for (int i = 0; i < acc.n && static_cast(i) < us.size(); ++i) { + configure_amg(acc.get(acc.self, i), us[static_cast(i)]->space().GetVDim()); + } } #ifdef SMITH_USE_PETSC diff --git a/src/smith/differentiable_numerics/tests/CMakeLists.txt b/src/smith/differentiable_numerics/tests/CMakeLists.txt index bbe95cbd5..680c5d077 100644 --- a/src/smith/differentiable_numerics/tests/CMakeLists.txt +++ b/src/smith/differentiable_numerics/tests/CMakeLists.txt @@ -16,6 +16,8 @@ set(differentiable_numerics_test_source test_thermal_static.cpp test_solid_static_with_internal_vars.cpp test_multiphysics_time_integrator.cpp + test_mixed_poisson.cpp + test_amg_vdim_block_preconditioner.cpp ) smith_add_tests( SOURCES ${differentiable_numerics_test_source} diff --git a/src/smith/differentiable_numerics/tests/test_amg_vdim_block_preconditioner.cpp b/src/smith/differentiable_numerics/tests/test_amg_vdim_block_preconditioner.cpp new file mode 100644 index 000000000..6256168b9 --- /dev/null +++ b/src/smith/differentiable_numerics/tests/test_amg_vdim_block_preconditioner.cpp @@ -0,0 +1,119 @@ +// Copyright (c) Lawrence Livermore National Security, LLC and +// other Smith Project Developers. See the top-level LICENSE file for +// details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +#include + +#include + +#include "mfem.hpp" +#include + +#include "smith/infrastructure/application_manager.hpp" +#include "smith/differentiable_numerics/nonlinear_block_solver.hpp" +#include "smith/numerics/equation_solver.hpp" +#include "smith/numerics/solver_config.hpp" +#include "smith/numerics/block_preconditioner.hpp" +#include "smith/numerics/solver_with_preconditioner.hpp" +#include "smith/physics/state/finite_element_state.hpp" + +namespace smith { +namespace { + +HYPRE_Int GetBoomerAMGNumFunctions(const mfem::HypreBoomerAMG& amg) +{ + HYPRE_Solver h = amg; // mfem::HypreBoomerAMG has operator HYPRE_Solver() + HYPRE_Int dim = -1; + HYPRE_BoomerAMGGetNumFunctions(h, &dim); + return dim; +} + +TEST(AMGVdimSetup, BlockDiagonalPreconditionerSetsPerBlockVdim) +{ + MPI_Comm comm = MPI_COMM_WORLD; + + mfem::Mesh mesh = mfem::Mesh::MakeCartesian2D(1, 1, mfem::Element::QUADRILATERAL, 1, 1.0, 1.0); + mfem::ParMesh pmesh(comm, mesh); + mesh.Clear(); + + const int order = 1; + mfem::H1_FECollection fec(order, pmesh.Dimension()); + + // Block 0: vdim=2 + mfem::ParFiniteElementSpace fes_v2(&pmesh, &fec, 2, smith::ordering); + // Block 1: vdim=1 + mfem::ParFiniteElementSpace fes_v1(&pmesh, &fec, 1, smith::ordering); + + auto u0 = std::make_shared(fes_v2, "u0"); + auto u1 = std::make_shared(fes_v1, "u1"); + std::vector us{u0, u1}; + + NonlinearSolverOptions nonlin_opts; + nonlin_opts.nonlin_solver = NonlinearSolver::Newton; + nonlin_opts.max_iterations = 1; + nonlin_opts.print_level = 0; + nonlin_opts.absolute_tol = 1e-12; + nonlin_opts.relative_tol = 1e-10; + + LinearSolverOptions lin_opts; + lin_opts.linear_solver = LinearSolver::GMRES; + lin_opts.preconditioner = Preconditioner::BlockDiagonal; + lin_opts.max_iterations = 1; + lin_opts.print_level = 0; + lin_opts.preconditioner_print_level = 0; + + LinearSolverOptions sub; + sub.linear_solver = LinearSolver::GMRES; + sub.preconditioner = Preconditioner::HypreAMG; + sub.max_iterations = 1; + sub.relative_tol = 0.99; + sub.absolute_tol = 0.0; + sub.print_level = 0; + sub.preconditioner_print_level = 0; + + lin_opts.sub_block_linear_solver_options.push_back(sub); + lin_opts.sub_block_linear_solver_options.push_back(sub); + + auto eq = std::make_unique(nonlin_opts, lin_opts, comm); + mfem::Solver* top_prec = &eq->preconditioner(); + + NonlinearBlockSolver nbs(std::move(eq), comm, nonlin_opts.absolute_tol, nonlin_opts.relative_tol, nonlin_opts, + lin_opts); + nbs.completeSetup(us); + + auto* block_diag = dynamic_cast(top_prec); + ASSERT_NE(block_diag, nullptr); + ASSERT_GE(block_diag->numSubSolvers(), 2); + + // Block 0 + { + auto* wrapped = dynamic_cast(block_diag->subSolver(0)); + ASSERT_NE(wrapped, nullptr); + + auto* amg = dynamic_cast(wrapped->preconditioner()); + ASSERT_NE(amg, nullptr); + EXPECT_EQ(GetBoomerAMGNumFunctions(*amg), 2); + } + + // Block 1 + { + auto* wrapped = dynamic_cast(block_diag->subSolver(1)); + ASSERT_NE(wrapped, nullptr); + + auto* amg = dynamic_cast(wrapped->preconditioner()); + ASSERT_NE(amg, nullptr); + EXPECT_EQ(GetBoomerAMGNumFunctions(*amg), 1); + } +} + +} // namespace +} // namespace smith + +int main(int argc, char** argv) +{ + ::testing::InitGoogleTest(&argc, argv); + smith::ApplicationManager applicationManager(argc, argv); + return RUN_ALL_TESTS(); +} diff --git a/src/smith/differentiable_numerics/tests/test_mixed_poisson.cpp b/src/smith/differentiable_numerics/tests/test_mixed_poisson.cpp new file mode 100644 index 000000000..45a021449 --- /dev/null +++ b/src/smith/differentiable_numerics/tests/test_mixed_poisson.cpp @@ -0,0 +1,281 @@ +#include + +#include "smith/infrastructure/application_manager.hpp" +#include "smith/numerics/equation_solver.hpp" +#include "smith/numerics/solver_config.hpp" +#include "smith/mesh_utils/mesh_utils.hpp" +#include "smith/physics/state/state_manager.hpp" +#include "smith/physics/boundary_conditions/boundary_condition_manager.hpp" +#include "smith/physics/functional_weak_form.hpp" +#include "smith/differentiable_numerics/field_state.hpp" +#include "smith/differentiable_numerics/state_advancer.hpp" +#include "smith/smith_config.hpp" +#include "smith/differentiable_numerics/nonlinear_block_solver.hpp" +#include "smith/differentiable_numerics/nonlinear_solve.hpp" +#include "smith/differentiable_numerics/paraview_writer.hpp" +#include "smith/numerics/block_preconditioner.hpp" + +#include "gretl/data_store.hpp" +#include "gretl/wang_checkpoint_strategy.hpp" + +using namespace smith; + +using ShapeDispSpace = H1<1, 2>; +// Taylor-Hood elements +using Space = H1<1>; +using VectorSpace = H1<2, 2>; + +enum class BlockSolverType +{ + Direct, + Iterative, + BoomerAMG +}; +enum class BlockPrecondType +{ + SchurLower, + SchurUpper, + SchurDiag, + SchurFull, + SchurFullCustom +}; + +struct BlockTestParams { + BlockSolverType solver_type; + BlockPrecondType precond_type; +}; + +std::string BlockParamNameGenerator(const ::testing::TestParamInfo& info) +{ + auto solver_to_str = [](BlockSolverType t) { + switch (t) { + case BlockSolverType::Direct: + return "Direct"; + case BlockSolverType::Iterative: + return "Iterative"; + case BlockSolverType::BoomerAMG: + return "BoomerAMG"; + } + return "Unknown"; + }; + auto precond_to_str = [](BlockPrecondType t) { + switch (t) { + case BlockPrecondType::SchurLower: + return "SchurLower"; + case BlockPrecondType::SchurUpper: + return "SchurUpper"; + case BlockPrecondType::SchurDiag: + return "SchurDiag"; + case BlockPrecondType::SchurFull: + return "SchurFull"; + case BlockPrecondType::SchurFullCustom: + return "SchurFullCustom"; + } + return "Unknown"; + }; + return std::string(solver_to_str(info.param.solver_type)) + "_" + precond_to_str(info.param.precond_type); +} + +class MeshFixture : public testing::Test { + protected: + double length = 1.0; + double width = 1.0; + int num_elements_x = 32; + int num_elements_y = 32; + double elem_size = length / num_elements_x; + + axom::sidre::DataStore datastore; + std::shared_ptr mesh; + + void SetUp() override + { + smith::StateManager::initialize(datastore, "mixed_poisson"); + + MPI_Barrier(MPI_COMM_WORLD); + int serial_refinement = 4; + int parallel_refinement = 0; + + std::string filename = SMITH_REPO_DIR "/data/meshes/square_attribute.mesh"; + + const std::string meshtag = "mesh"; + mesh = std::make_shared(smith::buildMeshFromFile(filename), meshtag, serial_refinement, + parallel_refinement); + } +}; + +class BlockPreconditionerTest : public MeshFixture, public ::testing::WithParamInterface {}; + +TEST_P(BlockPreconditionerTest, BlockSolve) +{ + const auto& test_params = GetParam(); + + std::string physics_name = "mixed_poisson"; + auto graph = std::make_shared(std::make_unique(100)); + auto shape_disp = createFieldState(*graph, ShapeDispSpace{}, physics_name + "_shape_displacement", mesh->tag()); + auto flux = createFieldState(*graph, VectorSpace{}, physics_name + "_flux", mesh->tag()); + auto potential = createFieldState(*graph, Space{}, physics_name + "_potential", mesh->tag()); + smith::FunctionalWeakForm<2, VectorSpace, smith::Parameters> con_form( + "constitutive_eqn", mesh, space(flux), spaces({flux, potential})); + smith::FunctionalWeakForm<2, Space, smith::Parameters> bal_form( + "balance_eqn", mesh, space(potential), spaces({flux, potential})); + + con_form.addBodyIntegral(DependsOn<0, 1>{}, mesh->entireBodyName(), + [](double /* t */, auto /* x */, auto SIGMA, auto U) { + auto sigma = get(SIGMA); + auto u = get(U); + // Need to wrap u in a tensor to convert grad(test_function) to . div(test_function) + using Scalar = decltype(u); + smith::tensor u_{}; + u_[0][0] = u; + u_[1][1] = u; + return smith::tuple{sigma, -u_}; + }); + + bal_form.addBodyIntegral(DependsOn<0>{}, mesh->entireBodyName(), [](double /* t */, auto X, auto SIGMA) { + const auto& x = get(X); + auto div_sigma = smith::tr(get(SIGMA)); + double pi = M_PI; + auto f = 2.0 * pi * pi * sin(pi * x[0]) * sin(pi * x[1]); + return smith::tuple{-f + div_sigma, smith::zero{}}; + }); + // u_exact = sin(M_PI * x(0)) * sin(M_PI * x(1)); + // sigma_exact + // pi * cos(pi * x(0)) * sin(pi * x(1)) + // pi * sin(pi * x(0)) * cos(pi * x(1)); + + auto flux_bc_manager = std::make_shared(mesh->mfemParMesh()); + auto potential_bc_manager = std::make_shared(mesh->mfemParMesh()); + + auto zero_bcs = std::make_shared([](const mfem::Vector&) { return 0.0; }); + potential_bc_manager->addEssential(std::set{1, 2, 3, 4}, zero_bcs, space(potential), 0); + + // Block Preconditioner Options + smith::LinearSolverOptions linear_options; + linear_options.linear_solver = smith::LinearSolver::GMRES; + linear_options.relative_tol = 1.0e-10; + linear_options.absolute_tol = 1.0e-14; + linear_options.max_iterations = 100; + linear_options.print_level = 1; + + // Parameter sweep: construct solvers according to test parameters + if (test_params.solver_type == BlockSolverType::Direct) { + smith::LinearSolverOptions direct_solver_options{.linear_solver = smith::LinearSolver::SuperLU}; + linear_options.sub_block_linear_solver_options.push_back(direct_solver_options); + linear_options.sub_block_linear_solver_options.push_back(direct_solver_options); + } else if (test_params.solver_type == BlockSolverType::Iterative) { + smith::LinearSolverOptions iter_solver_options = {.linear_solver = smith::LinearSolver::GMRES, + .preconditioner = smith::Preconditioner::HypreAMG, + .relative_tol = 1.0e-3, + .absolute_tol = 1.0e-6, + .max_iterations = 100, + .print_level = 1}; + linear_options.sub_block_linear_solver_options.push_back(iter_solver_options); + linear_options.sub_block_linear_solver_options.push_back(iter_solver_options); + } else if (test_params.solver_type == BlockSolverType::BoomerAMG) { + smith::LinearSolverOptions amg_solver_options; + amg_solver_options.linear_solver = smith::LinearSolver::None; + amg_solver_options.preconditioner = smith::Preconditioner::HypreAMG; + linear_options.sub_block_linear_solver_options.push_back(amg_solver_options); + linear_options.sub_block_linear_solver_options.push_back(amg_solver_options); + } + + // Need these here for the custom operator + auto time = graph->create_state(0.0); + auto dt = graph->create_state(0.025); + size_t cycle = 0; + std::vector params; + auto& flux_params = params; + auto& potential_params = params; + std::vector con_arguments{flux, potential}; + std::vector bal_arguments{flux, potential}; + + switch (test_params.precond_type) { + case BlockPrecondType::SchurLower: + linear_options.preconditioner = smith::Preconditioner::BlockSchur; + linear_options.block_schur_type = smith::BlockSchurType::Lower; + break; + case BlockPrecondType::SchurUpper: + linear_options.preconditioner = smith::Preconditioner::BlockSchur; + linear_options.block_schur_type = smith::BlockSchurType::Upper; + break; + case BlockPrecondType::SchurDiag: + linear_options.preconditioner = smith::Preconditioner::BlockSchur; + linear_options.block_schur_type = smith::BlockSchurType::Diagonal; + break; + case BlockPrecondType::SchurFull: + linear_options.preconditioner = smith::Preconditioner::BlockSchur; + linear_options.block_schur_type = smith::BlockSchurType::Full; + break; + case BlockPrecondType::SchurFullCustom: + linear_options.preconditioner = smith::Preconditioner::BlockSchur; + linear_options.block_schur_type = smith::BlockSchurType::Full; + /// linear_options.schur_approx_type = smith::BlockSchurType::Custom; + break; + } + + smith::NonlinearSolverOptions nonlin_opts; + nonlin_opts.nonlin_solver = smith::NonlinearSolver::Newton; + nonlin_opts.relative_tol = linear_options.relative_tol; + nonlin_opts.absolute_tol = linear_options.absolute_tol; + nonlin_opts.max_iterations = 1; + nonlin_opts.print_level = linear_options.print_level; + + auto nonlinear_block_solver = smith::buildNonlinearBlockSolver(nonlin_opts, linear_options, *mesh); + + auto sols = block_solve({&con_form, &bal_form}, {{0, 1}, {0, 1}}, shape_disp, {con_arguments, bal_arguments}, + {flux_params, potential_params}, smith::TimeInfo(time.get(), dt.get(), cycle), + nonlinear_block_solver.get(), {flux_bc_manager.get(), potential_bc_manager.get()}); + + auto pv_writer = smith::createParaviewWriter(*mesh, sols, physics_name); + pv_writer.write(0, 0.0, sols); + + // Calculate error versus exact solution + auto u_exact_fun = [](const mfem::Vector& X) -> double { + double x = X(0), y = X(1); + return std::sin(M_PI * x) * std::sin(M_PI * y); + }; + mfem::FunctionCoefficient u_exact_coef(u_exact_fun); + + auto sigma_exact_fun = [](const mfem::Vector& X, mfem::Vector& S) { + double x = X(0), y = X(1); + S.SetSize(2); + S(0) = -M_PI * std::cos(M_PI * x) * std::sin(M_PI * y); + S(1) = -M_PI * std::sin(M_PI * x) * std::cos(M_PI * y); + }; + mfem::VectorFunctionCoefficient sigma_exact_coef(2, sigma_exact_fun); + double u_err = smith::computeL2Error(*sols[1].get(), u_exact_coef); + double sigma_err = smith::computeL2Error(*sols[0].get(), sigma_exact_coef); + + EXPECT_LT(u_err, 5e-3); + EXPECT_LT(sigma_err, 5e-2); + SUCCEED(); +} + +INSTANTIATE_TEST_SUITE_P(BlockPrecondSweep, BlockPreconditionerTest, + ::testing::Values( + // Direct solvers + BlockTestParams{BlockSolverType::Direct, BlockPrecondType::SchurLower}, + BlockTestParams{BlockSolverType::Direct, BlockPrecondType::SchurUpper}, + BlockTestParams{BlockSolverType::Direct, BlockPrecondType::SchurDiag}, + BlockTestParams{BlockSolverType::Direct, BlockPrecondType::SchurFull}, + + // Iterative solvers + BlockTestParams{BlockSolverType::Iterative, BlockPrecondType::SchurLower}, + BlockTestParams{BlockSolverType::Iterative, BlockPrecondType::SchurUpper}, + BlockTestParams{BlockSolverType::Iterative, BlockPrecondType::SchurDiag}, + BlockTestParams{BlockSolverType::Iterative, BlockPrecondType::SchurFull}, + + // BoomerAMG solvers + BlockTestParams{BlockSolverType::BoomerAMG, BlockPrecondType::SchurLower}, + BlockTestParams{BlockSolverType::BoomerAMG, BlockPrecondType::SchurUpper}, + BlockTestParams{BlockSolverType::BoomerAMG, BlockPrecondType::SchurDiag}, + BlockTestParams{BlockSolverType::BoomerAMG, BlockPrecondType::SchurFull}, + BlockTestParams{BlockSolverType::BoomerAMG, BlockPrecondType::SchurFullCustom}), + BlockParamNameGenerator); + +int main(int argc, char* argv[]) +{ + ::testing::InitGoogleTest(&argc, argv); + smith::ApplicationManager applicationManager(argc, argv); + return RUN_ALL_TESTS(); +} diff --git a/src/smith/differentiable_numerics/tests/test_porous_heat_sink.cpp b/src/smith/differentiable_numerics/tests/test_porous_heat_sink.cpp index 394401693..83b856d8c 100644 --- a/src/smith/differentiable_numerics/tests/test_porous_heat_sink.cpp +++ b/src/smith/differentiable_numerics/tests/test_porous_heat_sink.cpp @@ -234,9 +234,8 @@ TEST_P(BlockPreconditionerTest, BlockSolve) linear_options.sub_block_linear_solver_options.push_back(iter_solver_options); } else if (test_params.solver_type == BlockSolverType::BoomerAMG) { smith::LinearSolverOptions amg_solver_options; - amg_solver_options.linear_solver = smith::LinearSolver::CG; + amg_solver_options.linear_solver = smith::LinearSolver::None; amg_solver_options.preconditioner = smith::Preconditioner::HypreAMG; - amg_solver_options.max_iterations = 1; // Since it's a preconditioner-only analog linear_options.sub_block_linear_solver_options.push_back(amg_solver_options); linear_options.sub_block_linear_solver_options.push_back(amg_solver_options); } diff --git a/src/smith/numerics/CMakeLists.txt b/src/smith/numerics/CMakeLists.txt index cc16a2ea0..1a01d3673 100644 --- a/src/smith/numerics/CMakeLists.txt +++ b/src/smith/numerics/CMakeLists.txt @@ -11,6 +11,7 @@ set(numerics_headers nonlinear_convergence.hpp odes.hpp solver_config.hpp + solver_with_preconditioner.hpp stdfunction_operator.hpp petsc_solvers.hpp trust_region_solver.hpp diff --git a/src/smith/numerics/block_preconditioner.hpp b/src/smith/numerics/block_preconditioner.hpp index 3b0bac38e..b56ae3d8b 100644 --- a/src/smith/numerics/block_preconditioner.hpp +++ b/src/smith/numerics/block_preconditioner.hpp @@ -57,6 +57,20 @@ class BlockDiagonalPreconditioner : public mfem::Solver { virtual ~BlockDiagonalPreconditioner(); + /** @brief Return the number of sub-solvers owned by this preconditioner. */ + int numSubSolvers() const { return num_blocks_; } + + /** + * @brief Access a sub-solver by index. + * @param i Sub-solver index in [0, numSubSolvers()). + * @return Pointer to the requested sub-solver (owned by this object). + */ + mfem::Solver* subSolver(int i) const + { + MFEM_VERIFY(i >= 0 && i < num_blocks_, "BlockDiagonalPreconditioner::subSolver index out of range"); + return mfem_solvers_[static_cast(i)].get(); + } + private: // Offsets for extracting block vector segments, populated by SetOperator(). mfem::Array block_offsets_; @@ -129,6 +143,20 @@ class BlockTriangularPreconditioner : public mfem::Solver { virtual ~BlockTriangularPreconditioner(); + /** @brief Return the number of sub-solvers owned by this preconditioner. */ + int numSubSolvers() const { return num_blocks_; } + + /** + * @brief Access a sub-solver by index. + * @param i Sub-solver index in [0, numSubSolvers()). + * @return Pointer to the requested sub-solver (owned by this object). + */ + mfem::Solver* subSolver(int i) const + { + MFEM_VERIFY(i >= 0 && i < num_blocks_, "BlockTriangularPreconditioner::subSolver index out of range"); + return mfem_solvers_[static_cast(i)].get(); + } + private: // Offsets for extracting block vector segments, populated by SetOperator(). mfem::Array block_offsets_; @@ -233,6 +261,21 @@ class BlockSchurPreconditioner : public mfem::Solver { virtual ~BlockSchurPreconditioner(); + /** @brief Return the number of sub-solvers owned by this preconditioner. */ + int numSubSolvers() const { return static_cast(mfem_solvers_.size()); } + + /** + * @brief Access a sub-solver by index. + * @param i Sub-solver index in [0, numSubSolvers()). + * @return Pointer to the requested sub-solver (owned by this object). + */ + mfem::Solver* subSolver(int i) const + { + MFEM_VERIFY(i >= 0 && i < static_cast(mfem_solvers_.size()), + "BlockSchurPreconditioner::subSolver index out of range"); + return mfem_solvers_[static_cast(i)].get(); + } + private: // Offsets for extracting block vector segments, populated by SetOperator(). mfem::Array block_offsets_; diff --git a/src/smith/numerics/equation_solver.cpp b/src/smith/numerics/equation_solver.cpp index 70ebd6a41..ea7431598 100644 --- a/src/smith/numerics/equation_solver.cpp +++ b/src/smith/numerics/equation_solver.cpp @@ -6,6 +6,7 @@ #include "smith/numerics/equation_solver.hpp" #include "smith/numerics/block_preconditioner.hpp" +#include "smith/numerics/solver_with_preconditioner.hpp" #include #include @@ -25,27 +26,35 @@ namespace smith { namespace { - -class SolverWithPreconditioner : public mfem::Solver { +/** + * @brief Simple solver wrapper that only applies a preconditioner. + */ +class PreconditionerOnlySolver : public mfem::IterativeSolver { public: - SolverWithPreconditioner(std::unique_ptr linear_solver, std::unique_ptr preconditioner) - : linear_solver_(std::move(linear_solver)), preconditioner_(std::move(preconditioner)) + PreconditionerOnlySolver(MPI_Comm mpi_comm) : mfem::IterativeSolver(mpi_comm) {} + + /// @overload + void Mult(const mfem::Vector& x, mfem::Vector& y) const override { - SLIC_ERROR_IF(!linear_solver_, "SolverWithPreconditioner requires a non-null linear solver"); + if (prec) { + prec->Mult(x, y); + } else { + y = x; + } } + /// @overload void SetOperator(const mfem::Operator& op) override { - height = op.Height(); + if (prec) { + prec->SetOperator(op); + } width = op.Width(); - linear_solver_->SetOperator(op); + height = op.Height(); } - void Mult(const mfem::Vector& x, mfem::Vector& y) const override { linear_solver_->Mult(x, y); } - private: - std::unique_ptr linear_solver_; - std::unique_ptr preconditioner_; + // Note: mfem::IterativeSolver already has a 'prec' member (mfem::Solver*) }; bool preconditionerSupportsBlockOperator(Preconditioner preconditioner) @@ -74,6 +83,7 @@ bool linearSolverSupportsBlockOperator(LinearSolver linear_solver) case LinearSolver::PetscCG: case LinearSolver::PetscGMRES: #endif + case LinearSolver::None: return true; default: return false; @@ -1355,6 +1365,9 @@ std::pair, std::unique_ptr> buildLin exit(1); break; #endif + case LinearSolver::None: + iter_lin_solver = std::make_unique(comm); + break; default: SLIC_ERROR_ROOT("Linear solver type not recognized."); exit(1); @@ -1485,7 +1498,7 @@ std::unique_ptr buildPreconditioner(LinearSolverOptions linear_opt std::vector> inner_solvers; for (const auto& opt : linear_opts.sub_block_linear_solver_options) { auto [lin, prec] = buildLinearSolverAndPreconditioner(opt, comm); - inner_solvers.push_back(std::make_unique(std::move(lin), std::move(prec))); + inner_solvers.push_back(std::make_unique(std::move(lin), std::move(prec))); } if (preconditioner == Preconditioner::BlockDiagonal) { diff --git a/src/smith/numerics/solver_config.hpp b/src/smith/numerics/solver_config.hpp index 6f30ef8d3..6cfd77c80 100644 --- a/src/smith/numerics/solver_config.hpp +++ b/src/smith/numerics/solver_config.hpp @@ -103,12 +103,13 @@ struct TimesteppingOptions { /// Linear solution method indicator enum class LinearSolver { - CG, /**< Conjugate gradient */ - GMRES, /**< Generalized minimal residual method */ - SuperLU, /**< SuperLU MPI-enabled direct nodal solver */ - Strumpack, /**< Strumpack MPI-enabled direct frontal solver*/ - PetscCG, /**< PETSc MPI-enabled conjugate gradient solver */ - PetscGMRES /**< PETSc MPI-enabled generalize minimal residual solver */ + CG, /**< Conjugate gradient */ + GMRES, /**< Generalized minimal residual method */ + SuperLU, /**< SuperLU MPI-enabled direct nodal solver */ + Strumpack, /**< Strumpack MPI-enabled direct frontal solver*/ + PetscCG, /**< PETSc MPI-enabled conjugate gradient solver */ + PetscGMRES, /**< PETSc MPI-enabled generalize minimal residual solver */ + None /**< Preconditioner application only, No linear solver Krylov iterations */ }; // _linear_solvers_end @@ -128,6 +129,8 @@ inline std::string linearName(const LinearSolver& s) return "PetscCG"; case LinearSolver::PetscGMRES: return "PetscGMRES"; + case LinearSolver::None: + return "None"; } // This cannot happen, but GCC doesn't know that return "UNKNOWN"; @@ -141,6 +144,7 @@ inline std::map linearSolverMap = { {"CG", LinearSolver::CG}, {"GMRES", LinearSolver::GMRES}, {"SuperLU", LinearSolver::SuperLU}, {"Strumpack", LinearSolver::Strumpack}, {"PetscCG", LinearSolver::PetscCG}, {"PetscGMRES", LinearSolver::PetscGMRES}, + {"None", LinearSolver::None}, }; // Add a custom list of strings? conduit node? diff --git a/src/smith/numerics/solver_with_preconditioner.hpp b/src/smith/numerics/solver_with_preconditioner.hpp new file mode 100644 index 000000000..938c514e1 --- /dev/null +++ b/src/smith/numerics/solver_with_preconditioner.hpp @@ -0,0 +1,55 @@ +// Copyright (c) Lawrence Livermore National Security, LLC and +// other Smith Project Developers. See the top-level LICENSE file for +// details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +#pragma once + +#include + +#include "mfem.hpp" + +namespace smith { + +/// @brief Simple wrapper that owns a linear solver and its preconditioner. +/// +/// This is used to keep a preconditioner alive when it is referenced by an +/// iterative solver (e.g. GMRES) via SetPreconditioner(). +class SolverWithPreconditioner : public mfem::Solver { + public: + /// @brief Construct from an owned linear solver and (optional) preconditioner. + /// @param linear_solver Owned linear solver (must be non-null). + /// @param preconditioner Owned preconditioner (may be null). + SolverWithPreconditioner(std::unique_ptr linear_solver, std::unique_ptr preconditioner) + : linear_solver_(std::move(linear_solver)), preconditioner_(std::move(preconditioner)) + { + MFEM_VERIFY(linear_solver_ != nullptr, "SolverWithPreconditioner requires a non-null linear solver"); + } + + /// @brief Set the operator on the underlying linear solver. + /// @param op Operator to be solved/applied. + void SetOperator(const mfem::Operator& op) override + { + height = op.Height(); + width = op.Width(); + linear_solver_->SetOperator(op); + } + + /// @brief Apply the underlying linear solver. + /// @param x Input vector. + /// @param y Output vector. + void Mult(const mfem::Vector& x, mfem::Vector& y) const override { linear_solver_->Mult(x, y); } + + /// @brief Non-owning access to the underlying linear solver. + mfem::Solver* linearSolver() const { return linear_solver_.get(); } + + /// @brief Non-owning access to the owned preconditioner (may be null). + mfem::Solver* preconditioner() const { return preconditioner_.get(); } + + private: + std::unique_ptr linear_solver_; + std::unique_ptr preconditioner_; +}; + +} // namespace smith diff --git a/src/smith/numerics/tests/CMakeLists.txt b/src/smith/numerics/tests/CMakeLists.txt index 10e693b21..f5d879b55 100644 --- a/src/smith/numerics/tests/CMakeLists.txt +++ b/src/smith/numerics/tests/CMakeLists.txt @@ -13,6 +13,7 @@ set(numerics_serial_test_sources test_block_preconditioner.cpp test_block_preconditioner_backend.cpp test_block_preconditioner_custom_operators.cpp + test_linear_solver_none.cpp ) smith_add_tests( SOURCES ${numerics_serial_test_sources} diff --git a/src/smith/numerics/tests/test_linear_solver_none.cpp b/src/smith/numerics/tests/test_linear_solver_none.cpp new file mode 100644 index 000000000..6d4994b52 --- /dev/null +++ b/src/smith/numerics/tests/test_linear_solver_none.cpp @@ -0,0 +1,111 @@ +// Copyright (c) Lawrence Livermore National Security, LLC and +// other Smith Project Developers. See the top-level LICENSE file for +// details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +#include "gtest/gtest.h" +#include "mfem.hpp" +#include "smith/numerics/equation_solver.hpp" +#include "smith/numerics/solver_config.hpp" +#include "smith/infrastructure/application_manager.hpp" + +namespace smith { + +/** + * @brief Simple identity-like operator: A*x = x + */ +class IdentityOperator : public mfem::Operator { + public: + IdentityOperator(int size) : mfem::Operator(size) {} + void Mult(const mfem::Vector& x, mfem::Vector& y) const override { y = x; } + mfem::Operator& GetGradient(const mfem::Vector& /*x*/) const override { return const_cast(*this); } +}; + +/** + * @brief Simple diagonal operator: A*x = d * x + */ +class DiagonalOperator : public mfem::Operator { + public: + DiagonalOperator(const mfem::Vector& d) : mfem::Operator(d.Size()), d_(d) {} + void Mult(const mfem::Vector& x, mfem::Vector& y) const override + { + for (int i = 0; i < height; i++) { + y(i) = d_(i) * x(i); + } + } + mfem::Operator& GetGradient(const mfem::Vector& /*x*/) const override { return const_cast(*this); } + + private: + const mfem::Vector& d_; +}; + +TEST(LinearSolverNone, Identity) +{ + int size = 10; + IdentityOperator op(size); + + LinearSolverOptions linear_opts; + linear_opts.linear_solver = LinearSolver::None; + linear_opts.preconditioner = Preconditioner::None; + + NonlinearSolverOptions nonlinear_opts; + nonlinear_opts.nonlin_solver = NonlinearSolver::Newton; + nonlinear_opts.max_iterations = 1; + + EquationSolver solver(nonlinear_opts, linear_opts, MPI_COMM_WORLD); + solver.setOperator(op); + + mfem::Vector x(size); + x = 1.0; // Initial guess + // Residual will be f(x) = x. + // x_new = x - [df/dx]^-1 * f(x) = x - I^-1 * x = 0. + + solver.solve(x); + + for (int i = 0; i < size; i++) { + EXPECT_NEAR(x(i), 0.0, 1e-12); + } +} + +TEST(LinearSolverNone, Jacobi) +{ + int size = 10; + mfem::Vector d(size); + d = 2.0; + DiagonalOperator op(d); + + LinearSolverOptions linear_opts; + linear_opts.linear_solver = LinearSolver::None; + linear_opts.preconditioner = Preconditioner::None; // We'll set this manually if needed, but wait. + // Actually, Preconditioner::None with LinearSolver::None should be Identity. + // Let's test that first. + + NonlinearSolverOptions nonlinear_opts; + nonlinear_opts.nonlin_solver = NonlinearSolver::Newton; + nonlinear_opts.max_iterations = 1; + + EquationSolver solver(nonlinear_opts, linear_opts, MPI_COMM_WORLD); + solver.setOperator(op); + + mfem::Vector x(size); + x = 1.0; + // f(x) = 2x. df/dx = 2I. + // If LinearSolver::None and Preconditioner::None, it uses identity: + // x_new = x - I * (2x) = x - 2x = -x. + + solver.solve(x); + + for (int i = 0; i < size; i++) { + EXPECT_NEAR(x(i), -1.0, 1e-12); + } +} + +} // namespace smith + +int main(int argc, char* argv[]) +{ + ::testing::InitGoogleTest(&argc, argv); + smith::ApplicationManager applicationManager(argc, argv); + return RUN_ALL_TESTS(); +}