www.mooseframework.org
PhysicsBasedPreconditioner.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
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 
11 
12 // MOOSE includes
14 #include "FEProblem.h"
15 #include "MooseEnum.h"
16 #include "MooseVariableFE.h"
17 #include "NonlinearSystem.h"
18 #include "PetscSupport.h"
19 
20 #include "libmesh/libmesh_common.h"
21 #include "libmesh/equation_systems.h"
22 #include "libmesh/nonlinear_implicit_system.h"
23 #include "libmesh/nonlinear_solver.h"
24 #include "libmesh/linear_implicit_system.h"
25 #include "libmesh/transient_system.h"
26 #include "libmesh/numeric_vector.h"
27 #include "libmesh/sparse_matrix.h"
28 #include "libmesh/string_to_enum.h"
29 #include "libmesh/coupling_matrix.h"
30 
32 
35 {
37 
38  params.addClassDescription("Physics-based preconditioner (PBP) allows individual physics to have "
39  "their own preconditioner.");
40 
41  params.addRequiredParam<std::vector<std::string>>(
42  "solve_order",
43  "The order the block rows will be solved in. Put the name of variables here "
44  "to stand for solving that variable's block row. A variable may appear more "
45  "than once (to create cylces if you like).");
46  params.addRequiredParam<std::vector<std::string>>("preconditioner", "TODO: docstring");
47 
48  return params;
49 }
50 
52  : MoosePreconditioner(params),
53  Preconditioner<Number>(MoosePreconditioner::_communicator),
54  _nl(_fe_problem.getNonlinearSystemBase(_nl_sys_num))
55 {
56  const auto & libmesh_system = _nl.system();
57  unsigned int num_systems = _nl.system().n_vars();
58  _systems.resize(num_systems);
59  _preconditioners.resize(num_systems);
60  _off_diag.resize(num_systems);
61  _off_diag_mats.resize(num_systems);
62  _pre_type.resize(num_systems);
63 
64  { // Setup the Coupling Matrix so MOOSE knows what we're doing
65  unsigned int n_vars = _nl.nVariables();
66 
67  // The coupling matrix is held and released by FEProblemBase, so it is not released in this
68  // object
69  std::unique_ptr<CouplingMatrix> cm = std::make_unique<CouplingMatrix>(n_vars);
70 
71  bool full = getParam<bool>("full");
72 
73  if (!full)
74  {
75  // put 1s on diagonal
76  for (unsigned int i = 0; i < n_vars; i++)
77  (*cm)(i, i) = 1;
78 
79  // off-diagonal entries
80  for (const auto & off_diag : getParam<NonlinearVariableName, NonlinearVariableName>(
81  "off_diag_row", "off_diag_column"))
82  {
83  const auto row = libmesh_system.variable_number(off_diag.first);
84  const auto column = libmesh_system.variable_number(off_diag.second);
85  (*cm)(row, column) = 1;
86  }
87 
88  // TODO: handle coupling entries between NL-vars and SCALAR-vars
89  }
90  else
91  {
92  for (unsigned int i = 0; i < n_vars; i++)
93  for (unsigned int j = 0; j < n_vars; j++)
94  (*cm)(i, j) = 1;
95  }
96 
97  setCouplingMatrix(std::move(cm));
98  }
99 
100  // PC types
101  const std::vector<std::string> & pc_types = getParam<std::vector<std::string>>("preconditioner");
102  for (unsigned int i = 0; i < num_systems; i++)
103  _pre_type[i] = Utility::string_to_enum<PreconditionerType>(pc_types[i]);
104 
105  // solve order
106  const std::vector<std::string> & solve_order = getParam<std::vector<std::string>>("solve_order");
107  _solve_order.resize(solve_order.size());
108  for (const auto i : index_range(solve_order))
109  _solve_order[i] = _nl.system().variable_number(solve_order[i]);
110 
111  // diag and off-diag systems
112  unsigned int n_vars = _nl.system().n_vars();
113 
114  // off-diagonal entries
115  std::vector<std::vector<unsigned int>> off_diag(n_vars);
116  if (isParamValid("off_diag_row") && isParamValid("off_diag_column"))
117  {
118  const std::vector<NonlinearVariableName> & odr =
119  getParam<std::vector<NonlinearVariableName>>("off_diag_row");
120  const std::vector<NonlinearVariableName> & odc =
121  getParam<std::vector<NonlinearVariableName>>("off_diag_column");
122 
123  for (const auto i : index_range(odr))
124  {
125  unsigned int row = _nl.system().variable_number(odr[i]);
126  unsigned int column = _nl.system().variable_number(odc[i]);
127  off_diag[row].push_back(column);
128  }
129  }
130  // Add all of the preconditioning systems
131  for (unsigned int var = 0; var < n_vars; var++)
132  addSystem(var, off_diag[var], _pre_type[var]);
133 
135 
137  mooseError("PBP must be used with JFNK solve type");
138 }
139 
141 
142 void
144  std::vector<unsigned int> off_diag,
145  PreconditionerType type)
146 {
147  std::string var_name = _nl.system().variable_name(var);
148 
149  LinearImplicitSystem & precond_system =
150  _fe_problem.es().add_system<LinearImplicitSystem>(var_name + "_system");
151  precond_system.assemble_before_solve = false;
152 
153  const std::set<SubdomainID> * active_subdomains = _nl.getVariableBlocks(var);
154  precond_system.add_variable(
155  var_name + "_prec", _nl.system().variable(var).type(), active_subdomains);
156 
157  _systems[var] = &precond_system;
158  _pre_type[var] = type;
159 
160  _off_diag_mats[var].resize(off_diag.size());
161  for (const auto i : index_range(off_diag))
162  {
163  // Add the matrix to hold the off-diagonal piece
164  _off_diag_mats[var][i] = &precond_system.add_matrix(_nl.system().variable_name(off_diag[i]));
165  }
166 
167  _off_diag[var] = off_diag;
168 }
169 
170 void
172 {
173  TIME_SECTION("init", 2, "Initializing PhysicsBasedPreconditioner");
174 
175  // Tell libMesh that this is initialized!
176  _is_initialized = true;
177 
178  const unsigned int num_systems = _systems.size();
179 
180  // If no order was specified, just solve them in increasing order
181  if (_solve_order.size() == 0)
182  {
183  _solve_order.resize(num_systems);
184  for (unsigned int i = 0; i < num_systems; i++)
185  _solve_order[i] = i;
186  }
187 
188  // Loop over variables
189  for (unsigned int system_var = 0; system_var < num_systems; system_var++)
190  {
191  LinearImplicitSystem & u_system = *_systems[system_var];
192 
193  if (!_preconditioners[system_var])
194  _preconditioners[system_var] =
196 
197  // we have to explicitly set the matrix in the preconditioner, because h-adaptivity could have
198  // changed it and we have to work with the current one
199  Preconditioner<Number> * preconditioner = _preconditioners[system_var].get();
200  preconditioner->set_matrix(u_system.get_system_matrix());
201  preconditioner->set_type(_pre_type[system_var]);
202 
203  preconditioner->init();
204  }
205 }
206 
207 void
209 {
210  const unsigned int num_systems = _systems.size();
211 
212  std::vector<JacobianBlock *> blocks;
213 
214  // Loop over variables
215  for (unsigned int system_var = 0; system_var < num_systems; system_var++)
216  {
217  LinearImplicitSystem & u_system = *_systems[system_var];
218 
219  {
220  JacobianBlock * block =
221  new JacobianBlock(u_system, u_system.get_system_matrix(), system_var, system_var);
222  blocks.push_back(block);
223  }
224 
225  for (const auto diag : index_range(_off_diag[system_var]))
226  {
227  unsigned int coupled_var = _off_diag[system_var][diag];
228  std::string coupled_name = _nl.system().variable_name(coupled_var);
229 
230  JacobianBlock * block =
231  new JacobianBlock(u_system, *_off_diag_mats[system_var][diag], system_var, coupled_var);
232  blocks.push_back(block);
233  }
234  }
235 
237 
238  // cleanup
239  for (auto & block : blocks)
240  delete block;
241 }
242 
243 void
245 {
246  TIME_SECTION("apply", 1, "Applying PhysicsBasedPreconditioner");
247 
248  const unsigned int num_systems = _systems.size();
249 
251 
252  // Zero out the solution vectors
253  for (unsigned int sys = 0; sys < num_systems; sys++)
254  _systems[sys]->solution->zero();
255 
256  // Loop over solve order
257  for (const auto i : index_range(_solve_order))
258  {
259  unsigned int system_var = _solve_order[i];
260 
261  LinearImplicitSystem & u_system = *_systems[system_var];
262 
263  // Copy rhs from the big system into the small one
265  mesh, _nl.system().number(), system_var, x, u_system.number(), 0, *u_system.rhs);
266 
267  // Modify the RHS by subtracting off the matvecs of the solutions for the other preconditioning
268  // systems with the off diagonal blocks in this system.
269  for (const auto diag : index_range(_off_diag[system_var]))
270  {
271  unsigned int coupled_var = _off_diag[system_var][diag];
272  LinearImplicitSystem & coupled_system = *_systems[coupled_var];
273  SparseMatrix<Number> & off_diag = *_off_diag_mats[system_var][diag];
274  NumericVector<Number> & rhs = *u_system.rhs;
275 
276  // This next bit computes rhs -= A*coupled_solution
277  // It does what it does because there is no vector_mult_sub()
278  rhs.close();
279  rhs.scale(-1.0);
280  rhs.close();
281  off_diag.vector_mult_add(rhs, *coupled_system.solution);
282  rhs.close();
283  rhs.scale(-1.0);
284  rhs.close();
285  }
286 
287  // If there is no off_diag, then u_system.rhs will not be closed.
288  // Thus, we need to close it right here
289  if (!_off_diag[system_var].size())
290  u_system.rhs->close();
291 
292  // Apply the preconditioner to the small system
293  _preconditioners[system_var]->apply(*u_system.rhs, *u_system.solution);
294 
295  // Copy solution from small system into the big one
296  // copyVarValues(mesh,system,0,*u_system.solution,0,system_var,y);
297  }
298 
299  // Copy the solutions out
300  for (unsigned int system_var = 0; system_var < num_systems; system_var++)
301  {
302  LinearImplicitSystem & u_system = *_systems[system_var];
303 
305  mesh, u_system.number(), 0, *u_system.solution, _nl.system().number(), system_var, y);
306  }
307 
308  y.close();
309 }
310 
311 void
313 {
314 }
virtual void setup()
This is called every time the "operator might have changed".
virtual const std::set< SubdomainID > * getVariableBlocks(unsigned int var_number)
Get the block where a variable of this system is defined.
Definition: SystemBase.C:153
Helper class for holding the preconditioning blocks to fill.
virtual void apply(const NumericVector< Number > &x, NumericVector< Number > &y)
Computes the preconditioned vector "y" based on input "x".
static void copyVarValues(MeshBase &mesh, const unsigned int from_system, const unsigned int from_var, const NumericVector< Number > &from_vector, const unsigned int to_system, const unsigned int to_var, NumericVector< Number > &to_vector)
Helper function for copying values associated with variables in vectors from two different systems...
std::vector< std::vector< unsigned int > > _off_diag
Holds which off diagonal blocks to compute.
SolverParams & solverParams()
Get the solver parameters.
virtual void attachPreconditioner(Preconditioner< Number > *preconditioner)=0
Attach a customized preconditioner that requires physics knowledge.
PreconditionerType type() const
MeshBase & mesh
NonlinearSystemBase & _nl
The nonlinear system this PBP is associated with (convenience reference)
The main MOOSE class responsible for handling user-defined parameters in almost every MOOSE system...
std::vector< PreconditionerType > _pre_type
Which preconditioner to use for each solve.
std::vector< LinearImplicitSystem * > _systems
List of linear system that build up the preconditioner.
void set_matrix(SparseMatrix< Number > &mat)
FEProblemBase & _fe_problem
Subproblem this preconditioner is part of.
const Parallel::Communicator & _communicator
virtual void clear()
Release all memory and clear data structures.
void addRequiredParam(const std::string &name, const std::string &doc_string)
This method adds a parameter and documentation string to the InputParameters object that will be extr...
Base class for MOOSE preconditioners.
virtual EquationSystems & es() override
virtual void init()
Initialize data structures if not done so already.
bool isParamValid(const std::string &name) const
Test if the supplied parameter is valid.
virtual unsigned int nVariables() const
Get the number of variables in this system.
Definition: SystemBase.C:840
virtual void computeJacobianBlocks(std::vector< JacobianBlock *> &blocks, const unsigned int nl_sys_num)
Computes several Jacobian blocks simultaneously, summing their contributions into smaller preconditio...
virtual void scale(const Number factor)=0
static InputParameters validParams()
Constructor.
unsigned int n_vars
PhysicsBasedPreconditioner(const InputParameters &params)
Jacobian-Free Newton Krylov.
Definition: MooseTypes.h:758
MooseMesh wraps a libMesh::Mesh object and enhances its capabilities by caching additional data and s...
Definition: MooseMesh.h:88
Implements a segregated solve preconditioner.
std::vector< std::vector< SparseMatrix< Number > * > > _off_diag_mats
Holds pointers to the off-diagonal matrices.
Moose::SolveType _type
Definition: SolverParams.h:19
void addSystem(unsigned int var, std::vector< unsigned int > off_diag, PreconditionerType type=AMG_PRECOND)
Add a diagonal system + possibly off-diagonals ones as well, also specifying type of preconditioning...
static InputParameters validParams()
unsigned int number() const
Gets the number of this system.
Definition: SystemBase.C:1125
virtual void close()=0
PreconditionerType
void vector_mult_add(NumericVector< Number > &dest, const NumericVector< Number > &arg) const
virtual System & system() override
Get the reference to the libMesh system.
std::vector< std::unique_ptr< Preconditioner< Number > > > _preconditioners
Holds one Preconditioner object per small system to solve.
virtual MooseMesh & mesh() override
std::vector< unsigned int > _solve_order
Holds the order the blocks are solved for.
void mooseError(Args &&... args) const
Emits an error prefixed with object name and type.
void addClassDescription(const std::string &doc_string)
This method adds a description of the class that will be displayed in the input file syntax dump...
void set_type(const PreconditionerType pct)
static std::unique_ptr< Preconditioner< Number > > build_preconditioner(const libMesh::Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package())
registerMooseObjectAliased("MooseApp", PhysicsBasedPreconditioner, "PBP")
Real Number
void setCouplingMatrix(std::unique_ptr< CouplingMatrix > cm)
Setup the coupling matrix on the finite element problem.
auto index_range(const T &sizable)