linear_solver.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan
// SPDX-License-Identifier: GNU Lesser General Public Version 2.1
 
#pragma once
 
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/solver_control.h>
 
#include <prismspf/core/conditional_ostreams.h>
#include <prismspf/core/group_solution_handler.h>
#include <prismspf/core/invm_manager.h>
#include <prismspf/core/timer.h>
#include <prismspf/core/type_enums.h>
#include <prismspf/core/types.h>
 
#include <prismspf/solvers/mf_operator.h>
#include <prismspf/solvers/solver_base.h>
 
#include <prismspf/user_inputs/linear_solve_parameters.h>
 
#include <prismspf/config.h>
 
PRISMS_PF_BEGIN_NAMESPACE
 
template <unsigned int dim, unsigned int degree, typename number>
class SolveContext;

template <unsigned int dim, unsigned int degree, typename number>
class LinearSolver : public SolverBase<dim, degree, number>
{
protected:
  using SolverBase<dim, degree, number>::solutions;
  using SolverBase<dim, degree, number>::solve_context;
  using SolverBase<dim, degree, number>::solve_block;
 
public:
  LinearSolver(SolveBlock                               _solve_block,
               const SolveContext<dim, degree, number> &_solve_context)
    : SolverBase<dim, degree, number>(_solve_block, _solve_context)
    , lin_params(
        solve_context->get_user_inputs().linear_solve_parameters.linear_solvers.at(
          solve_block.id))
  {}

  void
  init(const std::list<DependencyMap> &all_dependeny_sets) override
  {
    SolverBase<dim, degree, number>::init(all_dependeny_sets);
    unsigned int num_levels = solve_context->get_dof_manager().get_dof_handlers().size();
    rhs_vector.resize(num_levels);
    for (unsigned int relative_level = 0; relative_level < num_levels; ++relative_level)
      {
        rhs_vector[relative_level].reinit(
          solutions.get_solution_full_vector(relative_level));
      }
    // Initialize rhs_operators
    rhs_operators.reserve(num_levels);
    for (unsigned int relative_level = 0; relative_level < num_levels; ++relative_level)
      {
        rhs_operators.emplace_back(solve_context->get_pde_operator(),
                                   &PDEOperatorBase<dim, degree, number>::compute_rhs,
                                   solve_context->get_field_attributes(),
                                   solve_context->get_solution_indexer(),
                                   relative_level,
                                   solve_block.dependencies_rhs,
                                   solve_context->get_simulation_timer());
        rhs_operators[relative_level].initialize(solutions);
        rhs_operators[relative_level].set_scaling_diagonal(
          lin_params.tolerance_type != AbsoluteResidual,
          solve_context->get_invm_manager().get_invm_sqrt(
            solve_context->get_field_attributes(),
            solve_block.field_indices,
            relative_level));
      }
    // Initialize lhs_operators
    lhs_operators.reserve(num_levels);
    for (unsigned int relative_level = 0; relative_level < num_levels; ++relative_level)
      {
        lhs_operators.emplace_back(solve_context->get_pde_operator(),
                                   &PDEOperatorBase<dim, degree, number>::compute_lhs,
                                   solve_context->get_field_attributes(),
                                   solve_context->get_solution_indexer(),
                                   relative_level,
                                   solve_block.dependencies_lhs,
                                   solve_context->get_simulation_timer());
        lhs_operators[relative_level].initialize(solutions);
        lhs_operators[relative_level].set_scaling_diagonal(
          lin_params.tolerance_type != AbsoluteResidual,
          solve_context->get_invm_manager().get_invm_sqrt(
            solve_context->get_field_attributes(),
            solve_block.field_indices,
            relative_level));
      }
    linear_solver_control.set_max_steps(lin_params.max_iterations);
    linear_solver_control.set_tolerance(lin_params.tolerance * normalization_value());
    inhomogenous_values.reinit(solutions.get_solution_full_vector(0));
    solutions.apply_constraints(inhomogenous_values, 0);
    inhomogenous_rhs.reinit(solutions.get_solution_full_vector(0));
  }

