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 "VariableCondensationPreconditioner.h"
11 :
12 : // MOOSE includes
13 : #include "FEProblem.h"
14 : #include "MooseUtils.h"
15 : #include "MooseVariableFE.h"
16 : #include "NonlinearSystem.h"
17 : #include "ComputeJacobianBlocksThread.h"
18 : #include "MooseEnum.h"
19 :
20 : #include "libmesh/coupling_matrix.h"
21 : #include "libmesh/libmesh_common.h"
22 : #include "libmesh/equation_systems.h"
23 : #include "libmesh/nonlinear_implicit_system.h"
24 : #include "libmesh/nonlinear_solver.h"
25 : #include "libmesh/linear_implicit_system.h"
26 : #include "libmesh/transient_system.h"
27 : #include "libmesh/numeric_vector.h"
28 : #include "libmesh/sparse_matrix.h"
29 : #include "libmesh/string_to_enum.h"
30 : #include "libmesh/mesh_base.h"
31 : #include "libmesh/variable.h"
32 : #include "libmesh/petsc_matrix.h"
33 : #include "libmesh/parallel_object.h"
34 : #include "libmesh/boundary_info.h"
35 :
36 : #include <petscmat.h>
37 :
38 : using namespace libMesh;
39 :
40 : registerMooseObjectAliased("MooseApp", VariableCondensationPreconditioner, "VCP");
41 :
42 : InputParameters
43 14345 : VariableCondensationPreconditioner::validParams()
44 : {
45 14345 : InputParameters params = MoosePreconditioner::validParams();
46 :
47 14345 : params.addClassDescription(
48 : "Variable condensation preconditioner (VCP) condenses out specified variable(s) "
49 : "from the Jacobian matrix and produces a system of equations with less unkowns to "
50 : "be solved by the underlying preconditioners.");
51 :
52 14345 : params.addParam<std::vector<NonlinearVariableName>>(
53 : "coupled_groups",
54 : {},
55 : "List multiple space separated groups of comma separated variables. "
56 : "Off-diagonal jacobians will be generated for all pairs within a group.");
57 :
58 43035 : params.addParam<bool>(
59 : "is_lm_coupling_diagonal",
60 28690 : false,
61 : "Set to true if you are sure the coupling matrix between Lagrange multiplier variable and "
62 : "the coupled primal variable is strict diagonal. This will speedup the linear solve. "
63 : "Otherwise set to false to ensure linear solve accuracy.");
64 43035 : params.addParam<bool>(
65 : "adaptive_condensation",
66 28690 : true,
67 : "By default VCP will check the Jacobian and only condense the rows with zero diagonals. Set "
68 : "to false if you want to condense out all the specified variable dofs.");
69 14345 : params.addRequiredParam<std::vector<std::string>>("preconditioner", "Preconditioner type.");
70 14345 : params.addRequiredParam<std::vector<std::string>>(
71 : "lm_variable",
72 : "Name of the variable(s) that is to be condensed out. Usually "
73 : "this will be the Lagrange multiplier variable(s).");
74 14345 : params.addRequiredParam<std::vector<std::string>>(
75 : "primary_variable",
76 : "Name of the variable(s) that couples with the variable(s) specified in the `variable` "
77 : "block. Usually this is the primary variable that the Lagrange multiplier correspond to.");
78 14345 : return params;
79 0 : }
80 :
81 40 : VariableCondensationPreconditioner::VariableCondensationPreconditioner(
82 40 : const InputParameters & params)
83 : : MoosePreconditioner(params),
84 : Preconditioner<Number>(MoosePreconditioner::_communicator),
85 80 : _nl(_fe_problem.getNonlinearSystemBase(_nl_sys_num)),
86 40 : _mesh(_fe_problem.mesh()),
87 40 : _dofmap(_nl.system().get_dof_map()),
88 40 : _is_lm_coupling_diagonal(getParam<bool>("is_lm_coupling_diagonal")),
89 40 : _adaptive_condensation(getParam<bool>("adaptive_condensation")),
90 40 : _n_vars(_nl.nVariables()),
91 40 : _lm_var_names(getParam<std::vector<std::string>>("lm_variable")),
92 40 : _primary_var_names(getParam<std::vector<std::string>>("primary_variable")),
93 40 : _D(std::make_unique<PetscMatrix<Number>>(MoosePreconditioner::_communicator)),
94 40 : _M(std::make_unique<PetscMatrix<Number>>(MoosePreconditioner::_communicator)),
95 40 : _K(std::make_unique<PetscMatrix<Number>>(MoosePreconditioner::_communicator)),
96 40 : _dinv(nullptr),
97 40 : _J_condensed(std::make_unique<PetscMatrix<Number>>(MoosePreconditioner::_communicator)),
98 40 : _x_hat(NumericVector<Number>::build(MoosePreconditioner::_communicator)),
99 40 : _y_hat(NumericVector<Number>::build(MoosePreconditioner::_communicator)),
100 40 : _primary_rhs_vec(NumericVector<Number>::build(MoosePreconditioner::_communicator)),
101 40 : _lm_sol_vec(NumericVector<Number>::build(MoosePreconditioner::_communicator)),
102 40 : _need_condense(true),
103 40 : _init_timer(registerTimedSection("init", 2)),
104 200 : _apply_timer(registerTimedSection("apply", 1))
105 : {
106 40 : if (_lm_var_names.size() != _primary_var_names.size())
107 0 : paramError("coupled_variable", "coupled_variable should have the same size as the variable.");
108 :
109 40 : if (!_mesh.getMesh().is_replicated())
110 0 : mooseError("The VariableCondensationPreconditioner cannot be used with DistributedMesh");
111 :
112 : // get variable ids from the variable names
113 80 : for (const auto & var_name : _lm_var_names)
114 : {
115 40 : if (!_nl.system().has_variable(var_name))
116 0 : paramError("variable ", var_name, " does not exist in the system");
117 40 : const unsigned int id = _nl.system().variable_number(var_name);
118 40 : _lm_var_ids.push_back(id);
119 : }
120 :
121 : // get coupled variable ids from the coupled variable names
122 80 : for (const auto & var_name : _primary_var_names)
123 : {
124 40 : if (!_nl.system().has_variable(var_name))
125 0 : paramError("coupled_variable ", var_name, " does not exist in the system");
126 40 : const unsigned int id = _nl.system().variable_number(var_name);
127 40 : _primary_var_ids.push_back(id);
128 : }
129 :
130 : // PC type
131 40 : const std::vector<std::string> & pc_type = getParam<std::vector<std::string>>("preconditioner");
132 40 : if (pc_type.size() > 1)
133 0 : mooseWarning("We only use one preconditioner type in VCP, the ",
134 0 : pc_type[0],
135 : " preconditioner is utilized.");
136 40 : _pre_type = Utility::string_to_enum<PreconditionerType>(pc_type[0]);
137 :
138 : // The following obtains and sets the coupling matrix.
139 : // TODO: This part can be refactored together with what are in other classes, e.g.,
140 : // PhysicsBasedPreconditioner
141 40 : std::unique_ptr<CouplingMatrix> cm = std::make_unique<CouplingMatrix>(_n_vars);
142 40 : const bool full = getParam<bool>("full");
143 :
144 40 : if (!full)
145 : {
146 : // put 1s on diagonal
147 0 : for (const auto i : make_range(_n_vars))
148 0 : (*cm)(i, i) = 1;
149 :
150 : // off-diagonal entries from the off_diag_row and off_diag_column parameters
151 0 : std::vector<std::vector<unsigned int>> off_diag(_n_vars);
152 0 : if (isParamValid("off_diag_row") && isParamValid("off_diag_column"))
153 :
154 0 : for (const auto i : index_range(getParam<std::vector<NonlinearVariableName>>("off_diag_row")))
155 : {
156 : const unsigned int row =
157 0 : _nl.getVariable(0, getParam<std::vector<NonlinearVariableName>>("off_diag_row")[i])
158 0 : .number();
159 : const unsigned int column =
160 0 : _nl.getVariable(0, getParam<std::vector<NonlinearVariableName>>("off_diag_column")[i])
161 0 : .number();
162 0 : (*cm)(row, column) = 1;
163 : }
164 :
165 : // off-diagonal entries from the coupled_groups parameters
166 0 : for (const auto & coupled_group :
167 0 : getParam<std::vector<NonlinearVariableName>>("coupled_groups"))
168 : {
169 0 : std::vector<NonlinearVariableName> vars;
170 0 : MooseUtils::tokenize<NonlinearVariableName>(coupled_group, vars, 1, ",");
171 0 : for (unsigned int j : index_range(vars))
172 0 : for (unsigned int k = j + 1; k < vars.size(); ++k)
173 : {
174 0 : const unsigned int row = _nl.getVariable(0, vars[j]).number();
175 0 : const unsigned int column = _nl.getVariable(0, vars[k]).number();
176 0 : (*cm)(row, column) = 1;
177 0 : (*cm)(column, row) = 1;
178 : }
179 0 : }
180 0 : }
181 : else
182 : {
183 120 : for (unsigned int i = 0; i < _n_vars; i++)
184 240 : for (unsigned int j = 0; j < _n_vars; j++)
185 160 : (*cm)(i, j) = 1;
186 : }
187 :
188 40 : setCouplingMatrix(std::move(cm));
189 :
190 40 : _nl.attachPreconditioner(this);
191 40 : }
192 :
193 80 : VariableCondensationPreconditioner::~VariableCondensationPreconditioner() { this->clear(); }
194 :
195 : void
196 64 : VariableCondensationPreconditioner::getDofToCondense()
197 : {
198 : // clean the containers if we want to update the dofs
199 64 : _global_lm_dofs.clear();
200 64 : _lm_dofs.clear();
201 64 : _global_primary_dofs.clear();
202 64 : _primary_dofs.clear();
203 64 : _map_global_lm_primary.clear();
204 64 : _map_global_primary_order.clear();
205 :
206 : // TODO: this might not work for distributed mesh and needs to be improved
207 64 : NodeRange * active_nodes = _mesh.getActiveNodeRange();
208 :
209 : // loop through the variable ids
210 64 : std::vector<dof_id_type> di, cp_di;
211 128 : for (const auto & vn : index_range(_lm_var_ids))
212 1692 : for (const auto & node : *active_nodes)
213 : {
214 1628 : di.clear();
215 1628 : cp_di.clear();
216 1628 : const auto var_id = _lm_var_ids[vn];
217 : // get coupled variable id
218 1628 : const auto cp_var_id = _primary_var_ids[vn];
219 : // get var and cp_var dofs associated with this node
220 1628 : _dofmap.dof_indices(node, di, var_id);
221 : // skip when di is empty
222 1628 : if (di.empty())
223 1318 : continue;
224 310 : _dofmap.dof_indices(node, cp_di, cp_var_id);
225 310 : if (cp_di.size() != di.size())
226 0 : mooseError("variable and coupled variable do not have the same number of dof on node ",
227 0 : node->id(),
228 : ".");
229 620 : for (const auto & i : index_range(di))
230 : {
231 : // when we have adaptive condensation, skip when di does not contain any indices in
232 : // _zero_rows
233 310 : if (std::find(_zero_rows.begin(), _zero_rows.end(), di[i]) == _zero_rows.end() &&
234 0 : _adaptive_condensation)
235 0 : break;
236 310 : _global_lm_dofs.push_back(di[i]);
237 310 : if (_dofmap.local_index(di[i]))
238 234 : _lm_dofs.push_back(di[i]);
239 :
240 : // save the corresponding coupled dof indices
241 310 : _global_primary_dofs.push_back(cp_di[i]);
242 310 : if (_dofmap.local_index(cp_di[i]))
243 234 : _primary_dofs.push_back(cp_di[i]);
244 310 : _map_global_lm_primary.insert(std::make_pair(di[i], cp_di[i]));
245 : }
246 : }
247 :
248 : // check if we endup with none dof to condense
249 64 : if (_global_lm_dofs.empty())
250 : {
251 0 : _need_condense = false;
252 0 : _console << std::endl
253 : << "The variable(s) provided do not have a saddle-point character at this step. VCP "
254 0 : "will continue without condensing the dofs."
255 0 : << std::endl
256 0 : << std::endl;
257 : }
258 : else
259 64 : _need_condense = true;
260 :
261 64 : std::sort(_global_lm_dofs.begin(), _global_lm_dofs.end());
262 64 : std::sort(_lm_dofs.begin(), _lm_dofs.end());
263 :
264 64 : std::sort(_global_primary_dofs.begin(), _global_primary_dofs.end());
265 64 : std::sort(_primary_dofs.begin(), _primary_dofs.end());
266 :
267 374 : for (const auto & i : index_range(_global_lm_dofs))
268 : {
269 310 : auto it = _map_global_lm_primary.find(_global_lm_dofs[i]);
270 : mooseAssert(it != _map_global_lm_primary.end(), "Index does not exist in the map.");
271 310 : _map_global_primary_order.insert(std::make_pair(it->second, i));
272 : }
273 64 : }
274 :
275 : void
276 64 : VariableCondensationPreconditioner::getDofColRow()
277 : {
278 : // clean the containers if we want to update the dofs
279 64 : _global_rows.clear();
280 64 : _rows.clear();
281 64 : _global_cols.clear();
282 64 : _cols.clear();
283 64 : _global_rows_to_idx.clear();
284 64 : _rows_to_idx.clear();
285 64 : _global_cols_to_idx.clear();
286 64 : _cols_to_idx.clear();
287 : // row: all without primary variable dofs
288 : // col: all without lm variable dofs
289 2002 : for (dof_id_type i = 0; i < _dofmap.n_dofs(); ++i)
290 : {
291 1938 : if (std::find(_global_primary_dofs.begin(), _global_primary_dofs.end(), i) !=
292 3876 : _global_primary_dofs.end())
293 310 : continue;
294 :
295 1628 : _global_rows.push_back(i);
296 1628 : _global_rows_to_idx.insert(std::make_pair(i, _global_rows.size() - 1));
297 1628 : if (_dofmap.local_index(i))
298 : {
299 1220 : _rows.push_back(i);
300 1220 : _rows_to_idx.insert(std::make_pair(i, _global_rows_to_idx[i]));
301 : }
302 :
303 : // ensure the lm and primary correspondance, so that the condensed Jacobian has non-zero
304 : // diagonal if the dof corresponds to the lm variable, then find the corresponding primary
305 : // variable dof and add to _global_cols
306 1628 : if (_map_global_lm_primary.find(i) != _map_global_lm_primary.end())
307 : {
308 310 : auto primary_idx = _map_global_lm_primary[i];
309 310 : _global_cols.push_back(primary_idx);
310 310 : _global_cols_to_idx.insert(std::make_pair(primary_idx, _global_cols.size() - 1));
311 :
312 310 : if (_dofmap.local_index(primary_idx))
313 : {
314 234 : _cols.push_back(primary_idx);
315 234 : _cols_to_idx.insert(std::make_pair(primary_idx, _global_cols_to_idx[primary_idx]));
316 : }
317 : }
318 : else // if the dof does not correspond to the lm nor primary variable, just add to _global_cols
319 : {
320 1318 : _global_cols.push_back(i);
321 1318 : _global_cols_to_idx.insert(std::make_pair(i, _global_cols.size() - 1));
322 :
323 1318 : if (_dofmap.local_index(i))
324 : {
325 986 : _cols.push_back(i);
326 986 : _cols_to_idx.insert(std::make_pair(i, _global_cols_to_idx[i]));
327 : }
328 : }
329 : }
330 64 : }
331 :
332 : void
333 72 : VariableCondensationPreconditioner::init()
334 : {
335 72 : TIME_SECTION(_init_timer);
336 :
337 72 : if (!_preconditioner)
338 : _preconditioner =
339 36 : Preconditioner<Number>::build_preconditioner(MoosePreconditioner::_communicator);
340 :
341 72 : _is_initialized = true;
342 72 : }
343 :
344 : void
345 64 : VariableCondensationPreconditioner::condenseSystem()
346 : {
347 : // extract _M from the original matrix
348 64 : _matrix->create_submatrix(*_M, _rows, _lm_dofs);
349 :
350 : // get the row associated with the coupled primary variable
351 64 : _K->init(_global_primary_dofs.size(), _global_cols.size(), _primary_dofs.size(), _cols.size());
352 : // Note: enabling nonzero allocation may be expensive. Improved memeory pre-allocation will be
353 : // investigated in the future
354 64 : LibmeshPetscCallA(this->MoosePreconditioner::comm().get(),
355 : MatSetOption(_K->mat(), MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
356 : // here the _global_cols may not be sorted
357 64 : _matrix->create_submatrix_nosort(*_K, _global_primary_dofs, _global_cols);
358 :
359 64 : _matrix->create_submatrix(*_D, _primary_dofs, _lm_dofs);
360 :
361 : // clean dinv
362 64 : if (_dinv)
363 : {
364 28 : LibmeshPetscCallA(this->MoosePreconditioner::comm().get(), MatDestroy(&_dinv));
365 28 : _dinv = nullptr;
366 : }
367 :
368 : // Compute inverse of D
369 64 : if (_is_lm_coupling_diagonal)
370 : // when _D is strictly diagonal, we only need to compute the reciprocal number of the diagonal
371 : // entries
372 50 : computeDInverseDiag(_dinv);
373 : else
374 : // for general cases when _D is not necessarily strict diagonal, we compute the inverse of _D
375 : // using LU
376 14 : computeDInverse(_dinv);
377 :
378 : Mat MdinvK;
379 : // calculate MdinvK
380 64 : LibmeshPetscCallA(
381 : this->MoosePreconditioner::comm().get(),
382 : MatMatMatMult(_M->mat(), _dinv, _K->mat(), MAT_INITIAL_MATRIX, PETSC_DEFAULT, &MdinvK));
383 64 : PetscMatrix<Number> MDinv_K(MdinvK, MoosePreconditioner::_communicator);
384 :
385 : // Preallocate memory for _J_condensed
386 : // memory info is obtained from _matrix and MDinv_K
387 : // indices are from _global_rows, _global_cols
388 64 : auto pc_original_mat = cast_ptr<PetscMatrix<Number> *>(_matrix);
389 64 : preallocateCondensedJacobian(
390 64 : *_J_condensed, *pc_original_mat, _rows, _cols, _global_rows, _global_cols, MDinv_K);
391 :
392 : // Extract unchanged parts from _matrix and add changed parts (MDinv_K) to _J_condensed
393 64 : computeCondensedJacobian(*_J_condensed, *pc_original_mat, _global_rows, MDinv_K);
394 :
395 : // Destroy MdinvK here otherwise we will have memory leak
396 64 : LibmeshPetscCallA(this->MoosePreconditioner::comm().get(), MatDestroy(&MdinvK));
397 64 : }
398 :
399 : void
400 64 : VariableCondensationPreconditioner::computeCondensedJacobian(PetscMatrix<Number> & condensed_mat,
401 : PetscMatrix<Number> & original_mat,
402 : const std::vector<dof_id_type> & grows,
403 : PetscMatrix<Number> & block_mat)
404 : {
405 : // obtain entries from the original matrix
406 64 : PetscInt pc_ncols = 0, block_ncols = 0;
407 : const PetscInt *pc_cols, *block_cols;
408 : const PetscScalar *pc_vals, *block_vals;
409 :
410 : // containers for the data
411 64 : std::vector<PetscInt> sub_cols;
412 64 : std::vector<PetscScalar> sub_vals;
413 :
414 1692 : for (const auto & i : index_range(grows))
415 : {
416 1628 : PetscInt sub_rid[] = {static_cast<PetscInt>(i)};
417 1628 : PetscInt rid = grows[i];
418 1628 : if (grows[i] >= original_mat.row_start() && grows[i] < original_mat.row_stop())
419 : {
420 : // get one row of data from the original matrix
421 1220 : LibmeshPetscCallA(this->MoosePreconditioner::comm().get(),
422 : MatGetRow(original_mat.mat(), rid, &pc_ncols, &pc_cols, &pc_vals));
423 : // get corresponding row of data from the block matrix
424 1220 : LibmeshPetscCallA(this->MoosePreconditioner::comm().get(),
425 : MatGetRow(block_mat.mat(), i, &block_ncols, &block_cols, &block_vals));
426 : // extract data from certain cols, subtract the value from the block mat, and save the indices
427 : // and entries sub_cols and sub_vals
428 : // First, save the submatrix col index and value as a map
429 1220 : std::map<PetscInt, PetscScalar> pc_col_map;
430 15886 : for (PetscInt pc_idx = 0; pc_idx < pc_ncols; pc_idx++)
431 : {
432 : // save only if the col exists in the condensed matrix
433 14666 : if (_global_cols_to_idx.find(pc_cols[pc_idx]) != _global_cols_to_idx.end())
434 11338 : pc_col_map.insert(std::make_pair(_global_cols_to_idx[pc_cols[pc_idx]], pc_vals[pc_idx]));
435 : }
436 : // Second, check the block cols and calculate new entries for the condensed system
437 14328 : for (PetscInt block_idx = 0; block_idx < block_ncols; block_idx++)
438 : {
439 13108 : PetscInt block_col = block_cols[block_idx];
440 13108 : PetscScalar block_val = block_vals[block_idx];
441 : // if the block mat has nonzero at the same column, subtract value
442 : // otherwise, create a new key and save the negative value from the block matrix
443 13108 : if (pc_col_map.find(block_col) != pc_col_map.end())
444 4510 : pc_col_map[block_col] -= block_val;
445 : else
446 8598 : pc_col_map[block_col] = -block_val;
447 : }
448 :
449 : // Third, save keys in the sub_cols and values in the sub_vals
450 1220 : for (std::map<PetscInt, PetscScalar>::iterator it = pc_col_map.begin();
451 21156 : it != pc_col_map.end();
452 19936 : ++it)
453 : {
454 19936 : sub_cols.push_back(it->first);
455 19936 : sub_vals.push_back(it->second);
456 : }
457 :
458 : // Then, set values
459 1220 : LibmeshPetscCallA(this->MoosePreconditioner::comm().get(),
460 : MatSetValues(condensed_mat.mat(),
461 : 1,
462 : sub_rid,
463 : sub_vals.size(),
464 : sub_cols.data(),
465 : sub_vals.data(),
466 : INSERT_VALUES));
467 1220 : LibmeshPetscCallA(this->MoosePreconditioner::comm().get(),
468 : MatRestoreRow(original_mat.mat(), rid, &pc_ncols, &pc_cols, &pc_vals));
469 1220 : LibmeshPetscCallA(this->MoosePreconditioner::comm().get(),
470 : MatRestoreRow(block_mat.mat(), i, &block_ncols, &block_cols, &block_vals));
471 : // clear data for this row
472 1220 : sub_cols.clear();
473 1220 : sub_vals.clear();
474 1220 : }
475 : }
476 64 : condensed_mat.close();
477 64 : }
478 :
479 : void
480 64 : VariableCondensationPreconditioner::preallocateCondensedJacobian(
481 : PetscMatrix<Number> & condensed_mat,
482 : PetscMatrix<Number> & original_mat,
483 : const std::vector<dof_id_type> & rows,
484 : const std::vector<dof_id_type> & cols,
485 : const std::vector<dof_id_type> & grows,
486 : const std::vector<dof_id_type> & gcols,
487 : PetscMatrix<Number> & block_mat)
488 : {
489 : // quantities from the original matrix and the block matrix
490 64 : PetscInt ncols = 0, block_ncols = 0;
491 : const PetscInt * col_vals;
492 : const PetscInt * block_col_vals;
493 : const PetscScalar * vals;
494 : const PetscScalar * block_vals;
495 :
496 64 : std::vector<PetscInt> block_cols_to_org; // stores the nonzero column indices of the block
497 : // matrix w.r.t original matrix
498 : std::vector<PetscInt>
499 64 : merged_cols; // stores the nonzero column indices estimate of the condensed matrix
500 :
501 : // number of nonzeros in each row of the DIAGONAL and OFF-DIAGONAL portion of the local
502 : // condensed matrix
503 64 : std::vector<dof_id_type> n_nz, n_oz;
504 :
505 : // Get number of nonzeros from original_mat and block_mat for each row
506 1284 : for (const auto & row_id : _rows)
507 : {
508 : // get number of non-zeros in the original matrix
509 1220 : LibmeshPetscCallA(this->MoosePreconditioner::comm().get(),
510 : MatGetRow(original_mat.mat(), row_id, &ncols, &col_vals, &vals));
511 :
512 : // get number of non-zeros in the block matrix
513 : dof_id_type block_row_id; // row id in the block matrix
514 :
515 1220 : if (_global_rows_to_idx.find(row_id) != _global_rows_to_idx.end())
516 1220 : block_row_id = _global_rows_to_idx[row_id];
517 : else
518 0 : mooseError("DoF ", row_id, " does not exist in the rows of condensed_mat");
519 :
520 1220 : LibmeshPetscCallA(
521 : this->MoosePreconditioner::comm().get(),
522 : MatGetRow(block_mat.mat(), block_row_id, &block_ncols, &block_col_vals, &block_vals));
523 :
524 : // make sure the block index is transformed in terms of the original mat
525 1220 : block_cols_to_org.clear();
526 14328 : for (PetscInt i = 0; i < block_ncols; ++i)
527 : {
528 13108 : auto idx = gcols[block_col_vals[i]];
529 13108 : block_cols_to_org.push_back(idx);
530 : }
531 :
532 : // Now store nonzero column indices for the condensed Jacobian
533 : // merge `col_vals` and `block_cols_to_org` and save the common indices in `merged_cols`.
534 1220 : mergeArrays(col_vals, block_cols_to_org.data(), ncols, block_ncols, merged_cols);
535 :
536 1220 : LibmeshPetscCallA(
537 : this->MoosePreconditioner::comm().get(),
538 : MatRestoreRow(block_mat.mat(), block_row_id, &block_ncols, &block_col_vals, &block_vals));
539 :
540 1220 : LibmeshPetscCallA(this->MoosePreconditioner::comm().get(),
541 : MatRestoreRow(original_mat.mat(), row_id, &ncols, &col_vals, &vals));
542 :
543 : // Count the nnz for DIAGONAL and OFF-DIAGONAL parts
544 1220 : PetscInt row_n_nz = 0, row_n_oz = 0;
545 24484 : for (const auto & merged_col : merged_cols)
546 : {
547 : // find corresponding index in the block mat and skip the cols that do not exist in the
548 : // condensed system
549 23264 : if (_global_cols_to_idx.find(merged_col) == _global_cols_to_idx.end())
550 3328 : continue;
551 :
552 19936 : dof_id_type col_idx = _global_cols_to_idx[merged_col];
553 : // find the corresponding row index
554 19936 : dof_id_type row_idx = grows[col_idx];
555 : // check whether the index is local;
556 : // yes - DIAGONAL, no - OFF-DIAGONAL
557 19936 : if (_rows_to_idx.find(row_idx) != _rows_to_idx.end())
558 17312 : row_n_nz++;
559 : else
560 2624 : row_n_oz++;
561 : }
562 :
563 1220 : n_nz.push_back(cast_int<dof_id_type>(row_n_nz));
564 1220 : n_oz.push_back(cast_int<dof_id_type>(row_n_oz));
565 : }
566 : // Then initialize and allocate memory for the condensed system matrix
567 64 : condensed_mat.init(grows.size(), gcols.size(), rows.size(), cols.size(), n_nz, n_oz);
568 64 : }
569 :
570 : void
571 1220 : VariableCondensationPreconditioner::mergeArrays(const PetscInt * a,
572 : const PetscInt * b,
573 : const PetscInt & na,
574 : const PetscInt & nb,
575 : std::vector<PetscInt> & c)
576 : {
577 1220 : c.clear();
578 :
579 : // use map to store unique elements.
580 1220 : std::map<PetscInt, bool> mp;
581 :
582 : // Inserting values to a map.
583 15886 : for (const auto & i : make_range(na))
584 14666 : mp[a[i]] = true;
585 :
586 14328 : for (const auto & i : make_range(nb))
587 13108 : mp[b[i]] = true;
588 :
589 : // Save the merged values to c, if only the value also exist in gcols
590 24484 : for (const auto & i : mp)
591 23264 : c.push_back(i.first);
592 1220 : }
593 :
594 : void
595 64 : VariableCondensationPreconditioner::setup()
596 : {
597 64 : if (_adaptive_condensation)
598 64 : findZeroDiagonals(*_matrix, _zero_rows);
599 :
600 : // save dofs that are to be condensed out
601 64 : getDofToCondense();
602 :
603 : // solve the condensed system only when needed, otherwise solve the original system
604 64 : if (_need_condense)
605 : {
606 : // get condensed dofs for rows and cols
607 64 : getDofColRow();
608 :
609 64 : condenseSystem();
610 :
611 : // make sure diagonal entries are not empty
612 1284 : for (const auto & i : make_range(_J_condensed->row_start(), _J_condensed->row_stop()))
613 1220 : _J_condensed->add(i, i, 0.0);
614 64 : _J_condensed->close();
615 :
616 64 : _preconditioner->set_matrix(*_J_condensed);
617 : }
618 : else
619 0 : _preconditioner->set_matrix(*_matrix);
620 :
621 64 : _preconditioner->set_type(_pre_type);
622 64 : _preconditioner->init();
623 64 : }
624 :
625 : void
626 378 : VariableCondensationPreconditioner::apply(const NumericVector<Number> & y,
627 : NumericVector<Number> & x)
628 : {
629 378 : TIME_SECTION(_apply_timer);
630 :
631 378 : if (_need_condense)
632 : {
633 378 : getCondensedXY(y, x);
634 :
635 378 : _preconditioner->apply(*_y_hat, *_x_hat);
636 :
637 378 : computeCondensedVariables();
638 :
639 378 : getFullSolution(y, x);
640 : }
641 : else
642 : {
643 0 : _preconditioner->apply(y, x);
644 : }
645 378 : }
646 :
647 : void
648 378 : VariableCondensationPreconditioner::getCondensedXY(const NumericVector<Number> & y,
649 : NumericVector<Number> & x)
650 : {
651 : Mat mdinv;
652 : // calculate mdinv
653 378 : LibmeshPetscCallA(this->MoosePreconditioner::comm().get(),
654 : MatMatMult(_M->mat(), _dinv, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &mdinv));
655 378 : PetscMatrix<Number> MDinv(mdinv, MoosePreconditioner::_communicator);
656 :
657 378 : _x_hat->init(_J_condensed->n(), _J_condensed->local_n(), false, PARALLEL);
658 378 : _y_hat->init(_J_condensed->m(), _J_condensed->local_m(), false, PARALLEL);
659 :
660 378 : x.create_subvector(*_x_hat, _global_cols);
661 378 : y.create_subvector(*_y_hat, _global_rows);
662 :
663 378 : _primary_rhs_vec->init(MDinv.n(), MDinv.local_n(), false, PARALLEL);
664 :
665 : std::unique_ptr<NumericVector<Number>> mdinv_primary_rhs(
666 378 : NumericVector<Number>::build(MoosePreconditioner::_communicator));
667 378 : mdinv_primary_rhs->init(MDinv.m(), MDinv.local_m(), false, PARALLEL);
668 :
669 : // get _primary_rhs_vec from the original y
670 378 : y.create_subvector(*_primary_rhs_vec, _global_primary_dofs);
671 :
672 378 : MDinv.vector_mult(*mdinv_primary_rhs, *_primary_rhs_vec);
673 378 : mdinv_primary_rhs->close();
674 :
675 378 : (*_y_hat) -= (*mdinv_primary_rhs);
676 :
677 378 : _y_hat->close();
678 378 : _x_hat->close();
679 :
680 378 : LibmeshPetscCallA(this->MoosePreconditioner::comm().get(), MatDestroy(&mdinv));
681 378 : }
682 :
683 : void
684 378 : VariableCondensationPreconditioner::computeCondensedVariables()
685 : {
686 378 : _lm_sol_vec->clear();
687 :
688 378 : PetscMatrix<Number> Dinv(_dinv, MoosePreconditioner::_communicator);
689 :
690 378 : _lm_sol_vec->init(_D->m(), _D->local_m(), false, PARALLEL);
691 :
692 : std::unique_ptr<NumericVector<Number>> K_xhat(
693 378 : NumericVector<Number>::build(MoosePreconditioner::_communicator));
694 378 : K_xhat->init(_K->m(), _K->local_m(), false, PARALLEL);
695 378 : _K->vector_mult(*K_xhat, *_x_hat);
696 378 : K_xhat->close();
697 :
698 378 : (*_primary_rhs_vec) -= (*K_xhat);
699 378 : _primary_rhs_vec->close();
700 378 : Dinv.vector_mult(*_lm_sol_vec, *_primary_rhs_vec);
701 378 : _lm_sol_vec->close();
702 378 : }
703 :
704 : void
705 378 : VariableCondensationPreconditioner::getFullSolution(const NumericVector<Number> & /*y*/,
706 : NumericVector<Number> & x)
707 : {
708 378 : std::vector<dof_id_type> dof_indices;
709 378 : std::vector<Number> vals;
710 :
711 : // save values and indices from _x_hat and _lm_sol_vec
712 8222 : for (const auto & i : make_range(_x_hat->first_local_index(), _x_hat->last_local_index()))
713 : {
714 7844 : dof_indices.push_back(_global_cols[i]);
715 7844 : vals.push_back((*_x_hat)(i));
716 : }
717 :
718 378 : for (const auto & i :
719 2224 : make_range(_lm_sol_vec->first_local_index(), _lm_sol_vec->last_local_index()))
720 : {
721 1468 : dof_indices.push_back(_global_lm_dofs[i]);
722 1468 : vals.push_back((*_lm_sol_vec)(i));
723 : }
724 :
725 378 : x.insert(vals.data(), dof_indices);
726 378 : x.close();
727 378 : }
728 :
729 : void
730 64 : VariableCondensationPreconditioner::findZeroDiagonals(SparseMatrix<Number> & mat,
731 : std::vector<dof_id_type> & indices)
732 : {
733 64 : indices.clear();
734 : IS zerodiags, zerodiags_all;
735 : const PetscInt * petsc_idx;
736 : PetscInt nrows;
737 : // make sure we have a PETSc matrix
738 64 : auto * const petsc_mat = cast_ptr<PetscMatrix<Number> *>(&mat);
739 64 : LibmeshPetscCallA(this->MoosePreconditioner::comm().get(),
740 : MatFindZeroDiagonals(petsc_mat->mat(), &zerodiags));
741 : // synchronize all indices
742 64 : LibmeshPetscCallA(this->MoosePreconditioner::comm().get(),
743 : ISAllGather(zerodiags, &zerodiags_all));
744 64 : LibmeshPetscCallA(this->MoosePreconditioner::comm().get(),
745 : ISGetIndices(zerodiags_all, &petsc_idx));
746 64 : LibmeshPetscCallA(this->MoosePreconditioner::comm().get(), ISGetSize(zerodiags_all, &nrows));
747 :
748 374 : for (PetscInt i = 0; i < nrows; ++i)
749 310 : indices.push_back(petsc_idx[i]);
750 :
751 64 : LibmeshPetscCallA(this->MoosePreconditioner::comm().get(),
752 : ISRestoreIndices(zerodiags_all, &petsc_idx));
753 64 : LibmeshPetscCallA(this->MoosePreconditioner::comm().get(), ISDestroy(&zerodiags));
754 64 : LibmeshPetscCallA(this->MoosePreconditioner::comm().get(), ISDestroy(&zerodiags_all));
755 64 : }
756 :
757 : void
758 40 : VariableCondensationPreconditioner::clear()
759 : {
760 40 : if (_dinv != nullptr)
761 36 : LibmeshPetscCallA(this->MoosePreconditioner::comm().get(), MatDestroy(&_dinv));
762 40 : }
763 :
764 : void
765 14 : VariableCondensationPreconditioner::computeDInverse(Mat & dinv)
766 : {
767 : Mat F, I, dinv_dense;
768 : IS perm, iperm;
769 : MatFactorInfo info;
770 :
771 14 : LibmeshPetscCallA(
772 : this->MoosePreconditioner::comm().get(),
773 : MatCreateDense(
774 : PETSC_COMM_WORLD, _D->local_n(), _D->local_m(), _D->n(), _D->m(), NULL, &dinv_dense));
775 :
776 : // Create an identity matrix as the right-hand-side
777 14 : LibmeshPetscCallA(
778 : this->MoosePreconditioner::comm().get(),
779 : MatCreateDense(PETSC_COMM_WORLD, _D->local_m(), _D->local_m(), _D->m(), _D->m(), NULL, &I));
780 :
781 70 : for (unsigned int i = 0; i < _D->m(); ++i)
782 56 : LibmeshPetscCallA(this->MoosePreconditioner::comm().get(),
783 : MatSetValue(I, i, i, 1.0, INSERT_VALUES));
784 :
785 14 : LibmeshPetscCallA(this->MoosePreconditioner::comm().get(),
786 : MatAssemblyBegin(I, MAT_FINAL_ASSEMBLY));
787 14 : LibmeshPetscCallA(this->MoosePreconditioner::comm().get(), MatAssemblyEnd(I, MAT_FINAL_ASSEMBLY));
788 :
789 : // Factorize D
790 14 : LibmeshPetscCallA(this->MoosePreconditioner::comm().get(),
791 : MatGetOrdering(_D->mat(), MATORDERINGND, &perm, &iperm));
792 :
793 14 : LibmeshPetscCallA(this->MoosePreconditioner::comm().get(), MatFactorInfoInitialize(&info));
794 :
795 14 : LibmeshPetscCallA(this->MoosePreconditioner::comm().get(),
796 : MatGetFactor(_D->mat(), MATSOLVERSUPERLU_DIST, MAT_FACTOR_LU, &F));
797 :
798 14 : LibmeshPetscCallA(this->MoosePreconditioner::comm().get(),
799 : MatLUFactorSymbolic(F, _D->mat(), perm, iperm, &info));
800 :
801 14 : LibmeshPetscCallA(this->MoosePreconditioner::comm().get(),
802 : MatLUFactorNumeric(F, _D->mat(), &info));
803 :
804 : // Solve for Dinv
805 14 : LibmeshPetscCallA(this->MoosePreconditioner::comm().get(), MatMatSolve(F, I, dinv_dense));
806 :
807 14 : LibmeshPetscCallA(this->MoosePreconditioner::comm().get(),
808 : MatAssemblyBegin(dinv_dense, MAT_FINAL_ASSEMBLY));
809 14 : LibmeshPetscCallA(this->MoosePreconditioner::comm().get(),
810 : MatAssemblyEnd(dinv_dense, MAT_FINAL_ASSEMBLY));
811 :
812 : // copy value to dinv
813 14 : LibmeshPetscCallA(this->MoosePreconditioner::comm().get(),
814 : MatConvert(dinv_dense, MATAIJ, MAT_INITIAL_MATRIX, &dinv));
815 :
816 14 : LibmeshPetscCallA(this->MoosePreconditioner::comm().get(), MatDestroy(&dinv_dense));
817 :
818 14 : LibmeshPetscCallA(this->MoosePreconditioner::comm().get(), MatDestroy(&I));
819 14 : LibmeshPetscCallA(this->MoosePreconditioner::comm().get(), MatDestroy(&F));
820 14 : LibmeshPetscCallA(this->MoosePreconditioner::comm().get(), ISDestroy(&perm));
821 14 : LibmeshPetscCallA(this->MoosePreconditioner::comm().get(), ISDestroy(&iperm));
822 14 : }
823 :
824 : void
825 50 : VariableCondensationPreconditioner::computeDInverseDiag(Mat & dinv)
826 : {
827 50 : auto diag_D = NumericVector<Number>::build(MoosePreconditioner::_communicator);
828 : // Initialize dinv
829 50 : LibmeshPetscCallA(this->MoosePreconditioner::comm().get(),
830 : MatCreateAIJ(PETSC_COMM_WORLD,
831 : _D->local_n(),
832 : _D->local_m(),
833 : _D->n(),
834 : _D->m(),
835 : 1,
836 : NULL,
837 : 0,
838 : NULL,
839 : &dinv));
840 : // Allocate storage
841 50 : diag_D->init(_D->m(), _D->local_m(), false, PARALLEL);
842 : // Fill entries
843 244 : for (const auto & i : make_range(_D->row_start(), _D->row_stop()))
844 : {
845 194 : auto it = _map_global_primary_order.find(_global_primary_dofs[i]);
846 : mooseAssert(it != _map_global_primary_order.end(), "Index does not exist in the map.");
847 194 : diag_D->set(it->second, (*_D)(i, it->second));
848 : }
849 :
850 244 : for (const auto & i : make_range(_D->row_start(), _D->row_stop()))
851 : {
852 194 : if (MooseUtils::absoluteFuzzyEqual((*diag_D)(i), 0.0))
853 0 : mooseError("Trying to compute reciprocal of 0.");
854 194 : LibmeshPetscCallA(this->MoosePreconditioner::comm().get(),
855 : MatSetValue(dinv,
856 : i,
857 : _map_global_primary_order.at(_global_primary_dofs[i]),
858 : 1.0 / (*diag_D)(i),
859 : INSERT_VALUES));
860 : }
861 :
862 50 : LibmeshPetscCallA(this->MoosePreconditioner::comm().get(),
863 : MatAssemblyBegin(dinv, MAT_FINAL_ASSEMBLY));
864 50 : LibmeshPetscCallA(this->MoosePreconditioner::comm().get(),
865 : MatAssemblyEnd(dinv, MAT_FINAL_ASSEMBLY));
866 50 : }
|
