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 14497 : PhysicsBasedPreconditioner::validParams()
37 : {
38 14497 : InputParameters params = MoosePreconditioner::validParams();
39 :
40 14497 : params.addClassDescription("Physics-based preconditioner (PBP) allows individual physics to have "
41 : "their own preconditioner.");
42 :
43 14497 : 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 14497 : params.addRequiredParam<std::vector<std::string>>("preconditioner", "TODO: docstring");
49 :
50 14497 : return params;
51 0 : }
52 :
53 116 : PhysicsBasedPreconditioner::PhysicsBasedPreconditioner(const InputParameters & params)
54 : : MoosePreconditioner(params),
55 : Preconditioner<Number>(MoosePreconditioner::_communicator),
56 116 : _nl(_fe_problem.getNonlinearSystemBase(_nl_sys_num))
57 : {
58 116 : const auto & libmesh_system = _nl.system();
59 116 : unsigned int num_systems = _nl.system().n_vars();
60 116 : _systems.resize(num_systems);
61 116 : _preconditioners.resize(num_systems);
62 116 : _off_diag.resize(num_systems);
63 116 : _off_diag_mats.resize(num_systems);
64 116 : _pre_type.resize(num_systems);
65 :
66 : { // Setup the Coupling Matrix so MOOSE knows what we're doing
67 116 : 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 116 : std::unique_ptr<CouplingMatrix> cm = std::make_unique<CouplingMatrix>(n_vars);
72 :
73 116 : bool full = getParam<bool>("full");
74 :
75 116 : if (!full)
76 : {
77 : // put 1s on diagonal
78 444 : for (unsigned int i = 0; i < n_vars; i++)
79 328 : (*cm)(i, i) = 1;
80 :
81 : // off-diagonal entries
82 116 : for (const auto & off_diag : getParam<NonlinearVariableName, NonlinearVariableName>(
83 281 : "off_diag_row", "off_diag_column"))
84 : {
85 49 : const auto row = libmesh_system.variable_number(off_diag.first);
86 49 : const auto column = libmesh_system.variable_number(off_diag.second);
87 49 : (*cm)(row, column) = 1;
88 116 : }
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 116 : setCouplingMatrix(std::move(cm));
100 116 : }
101 :
102 : // PC types
103 116 : const std::vector<std::string> & pc_types = getParam<std::vector<std::string>>("preconditioner");
104 444 : for (unsigned int i = 0; i < num_systems; i++)
105 328 : _pre_type[i] = Utility::string_to_enum<PreconditionerType>(pc_types[i]);
106 :
107 : // solve order
108 116 : const std::vector<std::string> & solve_order = getParam<std::vector<std::string>>("solve_order");
109 116 : _solve_order.resize(solve_order.size());
110 444 : for (const auto i : index_range(solve_order))
111 328 : _solve_order[i] = _nl.system().variable_number(solve_order[i]);
112 :
113 : // diag and off-diag systems
114 116 : unsigned int n_vars = _nl.system().n_vars();
115 :
116 : // off-diagonal entries
117 116 : std::vector<std::vector<unsigned int>> off_diag(n_vars);
118 116 : if (isParamValid("off_diag_row") && isParamValid("off_diag_column"))
119 : {
120 : const std::vector<NonlinearVariableName> & odr =
121 42 : getParam<std::vector<NonlinearVariableName>>("off_diag_row");
122 : const std::vector<NonlinearVariableName> & odc =
123 42 : getParam<std::vector<NonlinearVariableName>>("off_diag_column");
124 :
125 91 : for (const auto i : index_range(odr))
126 : {
127 49 : unsigned int row = _nl.system().variable_number(odr[i]);
128 49 : unsigned int column = _nl.system().variable_number(odc[i]);
129 49 : off_diag[row].push_back(column);
130 : }
131 : }
132 : // Add all of the preconditioning systems
133 444 : for (unsigned int var = 0; var < n_vars; var++)
134 328 : addSystem(var, off_diag[var], _pre_type[var]);
135 :
136 116 : _nl.attachPreconditioner(this);
137 :
138 116 : if (_fe_problem.solverParams(_nl.number())._type != Moose::ST_JFNK)
139 0 : mooseError("PBP must be used with JFNK solve type");
140 116 : }
141 :
142 218 : PhysicsBasedPreconditioner::~PhysicsBasedPreconditioner() { this->clear(); }
143 :
144 : void
145 328 : PhysicsBasedPreconditioner::addSystem(unsigned int var,
146 : std::vector<unsigned int> off_diag,
147 : PreconditionerType type)
148 : {
149 328 : std::string var_name = _nl.system().variable_name(var);
150 :
151 : LinearImplicitSystem & precond_system =
152 328 : _fe_problem.es().add_system<LinearImplicitSystem>(var_name + "_system");
153 328 : precond_system.assemble_before_solve = false;
154 :
155 328 : const std::set<SubdomainID> * active_subdomains = _nl.getVariableBlocks(var);
156 328 : precond_system.add_variable(
157 656 : var_name + "_prec", _nl.system().variable(var).type(), active_subdomains);
158 :
159 328 : _systems[var] = &precond_system;
160 328 : _pre_type[var] = type;
161 :
162 328 : _off_diag_mats[var].resize(off_diag.size());
163 377 : for (const auto i : index_range(off_diag))
164 : {
165 : // Add the matrix to hold the off-diagonal piece
166 49 : _off_diag_mats[var][i] = &precond_system.add_matrix(_nl.system().variable_name(off_diag[i]));
167 : }
168 :
169 328 : _off_diag[var] = off_diag;
170 328 : }
171 :
172 : void
173 270 : PhysicsBasedPreconditioner::init()
174 : {
175 270 : TIME_SECTION("init", 2, "Initializing PhysicsBasedPreconditioner");
176 :
177 : // Tell libMesh that this is initialized!
178 270 : _is_initialized = true;
179 :
180 270 : const unsigned int num_systems = _systems.size();
181 :
182 : // If no order was specified, just solve them in increasing order
183 270 : 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 986 : for (unsigned int system_var = 0; system_var < num_systems; system_var++)
192 : {
193 716 : LinearImplicitSystem & u_system = *_systems[system_var];
194 :
195 716 : if (!_preconditioners[system_var])
196 300 : _preconditioners[system_var] =
197 600 : 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 716 : Preconditioner<Number> * preconditioner = _preconditioners[system_var].get();
202 716 : preconditioner->set_matrix(u_system.get_system_matrix());
203 716 : preconditioner->set_type(_pre_type[system_var]);
204 :
205 716 : preconditioner->init();
206 : }
207 270 : }
208 :
209 : void
210 569 : PhysicsBasedPreconditioner::setup()
211 : {
212 569 : const unsigned int num_systems = _systems.size();
213 :
214 569 : std::vector<JacobianBlock *> blocks;
215 :
216 : // Loop over variables
217 1883 : for (unsigned int system_var = 0; system_var < num_systems; system_var++)
218 : {
219 1314 : LinearImplicitSystem & u_system = *_systems[system_var];
220 :
221 : {
222 : JacobianBlock * block =
223 1314 : new JacobianBlock(u_system, u_system.get_system_matrix(), system_var, system_var);
224 1314 : blocks.push_back(block);
225 : }
226 :
227 1449 : for (const auto diag : index_range(_off_diag[system_var]))
228 : {
229 135 : unsigned int coupled_var = _off_diag[system_var][diag];
230 135 : std::string coupled_name = _nl.system().variable_name(coupled_var);
231 :
232 : JacobianBlock * block =
233 135 : new JacobianBlock(u_system, *_off_diag_mats[system_var][diag], system_var, coupled_var);
234 135 : blocks.push_back(block);
235 135 : }
236 : }
237 :
238 569 : _fe_problem.computeJacobianBlocks(blocks, _nl.number());
239 :
240 : // cleanup
241 2018 : for (auto & block : blocks)
242 1449 : delete block;
243 569 : }
244 :
245 : void
246 2424 : PhysicsBasedPreconditioner::apply(const NumericVector<Number> & x, NumericVector<Number> & y)
247 : {
248 2424 : TIME_SECTION("apply", 1, "Applying PhysicsBasedPreconditioner");
249 :
250 2424 : const unsigned int num_systems = _systems.size();
251 :
252 2424 : MooseMesh & mesh = _fe_problem.mesh();
253 :
254 : // Zero out the solution vectors
255 8152 : for (unsigned int sys = 0; sys < num_systems; sys++)
256 5728 : _systems[sys]->solution->zero();
257 :
258 : // Loop over solve order
259 8152 : for (const auto i : index_range(_solve_order))
260 : {
261 5728 : unsigned int system_var = _solve_order[i];
262 :
263 5728 : LinearImplicitSystem & u_system = *_systems[system_var];
264 :
265 : // Copy rhs from the big system into the small one
266 5728 : MoosePreconditioner::copyVarValues(
267 5728 : 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 6552 : for (const auto diag : index_range(_off_diag[system_var]))
272 : {
273 824 : unsigned int coupled_var = _off_diag[system_var][diag];
274 824 : LinearImplicitSystem & coupled_system = *_systems[coupled_var];
275 824 : SparseMatrix<Number> & off_diag = *_off_diag_mats[system_var][diag];
276 824 : 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 824 : rhs.close();
281 824 : rhs.scale(-1.0);
282 824 : rhs.close();
283 824 : off_diag.vector_mult_add(rhs, *coupled_system.solution);
284 824 : rhs.close();
285 824 : rhs.scale(-1.0);
286 824 : 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 5728 : if (!_off_diag[system_var].size())
292 4904 : u_system.rhs->close();
293 :
294 : // Apply the preconditioner to the small system
295 5728 : _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 8152 : for (unsigned int system_var = 0; system_var < num_systems; system_var++)
303 : {
304 5728 : LinearImplicitSystem & u_system = *_systems[system_var];
305 :
306 5728 : MoosePreconditioner::copyVarValues(
307 5728 : mesh, u_system.number(), 0, *u_system.solution, _nl.system().number(), system_var, y);
308 : }
309 :
310 2424 : y.close();
311 2424 : }
312 :
313 : void
314 109 : PhysicsBasedPreconditioner::clear()
315 : {
316 109 : }
|
