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 
33 template <>
36 {
38 
39  params.addRequiredParam<std::vector<std::string>>(
40  "solve_order",
41  "The order the block rows will be solved in. Put the name of variables here "
42  "to stand for solving that variable's block row. A variable may appear more "
43  "than once (to create cylces if you like).");
44  params.addRequiredParam<std::vector<std::string>>("preconditioner", "TODO: docstring");
45 
46  params.addParam<std::vector<std::string>>(
47  "off_diag_row",
48  "The off diagonal row you want to add into the matrix, it will be associated "
49  "with an off diagonal column from the same position in off_diag_colum.");
50  params.addParam<std::vector<std::string>>("off_diag_column",
51  "The off diagonal column you want to add into the "
52  "matrix, it will be associated with an off diagonal "
53  "row from the same position in off_diag_row.");
54 
55  return params;
56 }
57 
59  : MoosePreconditioner(params),
60  Preconditioner<Number>(MoosePreconditioner::_communicator),
61  _nl(_fe_problem.getNonlinearSystemBase()),
62  _init_timer(registerTimedSection("init", 2)),
63  _apply_timer(registerTimedSection("apply", 1))
64 {
65  unsigned int num_systems = _nl.system().n_vars();
66  _systems.resize(num_systems);
67  _preconditioners.resize(num_systems);
68  _off_diag.resize(num_systems);
69  _off_diag_mats.resize(num_systems);
70  _pre_type.resize(num_systems);
71 
72  { // Setup the Coupling Matrix so MOOSE knows what we're doing
74  unsigned int n_vars = nl.nVariables();
75 
76  // The coupling matrix is held and released by FEProblemBase, so it is not released in this
77  // object
78  std::unique_ptr<CouplingMatrix> cm = libmesh_make_unique<CouplingMatrix>(n_vars);
79 
80  bool full = false; // getParam<bool>("full"); // TODO: add a FULL option for PBP
81 
82  if (!full)
83  {
84  // put 1s on diagonal
85  for (unsigned int i = 0; i < n_vars; i++)
86  (*cm)(i, i) = 1;
87 
88  // off-diagonal entries
89  std::vector<std::vector<unsigned int>> off_diag(n_vars);
90  for (unsigned int i = 0; i < getParam<std::vector<std::string>>("off_diag_row").size(); i++)
91  {
92  unsigned int row =
93  nl.getVariable(0, getParam<std::vector<std::string>>("off_diag_row")[i]).number();
94  unsigned int column =
95  nl.getVariable(0, getParam<std::vector<std::string>>("off_diag_column")[i]).number();
96  (*cm)(row, column) = 1;
97  }
98 
99  // TODO: handle coupling entries between NL-vars and SCALAR-vars
100  }
101  else
102  {
103  for (unsigned int i = 0; i < n_vars; i++)
104  for (unsigned int j = 0; j < n_vars; j++)
105  (*cm)(i, j) = 1;
106  }
107 
108  _fe_problem.setCouplingMatrix(std::move(cm));
109  }
110 
111  // PC types
112  const std::vector<std::string> & pc_types = getParam<std::vector<std::string>>("preconditioner");
113  for (unsigned int i = 0; i < num_systems; i++)
114  _pre_type[i] = Utility::string_to_enum<PreconditionerType>(pc_types[i]);
115 
116  // solve order
117  const std::vector<std::string> & solve_order = getParam<std::vector<std::string>>("solve_order");
118  _solve_order.resize(solve_order.size());
119  for (unsigned int i = 0; i < solve_order.size(); i++)
120  _solve_order[i] = _nl.system().variable_number(solve_order[i]);
121 
122  // diag and off-diag systems
123  unsigned int n_vars = _nl.system().n_vars();
124 
125  // off-diagonal entries
126  const std::vector<std::string> & odr = getParam<std::vector<std::string>>("off_diag_row");
127  const std::vector<std::string> & odc = getParam<std::vector<std::string>>("off_diag_column");
128  std::vector<std::vector<unsigned int>> off_diag(n_vars);
129  for (unsigned int i = 0; i < odr.size(); i++)
130  {
131  unsigned int row = _nl.system().variable_number(odr[i]);
132  unsigned int column = _nl.system().variable_number(odc[i]);
133  off_diag[row].push_back(column);
134  }
135  // Add all of the preconditioning systems
136  for (unsigned int var = 0; var < n_vars; var++)
137  addSystem(var, off_diag[var], _pre_type[var]);
138 
139  _nl.nonlinearSolver()->attach_preconditioner(this);
140 
142  mooseError("PBP must be used with JFNK solve type");
143 }
144 
146 
147 void
149  std::vector<unsigned int> off_diag,
150  PreconditionerType type)
151 {
152  std::string var_name = _nl.system().variable_name(var);
153 
154  LinearImplicitSystem & precond_system =
155  _fe_problem.es().add_system<LinearImplicitSystem>(var_name + "_system");
156  precond_system.assemble_before_solve = false;
157 
158  const std::set<SubdomainID> * active_subdomains = _nl.getVariableBlocks(var);
159  precond_system.add_variable(
160  var_name + "_prec", _nl.system().variable(var).type(), active_subdomains);
161 
162  _systems[var] = &precond_system;
163  _pre_type[var] = type;
164 
165  _off_diag_mats[var].resize(off_diag.size());
166  for (unsigned int i = 0; i < off_diag.size(); i++)
167  {
168  // Add the matrix to hold the off-diagonal piece
169  _off_diag_mats[var][i] = &precond_system.add_matrix(_nl.system().variable_name(off_diag[i]));
170  }
171 
172  _off_diag[var] = off_diag;
173 }
174 
175 void
177 {
178  TIME_SECTION(_init_timer);
179 
180  // Tell libMesh that this is initialized!
181  _is_initialized = true;
182 
183  const unsigned int num_systems = _systems.size();
184 
185  // If no order was specified, just solve them in increasing order
186  if (_solve_order.size() == 0)
187  {
188  _solve_order.resize(num_systems);
189  for (unsigned int i = 0; i < num_systems; i++)
190  _solve_order[i] = i;
191  }
192 
193  // Loop over variables
194  for (unsigned int system_var = 0; system_var < num_systems; system_var++)
195  {
196  LinearImplicitSystem & u_system = *_systems[system_var];
197 
198  if (!_preconditioners[system_var])
199  _preconditioners[system_var] =
200  Preconditioner<Number>::build_preconditioner(MoosePreconditioner::_communicator);
201 
202  // we have to explicitly set the matrix in the preconditioner, because h-adaptivity could have
203  // changed it and we have to work with the current one
204  Preconditioner<Number> * preconditioner = _preconditioners[system_var].get();
205  preconditioner->set_matrix(*u_system.matrix);
206  preconditioner->set_type(_pre_type[system_var]);
207 
208  preconditioner->init();
209  }
210 }
211 
212 void
214 {
215  const unsigned int num_systems = _systems.size();
216 
217  std::vector<JacobianBlock *> blocks;
218 
219  // Loop over variables
220  for (unsigned int system_var = 0; system_var < num_systems; system_var++)
221  {
222  LinearImplicitSystem & u_system = *_systems[system_var];
223 
224  {
225  JacobianBlock * block = new JacobianBlock(u_system, *u_system.matrix, system_var, system_var);
226  blocks.push_back(block);
227  }
228 
229  for (unsigned int diag = 0; diag < _off_diag[system_var].size(); diag++)
230  {
231  unsigned int coupled_var = _off_diag[system_var][diag];
232  std::string coupled_name = _nl.system().variable_name(coupled_var);
233 
234  JacobianBlock * block =
235  new JacobianBlock(u_system, *_off_diag_mats[system_var][diag], system_var, coupled_var);
236  blocks.push_back(block);
237  }
238  }
239 
241 
242  // cleanup
243  for (auto & block : blocks)
244  delete block;
245 }
246 
247 void
248 PhysicsBasedPreconditioner::apply(const NumericVector<Number> & x, NumericVector<Number> & y)
249 {
250  TIME_SECTION(_apply_timer);
251 
252  const unsigned int num_systems = _systems.size();
253 
254  MooseMesh & mesh = _fe_problem.mesh();
255 
256  // Zero out the solution vectors
257  for (unsigned int sys = 0; sys < num_systems; sys++)
258  _systems[sys]->solution->zero();
259 
260  // Loop over solve order
261  for (unsigned int i = 0; i < _solve_order.size(); i++)
262  {
263  unsigned int system_var = _solve_order[i];
264 
265  LinearImplicitSystem & u_system = *_systems[system_var];
266 
267  // Copy rhs from the big system into the small one
269  mesh, _nl.system().number(), system_var, x, u_system.number(), 0, *u_system.rhs);
270 
271  // Modify the RHS by subtracting off the matvecs of the solutions for the other preconditioning
272  // systems with the off diagonal blocks in this system.
273  for (unsigned int diag = 0; diag < _off_diag[system_var].size(); diag++)
274  {
275  unsigned int coupled_var = _off_diag[system_var][diag];
276  LinearImplicitSystem & coupled_system = *_systems[coupled_var];
277  SparseMatrix<Number> & off_diag = *_off_diag_mats[system_var][diag];
278  NumericVector<Number> & rhs = *u_system.rhs;
279 
280  // This next bit computes rhs -= A*coupled_solution
281  // It does what it does because there is no vector_mult_sub()
282  rhs.close();
283  rhs.scale(-1.0);
284  rhs.close();
285  off_diag.vector_mult_add(rhs, *coupled_system.solution);
286  rhs.close();
287  rhs.scale(-1.0);
288  rhs.close();
289  }
290 
291  // Apply the preconditioner to the small system
292  _preconditioners[system_var]->apply(*u_system.rhs, *u_system.solution);
293 
294  // Copy solution from small system into the big one
295  // copyVarValues(mesh,system,0,*u_system.solution,0,system_var,y);
296  }
297 
298  // Copy the solutions out
299  for (unsigned int system_var = 0; system_var < num_systems; system_var++)
300  {
301  LinearImplicitSystem & u_system = *_systems[system_var];
302 
304  mesh, u_system.number(), 0, *u_system.solution, _nl.system().number(), system_var, y);
305  }
306 
307  y.close();
308 }
309 
310 void
312 {
313 }
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:168
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
unsigned int number() const
Get variable number coming from libMesh.
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 setCouplingMatrix(std::unique_ptr< CouplingMatrix > cm)
Set custom coupling matrix.
FEProblemBase & _fe_problem
Subproblem this preconditioner is part of.
void mooseError(Args &&... args) const
Definition: MooseObject.h:147
static PetscErrorCode Vec x
InputParameters validParams< MoosePreconditioner >()
virtual void clear()
Release all memory and clear data structures.
const T & getParam(const std::string &name) const
Retrieve a parameter for the object.
Definition: MooseObject.h:191
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 void computeJacobianBlocks(std::vector< JacobianBlock *> &blocks)
Computes several Jacobian blocks simultaneously, summing their contributions into smaller preconditio...
virtual EquationSystems & es() override
virtual void init()
Initialize data structures if not done so already.
const std::string & type() const
Get the type of this object.
Definition: MooseObject.h:53
virtual unsigned int nVariables() const
Get the number of variables in this system.
Definition: SystemBase.C:692
NonlinearSystemBase * nl
Nonlinear system to be solved.
MooseVariableFEBase & getVariable(THREAD_ID tid, const std::string &var_name)
Gets a reference to a variable of with specified name.
Definition: SystemBase.C:105
InputParameters validParams< PhysicsBasedPreconditioner >()
PhysicsBasedPreconditioner(const InputParameters &params)
Constructor.
Jacobian-Free Newton Krylov.
Definition: MooseTypes.h:592
MooseMesh wraps a libMesh::Mesh object and enhances its capabilities by caching additional data and s...
Definition: MooseMesh.h:74
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...
MatType type
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 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...
registerMooseObjectAliased("MooseApp", PhysicsBasedPreconditioner, "PBP")