  void
  reinit() override
  {
    SolverBase<dim, degree, number>::reinit();
    const unsigned int num_levels = rhs_vector.size();
    for (unsigned int relative_level = 0; relative_level < num_levels; ++relative_level)
      {
        rhs_vector[relative_level].reinit(
          solutions.get_solution_full_vector(relative_level));
      }
    inhomogenous_values.reinit(solutions.get_solution_full_vector(0));
    solutions.apply_constraints(inhomogenous_values, 0);
    inhomogenous_rhs.reinit(solutions.get_solution_full_vector(0));
  }

  void
  solve_level(unsigned int relative_level) override
  {
    // Zero out the ghosts
    Timer::start_section("Zero ghosts");
    solutions.zero_out_ghosts(relative_level);
    Timer::end_section("Zero ghosts");
 
    // Set up rhs vector
    rhs_operators[relative_level].compute_operator(rhs_vector[relative_level]);
    if (relative_level == 0)
      {
        lhs_operators[0].read_plain = true;
        lhs_operators[0].compute_operator(inhomogenous_rhs, inhomogenous_values);
        lhs_operators[0].read_plain = false;
        rhs_vector[0] -= inhomogenous_rhs;
      }
    // Linear solve
    do_linear_solve(rhs_vector[relative_level],
                    lhs_operators[relative_level],
                    solutions.get_solution_full_vector(relative_level));
 
    if (relative_level == 0)
      {
        solutions.get_solution_full_vector(0) += inhomogenous_values;
      }
    // Apply constraints
    solutions.apply_constraints(relative_level);
 
    // Update the ghosts
    Timer::start_section("Update ghosts");
    solutions.update_ghosts(relative_level);
    Timer::end_section("Update ghosts");
  }
 
  int
  do_linear_solve(BlockVector<number>             &b_vector,
                  MFOperator<dim, degree, number> &lhs_operator,
                  BlockVector<number>             &x_vector)
  {
    // Linear solve
    try
      {
        dealii::SolverCG<BlockVector<number>> cg_solver(linear_solver_control);
        cg_solver.solve(lhs_operator, x_vector, b_vector, dealii::PreconditionIdentity());
        if (solve_context->get_user_inputs().output_parameters.should_output(
              solve_context->get_simulation_timer().get_increment()))
          {
            ConditionalOStreams::pout_summary()
              << " Linear solve final residual : "
              << linear_solver_control.last_value() / normalization_value()
              << " Linear steps: " << linear_solver_control.last_step() << "\n"
              << std::flush;
          }
      }
    catch (...) // TODO: more specific catch
      {
        ConditionalOStreams::pout_base()
          << "[Increment " << solve_context->get_simulation_timer().get_increment()
          << "] "
          << "Warning: linear solver did not converge as per set tolerances before "
          << lin_params.max_iterations << " iterations.\n";
      }
    return linear_solver_control.last_step();
  }
 
protected:
  std::vector<MFOperator<dim, degree, number>> rhs_operators;
 
  std::vector<MFOperator<dim, degree, number>> lhs_operators;
  std::vector<BlockVector<number>>             rhs_vector;
 
  double
  normalization_value()
  {
    SolverToleranceType type  = lin_params.tolerance_type;
    double              value = 1.0;
    if (type == RMSEPerField || type == RMSETotal)
      {
        value *= std::sqrt(solve_context->get_triangulation_manager().get_volume());
      }
    if (type == RMSEPerField || type == IntegratedPerField)
      {
        value *= std::sqrt(double(solve_block.field_indices.size()));
      }
    return value;
  }
 
private:
  LinearSolverParameters lin_params;

  dealii::SolverControl linear_solver_control;

  BlockVector<number> inhomogenous_values;

  BlockVector<number> inhomogenous_rhs;
};
 
PRISMS_PF_END_NAMESPACE