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 "PhysicsBasedPreconditioner.h"
11 :
12 : // MOOSE includes
13 : #include "ComputeJacobianBlocksThread.h"
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 :
31 : using namespace libMesh;
32 :
33 : registerMooseObjectAliased("MooseApp", PhysicsBasedPreconditioner, "PBP");
34 :
35 : InputParameters
36 14505 : PhysicsBasedPreconditioner::validParams()
37 : {
38 14505 : InputParameters params = MoosePreconditioner::validParams();
39 :
40 14505 : params.addClassDescription("Physics-based preconditioner (PBP) allows individual physics to have "
41 : "their own preconditioner.");
42 :
43 14505 : params.addRequiredParam<std::vector<std::string>>(
44 : "solve_order",
45 : "The order the block rows will be solved in. Put the name of variables here "
46 : "to stand for solving that variable's block row. A variable may appear more "
47 : "than once (to create cylces if you like).");
48 14505 : params.addRequiredParam<std::vector<std::string>>("preconditioner", "TODO: docstring");
49 :
50 14505 : return params;
51 0 : }
52 :
53 120 : PhysicsBasedPreconditioner::PhysicsBasedPreconditioner(const InputParameters & params)
54 : : MoosePreconditioner(params),
55 : Preconditioner<Number>(MoosePreconditioner::_communicator),
56 120 : _nl(_fe_problem.getNonlinearSystemBase(_nl_sys_num))
57 : {
58 120 : const auto & libmesh_system = _nl.system();
59 120 : unsigned int num_systems = _nl.system().n_vars();
60 120 : _systems.resize(num_systems);
61 120 : _preconditioners.resize(num_systems);
62 120 : _off_diag.resize(num_systems);
63 120 : _off_diag_mats.resize(num_systems);
64 120 : _pre_type.resize(num_systems);
65 :
66 : { // Setup the Coupling Matrix so MOOSE knows what we're doing
67 120 : unsigned int n_vars = _nl.nVariables();
68 :
69 : // The coupling matrix is held and released by FEProblemBase, so it is not released in this
70 : // object
71 120 : std::unique_ptr<CouplingMatrix> cm = std::make_unique<CouplingMatrix>(n_vars);
72 :
73 120 : bool full = getParam<bool>("full");
74 :
75 120 : if (!full)
76 : {
77 : // put 1s on diagonal
78 464 : for (unsigned int i = 0; i < n_vars; i++)
79 344 : (*cm)(i, i) = 1;
80 :
81 : // off-diagonal entries
82 120 : for (const auto & off_diag : getParam<NonlinearVariableName, NonlinearVariableName>(
83 291 : "off_diag_row", "off_diag_column"))
84 : {
85 51 : const auto row = libmesh_system.variable_number(off_diag.first);
86 51 : const auto column = libmesh_system.variable_number(off_diag.second);
87 51 : (*cm)(row, column) = 1;
88 120 : }
89 :
90 : // TODO: handle coupling entries between NL-vars and SCALAR-vars
91 : }
92 : else
93 : {
94 0 : for (unsigned int i = 0; i < n_vars; i++)
95 0 : for (unsigned int j = 0; j < n_vars; j++)
96 0 : (*cm)(i, j) = 1;
97 : }
98 :
99 120 : setCouplingMatrix(std::move(cm));
100 120 : }
101 :
102 : // PC types
103 120 : const std::vector<std::string> & pc_types = getParam<std::vector<std::string>>("preconditioner");
104 464 : for (unsigned int i = 0; i < num_systems; i++)
105 344 : _pre_type[i] = Utility::string_to_enum<PreconditionerType>(pc_types[i]);
106 :
107 : // solve order
108 120 : const std::vector<std::string> & solve_order = getParam<std::vector<std::string>>("solve_order");
109 120 : _solve_order.resize(solve_order.size());
110 464 : for (const auto i : index_range(solve_order))
111 344 : _solve_order[i] = _nl.system().variable_number(solve_order[i]);
112 :
113 : // diag and off-diag systems
114 120 : unsigned int n_vars = _nl.system().n_vars();
115 :
116 : // off-diagonal entries
117 120 : std::vector<std::vector<unsigned int>> off_diag(n_vars);
118 120 : if (isParamValid("off_diag_row") && isParamValid("off_diag_column"))
119 : {
120 : const std::vector<NonlinearVariableName> & odr =
121 44 : getParam<std::vector<NonlinearVariableName>>("off_diag_row");
122 : const std::vector<NonlinearVariableName> & odc =
123 44 : getParam<std::vector<NonlinearVariableName>>("off_diag_column");
124 :
125 95 : for (const auto i : index_range(odr))
126 : {
127 51 : unsigned int row = _nl.system().variable_number(odr[i]);
128 51 : unsigned int column = _nl.system().variable_number(odc[i]);
129 51 : off_diag[row].push_back(column);
130 : }
131 : }
132 : // Add all of the preconditioning systems
133 464 : for (unsigned int var = 0; var < n_vars; var++)
134 344 : addSystem(var, off_diag[var], _pre_type[var]);
135 :
136 120 : _nl.attachPreconditioner(this);
137 :
138 120 : if (_fe_problem.solverParams(_nl.number())._type != Moose::ST_JFNK)
139 0 : mooseError("PBP must be used with JFNK solve type");
140 120 : }
141 :
142 226 : PhysicsBasedPreconditioner::~PhysicsBasedPreconditioner() { this->clear(); }
143 :
144 : void
145 344 : PhysicsBasedPreconditioner::addSystem(unsigned int var,
146 : std::vector<unsigned int> off_diag,
147 : PreconditionerType type)
148 : {
149 344 : std::string var_name = _nl.system().variable_name(var);
150 :
151 : LinearImplicitSystem & precond_system =
152 344 : _fe_problem.es().add_system<LinearImplicitSystem>(var_name + "_system");
153 344 : precond_system.assemble_before_solve = false;
154 :
155 344 : const std::set<SubdomainID> * active_subdomains = _nl.getVariableBlocks(var);
156 344 : precond_system.add_variable(
157 688 : var_name + "_prec", _nl.system().variable(var).type(), active_subdomains);
158 :
159 344 : _systems[var] = &precond_system;
160 344 : _pre_type[var] = type;
161 :
162 344 : _off_diag_mats[var].resize(off_diag.size());
163 395 : for (const auto i : index_range(off_diag))
164 : {
165 : // Add the matrix to hold the off-diagonal piece
166 51 : _off_diag_mats[var][i] = &precond_system.add_matrix(_nl.system().variable_name(off_diag[i]));
167 : }
168 :
169 344 : _off_diag[var] = off_diag;
170 344 : }
171 :
172 : void
173 290 : PhysicsBasedPreconditioner::init()
174 : {
175 290 : TIME_SECTION("init", 2, "Initializing PhysicsBasedPreconditioner");
176 :
177 : // Tell libMesh that this is initialized!
178 290 : _is_initialized = true;
179 :
180 290 : const unsigned int num_systems = _systems.size();
181 :
182 : // If no order was specified, just solve them in increasing order
183 290 : if (_solve_order.size() == 0)
184 : {
185 0 : _solve_order.resize(num_systems);
186 0 : for (unsigned int i = 0; i < num_systems; i++)
187 0 : _solve_order[i] = i;
188 : }
189 :
190 : // Loop over variables
191 1062 : for (unsigned int system_var = 0; system_var < num_systems; system_var++)
192 : {
193 772 : LinearImplicitSystem & u_system = *_systems[system_var];
194 :
195 772 : if (!_preconditioners[system_var])
196 316 : _preconditioners[system_var] =
197 632 : Preconditioner<Number>::build_preconditioner(MoosePreconditioner::_communicator);
198 :
199 : // we have to explicitly set the matrix in the preconditioner, because h-adaptivity could have
200 : // changed it and we have to work with the current one
201 772 : Preconditioner<Number> * preconditioner = _preconditioners[system_var].get();
202 772 : preconditioner->set_matrix(u_system.get_system_matrix());
203 772 : preconditioner->set_type(_pre_type[system_var]);
204 :
205 772 : preconditioner->init();
206 : }
207 290 : }
208 :
209 : void
210 588 : PhysicsBasedPreconditioner::setup()
211 : {
212 588 : const unsigned int num_systems = _systems.size();
213 :
214 588 : std::vector<JacobianBlock *> blocks;
215 :
216 : // Loop over variables
217 1956 : for (unsigned int system_var = 0; system_var < num_systems; system_var++)
218 : {
219 1368 : LinearImplicitSystem & u_system = *_systems[system_var];
220 :
221 : {
222 : JacobianBlock * block =
223 1368 : new JacobianBlock(u_system, u_system.get_system_matrix(), system_var, system_var);
224 1368 : blocks.push_back(block);
225 : }
226 :
227 1512 : for (const auto diag : index_range(_off_diag[system_var]))
228 : {
229 144 : unsigned int coupled_var = _off_diag[system_var][diag];
230 144 : std::string coupled_name = _nl.system().variable_name(coupled_var);
231 :
232 : JacobianBlock * block =
233 144 : new JacobianBlock(u_system, *_off_diag_mats[system_var][diag], system_var, coupled_var);
234 144 : blocks.push_back(block);
235 144 : }
236 : }
237 :
238 588 : _fe_problem.computeJacobianBlocks(blocks, _nl.number());
239 :
240 : // cleanup
241 2100 : for (auto & block : blocks)
242 1512 : delete block;
243 588 : }
244 :
245 : void
246 2536 : PhysicsBasedPreconditioner::apply(const NumericVector<Number> & x, NumericVector<Number> & y)
247 : {
248 2536 : TIME_SECTION("apply", 1, "Applying PhysicsBasedPreconditioner");
249 :
250 2536 : const unsigned int num_systems = _systems.size();
251 :
252 2536 : MooseMesh & mesh = _fe_problem.mesh();
253 :
254 : // Zero out the solution vectors
255 8568 : for (unsigned int sys = 0; sys < num_systems; sys++)
256 6032 : _systems[sys]->solution->zero();
257 :
258 : // Loop over solve order
259 8568 : for (const auto i : index_range(_solve_order))
260 : {
261 6032 : unsigned int system_var = _solve_order[i];
262 :
263 6032 : LinearImplicitSystem & u_system = *_systems[system_var];
264 :
265 : // Copy rhs from the big system into the small one
266 6032 : MoosePreconditioner::copyVarValues(
267 6032 : mesh, _nl.system().number(), system_var, x, u_system.number(), 0, *u_system.rhs);
268 :
269 : // Modify the RHS by subtracting off the matvecs of the solutions for the other preconditioning
270 : // systems with the off diagonal blocks in this system.
271 6918 : for (const auto diag : index_range(_off_diag[system_var]))
272 : {
273 886 : unsigned int coupled_var = _off_diag[system_var][diag];
274 886 : LinearImplicitSystem & coupled_system = *_systems[coupled_var];
275 886 : SparseMatrix<Number> & off_diag = *_off_diag_mats[system_var][diag];
276 886 : NumericVector<Number> & rhs = *u_system.rhs;
277 :
278 : // This next bit computes rhs -= A*coupled_solution
279 : // It does what it does because there is no vector_mult_sub()
280 886 : rhs.close();
281 886 : rhs.scale(-1.0);
282 886 : rhs.close();
283 886 : off_diag.vector_mult_add(rhs, *coupled_system.solution);
284 886 : rhs.close();
285 886 : rhs.scale(-1.0);
286 886 : rhs.close();
287 : }
288 :
289 : // If there is no off_diag, then u_system.rhs will not be closed.
290 : // Thus, we need to close it right here
291 6032 : if (!_off_diag[system_var].size())
292 5146 : u_system.rhs->close();
293 :
294 : // Apply the preconditioner to the small system
295 6032 : _preconditioners[system_var]->apply(*u_system.rhs, *u_system.solution);
296 :
297 : // Copy solution from small system into the big one
298 : // copyVarValues(mesh,system,0,*u_system.solution,0,system_var,y);
299 : }
300 :
301 : // Copy the solutions out
302 8568 : for (unsigned int system_var = 0; system_var < num_systems; system_var++)
303 : {
304 6032 : LinearImplicitSystem & u_system = *_systems[system_var];
305 :
306 6032 : MoosePreconditioner::copyVarValues(
307 6032 : mesh, u_system.number(), 0, *u_system.solution, _nl.system().number(), system_var, y);
308 : }
309 :
310 2536 : y.close();
311 2536 : }
312 :
313 : void
314 113 : PhysicsBasedPreconditioner::clear()
315 : {
316 113 : }
|
