www.mooseframework.org
PhysicsBasedPreconditioner.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* DO NOT MODIFY THIS HEADER */
3 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
4 /* */
5 /* (c) 2010 Battelle Energy Alliance, LLC */
6 /* ALL RIGHTS RESERVED */
7 /* */
8 /* Prepared by Battelle Energy Alliance, LLC */
9 /* Under Contract No. DE-AC07-05ID14517 */
10 /* With the U. S. Department of Energy */
11 /* */
12 /* See COPYRIGHT for full restrictions */
13 /****************************************************************/
14 
16 
17 // MOOSE includes
19 #include "FEProblem.h"
20 #include "MooseEnum.h"
21 #include "MooseVariable.h"
22 #include "NonlinearSystem.h"
23 #include "PetscSupport.h"
24 
25 #include "libmesh/libmesh_common.h"
26 #include "libmesh/equation_systems.h"
27 #include "libmesh/nonlinear_implicit_system.h"
28 #include "libmesh/nonlinear_solver.h"
29 #include "libmesh/linear_implicit_system.h"
30 #include "libmesh/transient_system.h"
31 #include "libmesh/numeric_vector.h"
32 #include "libmesh/sparse_matrix.h"
33 #include "libmesh/string_to_enum.h"
34 #include "libmesh/coupling_matrix.h"
35 
36 template <>
39 {
41 
42  params.addRequiredParam<std::vector<std::string>>(
43  "solve_order",
44  "The order the block rows will be solved in. Put the name of variables here "
45  "to stand for solving that variable's block row. A variable may appear more "
46  "than once (to create cylces if you like).");
47  params.addRequiredParam<std::vector<std::string>>("preconditioner", "TODO: docstring");
48 
49  params.addParam<std::vector<std::string>>(
50  "off_diag_row",
51  "The off diagonal row you want to add into the matrix, it will be associated "
52  "with an off diagonal column from the same position in off_diag_colum.");
53  params.addParam<std::vector<std::string>>("off_diag_column",
54  "The off diagonal column you want to add into the "
55  "matrix, it will be associated with an off diagonal "
56  "row from the same position in off_diag_row.");
57 
58  return params;
59 }
60 
62  : MoosePreconditioner(params),
63  Preconditioner<Number>(MoosePreconditioner::_communicator),
64  _nl(_fe_problem.getNonlinearSystemBase())
65 {
66  unsigned int num_systems = _nl.system().n_vars();
67  _systems.resize(num_systems);
68  _preconditioners.resize(num_systems);
69  _off_diag.resize(num_systems);
70  _off_diag_mats.resize(num_systems);
71  _pre_type.resize(num_systems);
72 
73  { // Setup the Coupling Matrix so MOOSE knows what we're doing
75  unsigned int n_vars = nl.nVariables();
76 
77  // The coupling matrix is held and released by FEProblemBase, so it is not released in this
78  // object
79  std::unique_ptr<CouplingMatrix> cm = libmesh_make_unique<CouplingMatrix>(n_vars);
80 
81  bool full = false; // getParam<bool>("full"); // TODO: add a FULL option for PBP
82 
83  if (!full)
84  {
85  // put 1s on diagonal
86  for (unsigned int i = 0; i < n_vars; i++)
87  (*cm)(i, i) = 1;
88 
89  // off-diagonal entries
90  std::vector<std::vector<unsigned int>> off_diag(n_vars);
91  for (unsigned int i = 0; i < getParam<std::vector<std::string>>("off_diag_row").size(); i++)
92  {
93  unsigned int row =
94  nl.getVariable(0, getParam<std::vector<std::string>>("off_diag_row")[i]).number();
95  unsigned int column =
96  nl.getVariable(0, getParam<std::vector<std::string>>("off_diag_column")[i]).number();
97  (*cm)(row, column) = 1;
98  }
99 
100  // TODO: handle coupling entries between NL-vars and SCALAR-vars
101  }
102  else
103  {
104  for (unsigned int i = 0; i < n_vars; i++)
105  for (unsigned int j = 0; j < n_vars; j++)
106  (*cm)(i, j) = 1;
107  }
108 
109  _fe_problem.setCouplingMatrix(std::move(cm));
110  }
111 
112  // PC types
113  const std::vector<std::string> & pc_types = getParam<std::vector<std::string>>("preconditioner");
114  for (unsigned int i = 0; i < num_systems; i++)
115  _pre_type[i] = Utility::string_to_enum<PreconditionerType>(pc_types[i]);
116 
117  // solve order
118  const std::vector<std::string> & solve_order = getParam<std::vector<std::string>>("solve_order");
119  _solve_order.resize(solve_order.size());
120  for (unsigned int i = 0; i < solve_order.size(); i++)
121  _solve_order[i] = _nl.system().variable_number(solve_order[i]);
122 
123  // diag and off-diag systems
124  unsigned int n_vars = _nl.system().n_vars();
125 
126  // off-diagonal entries
127  const std::vector<std::string> & odr = getParam<std::vector<std::string>>("off_diag_row");
128  const std::vector<std::string> & odc = getParam<std::vector<std::string>>("off_diag_column");
129  std::vector<std::vector<unsigned int>> off_diag(n_vars);
130  for (unsigned int i = 0; i < odr.size(); i++)
131  {
132  unsigned int row = _nl.system().variable_number(odr[i]);
133  unsigned int column = _nl.system().variable_number(odc[i]);
134  off_diag[row].push_back(column);
135  }
136  // Add all of the preconditioning systems
137  for (unsigned int var = 0; var < n_vars; var++)
138  addSystem(var, off_diag[var], _pre_type[var]);
139 
140  _nl.nonlinearSolver()->attach_preconditioner(this);
141 
143  mooseError("PBP must be used with JFNK solve type");
144 }
145 
147 {
148  this->clear();
149 
150  for (auto & pc : _preconditioners)
151  delete pc;
152 }
153 
154 void
156  std::vector<unsigned int> off_diag,
157  PreconditionerType type)
158 {
159  std::string var_name = _nl.system().variable_name(var);
160 
161  LinearImplicitSystem & precond_system =
162  _fe_problem.es().add_system<LinearImplicitSystem>(var_name + "_system");
163  precond_system.assemble_before_solve = false;
164 
165  const std::set<SubdomainID> * active_subdomains = _nl.getVariableBlocks(var);
166  precond_system.add_variable(
167  var_name + "_prec", _nl.system().variable(var).type(), active_subdomains);
168 
169  _systems[var] = &precond_system;
170  _pre_type[var] = type;
171 
172  _off_diag_mats[var].resize(off_diag.size());
173  for (unsigned int i = 0; i < off_diag.size(); i++)
174  {
175  // Add the matrix to hold the off-diagonal piece
176  _off_diag_mats[var][i] = &precond_system.add_matrix(_nl.system().variable_name(off_diag[i]));
177  }
178 
179  _off_diag[var] = off_diag;
180 }
181 
182 void
184 {
185  Moose::perf_log.push("init()", "PhysicsBasedPreconditioner");
186 
187  // Tell libMesh that this is initialized!
188  _is_initialized = true;
189 
190  const unsigned int num_systems = _systems.size();
191 
192  // If no order was specified, just solve them in increasing order
193  if (_solve_order.size() == 0)
194  {
195  _solve_order.resize(num_systems);
196  for (unsigned int i = 0; i < num_systems; i++)
197  _solve_order[i] = i;
198  }
199 
200  // Loop over variables
201  for (unsigned int system_var = 0; system_var < num_systems; system_var++)
202  {
203  LinearImplicitSystem & u_system = *_systems[system_var];
204 
205  if (!_preconditioners[system_var])
206  _preconditioners[system_var] =
207  Preconditioner<Number>::build(MoosePreconditioner::_communicator);
208 
209  // we have to explicitly set the matrix in the preconditioner, because h-adaptivity could have
210  // changed it and we have to work with the current one
211  Preconditioner<Number> * preconditioner = _preconditioners[system_var];
212  preconditioner->set_matrix(*u_system.matrix);
213  preconditioner->set_type(_pre_type[system_var]);
214 
215  preconditioner->init();
216  }
217 
218  Moose::perf_log.pop("init()", "PhysicsBasedPreconditioner");
219 }
220 
221 void
223 {
224  const unsigned int num_systems = _systems.size();
225 
226  std::vector<JacobianBlock *> blocks;
227 
228  // Loop over variables
229  for (unsigned int system_var = 0; system_var < num_systems; system_var++)
230  {
231  LinearImplicitSystem & u_system = *_systems[system_var];
232 
233  {
234  JacobianBlock * block = new JacobianBlock(u_system, *u_system.matrix, system_var, system_var);
235  blocks.push_back(block);
236  }
237 
238  for (unsigned int diag = 0; diag < _off_diag[system_var].size(); diag++)
239  {
240  unsigned int coupled_var = _off_diag[system_var][diag];
241  std::string coupled_name = _nl.system().variable_name(coupled_var);
242 
243  JacobianBlock * block =
244  new JacobianBlock(u_system, *_off_diag_mats[system_var][diag], system_var, coupled_var);
245  blocks.push_back(block);
246  }
247  }
248 
250 
251  // cleanup
252  for (auto & block : blocks)
253  delete block;
254 }
255 
256 void
257 PhysicsBasedPreconditioner::apply(const NumericVector<Number> & x, NumericVector<Number> & y)
258 {
259  Moose::perf_log.push("apply()", "PhysicsBasedPreconditioner");
260 
261  const unsigned int num_systems = _systems.size();
262 
263  MooseMesh & mesh = _fe_problem.mesh();
264 
265  // Zero out the solution vectors
266  for (unsigned int sys = 0; sys < num_systems; sys++)
267  _systems[sys]->solution->zero();
268 
269  // Loop over solve order
270  for (unsigned int i = 0; i < _solve_order.size(); i++)
271  {
272  unsigned int system_var = _solve_order[i];
273 
274  LinearImplicitSystem & u_system = *_systems[system_var];
275 
276  // Copy rhs from the big system into the small one
278  mesh, _nl.system().number(), system_var, x, u_system.number(), 0, *u_system.rhs);
279 
280  // Modify the RHS by subtracting off the matvecs of the solutions for the other preconditioning
281  // systems with the off diagonal blocks in this system.
282  for (unsigned int diag = 0; diag < _off_diag[system_var].size(); diag++)
283  {
284  unsigned int coupled_var = _off_diag[system_var][diag];
285  LinearImplicitSystem & coupled_system = *_systems[coupled_var];
286  SparseMatrix<Number> & off_diag = *_off_diag_mats[system_var][diag];
287  NumericVector<Number> & rhs = *u_system.rhs;
288 
289  // This next bit computes rhs -= A*coupled_solution
290  // It does what it does because there is no vector_mult_sub()
291  rhs.close();
292  rhs.scale(-1.0);
293  rhs.close();
294  off_diag.vector_mult_add(rhs, *coupled_system.solution);
295  rhs.close();
296  rhs.scale(-1.0);
297  rhs.close();
298  }
299 
300  // Apply the preconditioner to the small system
301  _preconditioners[system_var]->apply(*u_system.rhs, *u_system.solution);
302 
303  // Copy solution from small system into the big one
304  // copyVarValues(mesh,system,0,*u_system.solution,0,system_var,y);
305  }
306 
307  // Copy the solutions out
308  for (unsigned int system_var = 0; system_var < num_systems; system_var++)
309  {
310  LinearImplicitSystem & u_system = *_systems[system_var];
311 
313  mesh, u_system.number(), 0, *u_system.solution, _nl.system().number(), system_var, y);
314  }
315 
316  y.close();
317 
318  Moose::perf_log.pop("apply()", "PhysicsBasedPreconditioner");
319 }
320 
321 void
323 {
324 }
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:140
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.
NonlinearSystemBase & getNonlinearSystemBase()
virtual NonlinearSolver< Number > * nonlinearSolver()=0
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< Preconditioner< Number > * > _preconditioners
Holds one Preconditioner object per small system to solve.
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 setCouplingMatrix(std::unique_ptr< CouplingMatrix > cm)
Set custom coupling matrix.
FEProblemBase & _fe_problem
Subproblem this preconditioner is part of.
static PetscErrorCode Vec Mat Mat pc
static PetscErrorCode Vec x
InputParameters validParams< MoosePreconditioner >()
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.
NonlinearSystemBase * nl
Nonlinear system to be solved.
InputParameters validParams< PhysicsBasedPreconditioner >()
virtual MooseVariable & getVariable(THREAD_ID tid, const std::string &var_name)
Gets a reference to a variable of with specified name.
Definition: SystemBase.C:103
PhysicsBasedPreconditioner(const InputParameters &params)
Constructor.
Jacobian-Free Newton Krylov.
Definition: MooseTypes.h:248
MooseMesh wraps a libMesh::Mesh object and enhances its capabilities by caching additional data and s...
Definition: MooseMesh.h:74
PerfLog perf_log
Perflog to be used by applications.
std::vector< std::vector< SparseMatrix< Number > * > > _off_diag_mats
Holds pointers to the off-diagonal matrices.
Moose::SolveType _type
Definition: SolverParams.h:25
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...
const T & getParam(const std::string &name) const
Retrieve a parameter for the object.
Definition: MooseObject.h:122
MatType type
unsigned int number() const
Get variable number coming from libMesh.
virtual System & system() override
Get the reference to the libMesh system.
virtual MooseMesh & mesh() override
std::vector< unsigned int > _solve_order
Holds the order the blocks are solved for.
void addParam(const std::string &name, const S &value, const std::string &doc_string)
These methods add an option parameter and a documentation string to the InputParameters object...
void mooseError(Args &&...args) const
Definition: MooseObject.h:80
virtual void computeJacobianBlocks(std::vector< JacobianBlock * > &blocks)
Computes several Jacobian blocks simultaneously, summing their contributions into smaller preconditio...
virtual unsigned int nVariables()
Get the number of variables in this system.
Definition: SystemBase.C:580