Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "libmesh/petsc_macro.h"
      11             : #include "FieldSplitPreconditioner.h"
      12             : 
      13             : // MOOSE includes
      14             : #include "FEProblem.h"
      15             : #include "MooseEnum.h"
      16             : #include "MooseVariableFE.h"
      17             : #include "NonlinearSystem.h"
      18             : #include "PetscSupport.h"
      19             : #include "MoosePreconditioner.h"
      20             : #include "MooseStaticCondensationPreconditioner.h"
      21             : #include "Split.h"
      22             : #include "PetscDMMoose.h"
      23             : 
      24             : #include "libmesh/libmesh_common.h"
      25             : #include "libmesh/petsc_nonlinear_solver.h"
      26             : #include "libmesh/coupling_matrix.h"
      27             : 
      28             : using namespace libMesh;
      29             : 
      30             : template <typename Base>
      31             : InputParameters
      32       29768 : FieldSplitPreconditionerTempl<Base>::validParams()
      33             : {
      34       29768 :   InputParameters params = Base::validParams();
      35       59536 :   params.addClassDescription("Preconditioner designed to map onto PETSc's PCFieldSplit.");
      36             : 
      37      119072 :   params.addRequiredParam<std::string>(
      38             :       "topsplit", "Entrance to splits, the top split will specify how splits will go.");
      39             :   // We should use full coupling Jacobian matrix by default
      40       59536 :   params.addParam<bool>("full",
      41       59536 :                         true,
      42             :                         "Set to true if you want the full set of couplings between variables "
      43             :                         "simply for convenience so you don't have to set every off_diag_row "
      44             :                         "and off_diag_column combination.");
      45       29768 :   return params;
      46           0 : }
      47             : 
      48             : template <typename Base>
      49         139 : FieldSplitPreconditionerTempl<Base>::FieldSplitPreconditionerTempl(
      50             :     const InputParameters & parameters)
      51             :   : Base(parameters),
      52         278 :     _nl(this->_fe_problem.getNonlinearSystemBase(this->_nl_sys_num)),
      53         278 :     _decomposition_split(this->template getParam<std::string>("topsplit"))
      54             : {
      55         139 :   if (libMesh::default_solver_package() != libMesh::PETSC_SOLVERS)
      56           0 :     mooseError("The field split preconditioner can only be used with PETSc");
      57             : 
      58             :   // number of variables
      59         139 :   unsigned int n_vars = _nl.nVariables();
      60             :   // if we want to construct a full Jacobian?
      61             :   // it is recommended to have a full Jacobian for using
      62             :   // the fieldSplit preconditioner
      63         278 :   bool full = this->template getParam<bool>("full");
      64             : 
      65             :   // how variables couple
      66         139 :   std::unique_ptr<CouplingMatrix> cm = std::make_unique<CouplingMatrix>(n_vars);
      67         139 :   if (!full)
      68             :   {
      69           0 :     if (this->isParamValid("off_diag_row") && this->isParamValid("off_diag_column"))
      70             :     {
      71             : 
      72           0 :       const auto off_diag_rows =
      73             :           this->template getParam<std::vector<NonlinearVariableName>>("off_diag_row");
      74           0 :       const auto off_diag_columns =
      75             :           this->template getParam<std::vector<NonlinearVariableName>>("off_diag_column");
      76             : 
      77             :       // put 1s on diagonal
      78           0 :       for (unsigned int i = 0; i < n_vars; i++)
      79           0 :         (*cm)(i, i) = 1;
      80             : 
      81             :       // off-diagonal entries
      82           0 :       std::vector<std::vector<unsigned int>> off_diag(n_vars);
      83           0 :       if (off_diag_rows.size() * off_diag_columns.size() != 0 &&
      84           0 :           off_diag_rows.size() == off_diag_columns.size())
      85           0 :         for (const auto i : index_range(off_diag_rows))
      86             :         {
      87           0 :           unsigned int row = _nl.getVariable(0, off_diag_rows[i]).number();
      88           0 :           unsigned int column = _nl.getVariable(0, off_diag_columns[i]).number();
      89           0 :           (*cm)(row, column) = 1;
      90             :         }
      91           0 :     }
      92             :   }
      93             :   else
      94             :   {
      95         434 :     for (unsigned int i = 0; i < n_vars; i++)
      96         936 :       for (unsigned int j = 0; j < n_vars; j++)
      97         641 :         (*cm)(i, j) = 1; // full coupling
      98             :   }
      99         139 :   this->setCouplingMatrix(std::move(cm));
     100             : 
     101             :   // turn on a flag
     102         139 :   _nl.useFieldSplitPreconditioner(this);
     103         139 : }
     104             : 
     105             : template <typename Base>
     106             : void
     107         129 : FieldSplitPreconditionerTempl<Base>::createMooseDM(DM * dm)
     108             : {
     109         129 :   LibmeshPetscCallA(
     110             :       _nl.comm().get(),
     111             :       DMCreateMoose(_nl.comm().get(), _nl, dofMap(), system(), _decomposition_split, dm));
     112         129 :   LibmeshPetscCallA(_nl.comm().get(),
     113             :                     PetscObjectSetOptionsPrefix((PetscObject)*dm, prefix().c_str()));
     114         129 :   LibmeshPetscCall(DMSetFromOptions(*dm));
     115         129 :   LibmeshPetscCall(DMSetUp(*dm));
     116         129 : }
     117             : 
     118             : registerMooseObjectAliased("MooseApp", FieldSplitPreconditioner, "FSP");
     119             : 
     120             : InputParameters
     121       14997 : FieldSplitPreconditioner::validParams()
     122             : {
     123       14997 :   return FieldSplitPreconditionerTempl<MoosePreconditioner>::validParams();
     124             : }
     125             : 
     126         126 : FieldSplitPreconditioner::FieldSplitPreconditioner(const InputParameters & params)
     127         126 :   : FieldSplitPreconditionerTempl<MoosePreconditioner>(params)
     128             : {
     129         126 :   std::shared_ptr<Split> top_split = _nl.getSplit(_decomposition_split);
     130         126 :   top_split->setup(_nl, _nl.prefix());
     131         126 : }
     132             : 
     133             : void
     134         117 : FieldSplitPreconditioner::setupDM()
     135             : {
     136             :   PetscBool ismoose;
     137         117 :   DM dm = LIBMESH_PETSC_NULLPTR;
     138             : 
     139             :   // Initialize the part of the DM package that's packaged with Moose; in the PETSc source tree this
     140             :   // call would be in DMInitializePackage()
     141         117 :   LibmeshPetscCall(DMMooseRegisterAll());
     142             :   // Create and set up the DM that will consume the split options and deal with block matrices.
     143             :   auto * const petsc_solver =
     144         117 :       libMesh::cast_ptr<PetscNonlinearSolver<Number> *>(_nl.nonlinearSolver());
     145         117 :   SNES snes = petsc_solver->snes(prefix().c_str());
     146             :   // if there exists a DMMoose object, not to recreate a new one
     147         117 :   LibmeshPetscCall(SNESGetDM(snes, &dm));
     148         117 :   if (dm)
     149             :   {
     150         117 :     LibmeshPetscCall(PetscObjectTypeCompare((PetscObject)dm, DMMOOSE, &ismoose));
     151         117 :     if (ismoose)
     152           0 :       return;
     153             :   }
     154         117 :   createMooseDM(&dm);
     155         117 :   LibmeshPetscCall(SNESSetDM(snes, dm));
     156         117 :   LibmeshPetscCall(DMDestroy(&dm));
     157             : }
     158             : 
     159             : const DofMapBase &
     160         117 : FieldSplitPreconditioner::dofMap() const
     161             : {
     162         117 :   return _nl.dofMap();
     163             : }
     164             : 
     165             : const System &
     166         117 : FieldSplitPreconditioner::system() const
     167             : {
     168         117 :   return _nl.system();
     169             : }
     170             : 
     171             : std::string
     172         234 : FieldSplitPreconditioner::prefix() const
     173             : {
     174         234 :   return _nl.prefix();
     175             : }
     176             : 
     177             : KSP
     178           0 : FieldSplitPreconditioner::getKSP()
     179             : {
     180             :   KSP ksp;
     181           0 :   auto snes = _nl.getSNES();
     182           0 :   LibmeshPetscCall(SNESGetKSP(snes, &ksp));
     183           0 :   return ksp;
     184             : }
     185             : 
     186             : template class FieldSplitPreconditionerTempl<MooseStaticCondensationPreconditioner>;