nonlinear_block_solver.hpp

Go to the documentation of this file.

// 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 <memory>
#include <functional>
#include <optional>
#include <vector>
#include <mpi.h>
 
#include "smith/numerics/solver_config.hpp"
#include "smith/numerics/nonlinear_convergence.hpp"
 
namespace mfem {
class Solver;
class Vector;
class HypreParMatrix;
class BlockOperator;
}  // namespace mfem
 
namespace smith {
 
class EquationSolver;
class BoundaryConditionManager;
class FiniteElementState;
class FiniteElementDual;
class Mesh;
struct NonlinearSolverOptions;
struct LinearSolverOptions;
 
class NonlinearBlockSolverBase {
 public:
  virtual ~NonlinearBlockSolverBase() {}
 
  using FieldT = FiniteElementState;                        
  using FieldPtr = std::shared_ptr<FieldT>;                 
  using FieldD = FiniteElementDual;                         
  using DualPtr = std::shared_ptr<FieldD>;                  
  using MatrixPtr = std::unique_ptr<mfem::HypreParMatrix>;  
 
  virtual std::vector<FieldPtr> solve(
      const std::vector<FieldPtr>& u_guesses,
      std::function<std::vector<mfem::Vector>(const std::vector<FieldPtr>&)> residuals,
      std::function<std::vector<std::vector<MatrixPtr>>(const std::vector<FieldPtr>&)> jacobians) const = 0;
 
  virtual std::vector<FieldPtr> solveAdjoint(const std::vector<DualPtr>& u_bars,
                                             std::vector<std::vector<MatrixPtr>>& jacobian_transposed) const = 0;
 
  ConvergenceStatus convergenceStatus(double tolerance_multiplier, const std::vector<mfem::Vector>& residuals) const
  {
    return convergenceStatus(tolerance_multiplier, residuals, convergence_context_);
  }
 
  virtual ConvergenceStatus convergenceStatus(double tolerance_multiplier, const std::vector<mfem::Vector>& residuals,
                                              NonlinearConvergenceContext& context) const = 0;
 
  virtual void primeConvergenceContext(const std::vector<mfem::Vector>& residuals,
                                       NonlinearConvergenceContext& context) const = 0;
 
  bool checkConvergence(double tolerance_multiplier, const std::vector<mfem::Vector>& residuals) const
  {
    return convergenceStatus(tolerance_multiplier, residuals).converged;
  }
 
  virtual void resetConvergenceState() const { convergence_context_.reset(); }
 
  virtual void clearMemory() const {}
 
  virtual void setInnerToleranceMultiplier(double multiplier) = 0;
 
 protected:
  mutable bool is_setup_ = false;  
  mutable NonlinearConvergenceContext convergence_context_ = {};  
};
 
class NonlinearBlockSolver : public NonlinearBlockSolverBase {
 public:
  using NonlinearBlockSolverBase::convergenceStatus;
 
  NonlinearBlockSolver(std::unique_ptr<EquationSolver> s, MPI_Comm comm, double abs_tol = 1e-12, double rel_tol = 1e-8,
                       std::optional<NonlinearSolverOptions> retained_nonlinear_options = std::nullopt,
                       std::optional<LinearSolverOptions> retained_linear_options = std::nullopt);
 
  ConvergenceStatus convergenceStatus(double tolerance_multiplier, const std::vector<mfem::Vector>& residuals,
                                      NonlinearConvergenceContext& context) const override;
 
  void primeConvergenceContext(const std::vector<mfem::Vector>& residuals,
                               NonlinearConvergenceContext& context) const override;
 
  std::vector<FieldPtr> solve(
      const std::vector<FieldPtr>& u_guesses,
      std::function<std::vector<mfem::Vector>(const std::vector<FieldPtr>&)> residuals,
      std::function<std::vector<std::vector<MatrixPtr>>(const std::vector<FieldPtr>&)> jacobians) const override;
 
  std::vector<FieldPtr> solveAdjoint(const std::vector<DualPtr>& u_bars,
                                     std::vector<std::vector<MatrixPtr>>& jacobian_transposed) const override;
 
  void completeSetup(const std::vector<FieldPtr>& us) const;
 
  void setInnerToleranceMultiplier(double multiplier) override { inner_tol_multiplier_ = multiplier; }
 
  std::shared_ptr<NonlinearBlockSolver> cloneFresh() const;
 
  mutable std::unique_ptr<mfem::BlockOperator>
      block_jac_;  
  mutable std::vector<std::vector<MatrixPtr>>
      matrix_of_jacs_;  
 
  mutable std::unique_ptr<EquationSolver>
      nonlinear_solver_;  
 
  MPI_Comm comm_;                      
  double abs_tol_;                     
  double rel_tol_;                     
  double inner_tol_multiplier_ = 1.0;  
  std::optional<NonlinearSolverOptions> retained_nonlinear_options_ = std::nullopt;  
  std::optional<LinearSolverOptions> retained_linear_options_ = std::nullopt;        
};
 
std::shared_ptr<NonlinearBlockSolver> buildNonlinearBlockSolver(NonlinearSolverOptions nonlinear_opts,
                                                                LinearSolverOptions linear_opts,
                                                                const smith::Mesh& mesh);
 
}  // namespace smith