Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 : #include "libmesh/enum_parallel_type.h"
19 : #include "libmesh/static_condensation.h"
20 :
21 : #if defined(LIBMESH_HAVE_EIGEN) && defined(LIBMESH_HAVE_PETSC)
22 :
23 : #include "libmesh/mesh_base.h"
24 : #include "libmesh/dof_map.h"
25 : #include "libmesh/elem.h"
26 : #include "libmesh/int_range.h"
27 : #include "libmesh/numeric_vector.h"
28 : #include "libmesh/linear_solver.h"
29 : #include "libmesh/static_condensation_preconditioner.h"
30 : #include "libmesh/system.h"
31 : #include "libmesh/petsc_matrix.h"
32 : #include "libmesh/equation_systems.h"
33 : #include "libmesh/static_condensation_dof_map.h"
34 : #include "timpi/parallel_sync.h"
35 : #include <unordered_set>
36 :
37 : namespace libMesh
38 : {
39 630 : StaticCondensation::StaticCondensation(const MeshBase & mesh,
40 : System & system,
41 : const DofMap & full_dof_map,
42 630 : StaticCondensationDofMap & reduced_dof_map)
43 : : PetscMatrixShellMatrix<Number>(full_dof_map.comm()),
44 594 : _mesh(mesh),
45 594 : _system(system),
46 594 : _full_dof_map(full_dof_map),
47 594 : _reduced_dof_map(reduced_dof_map),
48 594 : _current_elem_id(DofObject::invalid_id),
49 594 : _sc_is_initialized(false),
50 594 : _have_cached_values(false),
51 594 : _parallel_type(INVALID_PARALLELIZATION),
52 630 : _uncondensed_dofs_only(false)
53 : {
54 612 : _size_one_mat.resize(1, 1);
55 1224 : _scp = std::make_unique<StaticCondensationPreconditioner>(*this);
56 630 : }
57 :
58 1188 : StaticCondensation::~StaticCondensation() = default;
59 :
60 0 : SparseMatrix<Number> & StaticCondensation::operator=(const SparseMatrix<Number> &)
61 : {
62 0 : libmesh_not_implemented();
63 : }
64 :
65 0 : std::unique_ptr<SparseMatrix<Number>> StaticCondensation::zero_clone() const
66 : {
67 0 : libmesh_not_implemented();
68 : }
69 :
70 0 : std::unique_ptr<SparseMatrix<Number>> StaticCondensation::clone() const
71 : {
72 0 : libmesh_not_implemented();
73 : }
74 :
75 4130 : void StaticCondensation::clear() noexcept
76 : {
77 4130 : PetscMatrixShellMatrix<Number>::clear();
78 :
79 118 : _elem_to_matrix_data.clear();
80 118 : _reduced_sys_mat.reset();
81 118 : _reduced_sol.reset();
82 118 : _reduced_rhs.reset();
83 118 : _reduced_solver.reset();
84 4130 : _current_elem_id = DofObject::invalid_id;
85 4130 : _have_cached_values = false;
86 4130 : _sc_is_initialized = false;
87 4130 : }
88 :
89 0 : void StaticCondensation::init(const numeric_index_type m,
90 : const numeric_index_type n,
91 : const numeric_index_type m_l,
92 : const numeric_index_type n_l,
93 : const numeric_index_type nnz,
94 : const numeric_index_type noz,
95 : const numeric_index_type blocksize)
96 : {
97 0 : if (!this->initialized())
98 : {
99 0 : PetscMatrixShellMatrix<Number>::init(m, n, m_l, n_l, nnz, noz, blocksize);
100 0 : _parallel_type = ((m == m_l) && (this->n_processors() > 1)) ? SERIAL : PARALLEL;
101 0 : this->init();
102 : }
103 0 : }
104 :
105 4130 : void StaticCondensation::init(const ParallelType type)
106 : {
107 4130 : if (!this->initialized())
108 : {
109 4130 : PetscMatrixShellMatrix<Number>::init(type);
110 4130 : _parallel_type = type;
111 4130 : this->init();
112 : }
113 4130 : }
114 :
115 33084 : bool StaticCondensation::initialized() const
116 : {
117 33084 : return PetscMatrixShellMatrix<Number>::initialized() && _sc_is_initialized;
118 : }
119 :
120 15750 : void StaticCondensation::init()
121 : {
122 15750 : if (this->initialized())
123 332 : return;
124 :
125 118 : libmesh_assert(_reduced_dof_map.initialized());
126 :
127 : // This API is public, so it can be called without going through the other init overloads which
128 : // would give us an indication of what kind of parallel type the user wants. If we've gotten no
129 : // indication so far, we will default to parallel
130 4130 : if (_parallel_type == INVALID_PARALLELIZATION)
131 0 : _parallel_type = PARALLEL;
132 118 : libmesh_assert(_parallel_type != SERIAL);
133 :
134 190059 : for (const auto & [elem_id, dof_data] : _reduced_dof_map._elem_to_dof_data)
135 : {
136 15743 : auto & matrix_data = _elem_to_matrix_data[elem_id];
137 :
138 15743 : const auto condensed_dof_size = dof_data.condensed_global_to_local_map.size();
139 15743 : const auto uncondensed_dof_size = dof_data.uncondensed_global_to_local_map.size();
140 :
141 185929 : matrix_data.Acc.setZero(condensed_dof_size, condensed_dof_size);
142 185929 : matrix_data.Acu.setZero(condensed_dof_size, uncondensed_dof_size);
143 185929 : matrix_data.Auc.setZero(uncondensed_dof_size, condensed_dof_size);
144 185929 : matrix_data.Auu.setZero(uncondensed_dof_size, uncondensed_dof_size);
145 : }
146 :
147 : //
148 : // Build the reduced system data
149 : //
150 4130 : const auto n = _reduced_dof_map.n_dofs();
151 : const auto n_local =
152 4130 : (_parallel_type == SERIAL) ? _reduced_dof_map.n_dofs() : _reduced_dof_map.n_local_dofs();
153 8142 : _reduced_solver = LinearSolver<Number>::build(this->comm());
154 4248 : _reduced_solver->init((_system.name() + "_condensed_").c_str());
155 8142 : _reduced_rhs = NumericVector<Number>::build(this->comm());
156 : // Init the RHS vector so we can conveniently get processor row offsets
157 4130 : _reduced_rhs->init(n, n_local);
158 :
159 : // Initialize the reduced system matrix
160 4130 : _sp = _reduced_dof_map._reduced_sp.get();
161 8142 : _reduced_sys_mat = SparseMatrix<Number>::build(this->comm());
162 4130 : if (auto * const petsc_mat = dynamic_cast<PetscMatrix<Number> *>(_reduced_sys_mat.get()))
163 : {
164 : // Optimization for PETSc. This is critical for problems in which there are SCALAR dofs that
165 : // introduce dense rows to avoid allocating a dense matrix
166 4130 : petsc_mat->init(
167 4130 : n, n, n_local, n_local, _reduced_dof_map._reduced_nnz, _reduced_dof_map._reduced_noz);
168 : }
169 : else
170 : {
171 0 : const auto & nnz = _sp->get_n_nz();
172 0 : const auto & noz = _sp->get_n_oz();
173 0 : const auto nz = nnz.empty() ? dof_id_type(0) : *std::max_element(nnz.begin(), nnz.end());
174 0 : const auto oz = noz.empty() ? dof_id_type(0) : *std::max_element(noz.begin(), noz.end());
175 0 : _reduced_sys_mat->init(n, n, n_local, n_local, nz, oz);
176 : }
177 :
178 : // Build ghosted full solution vector. Note that this is, in general, *not equal* to the system
179 : // solution, e.g. this may correspond to the solution for the Newton *update*
180 8142 : _ghosted_full_sol = _system.current_local_solution->clone();
181 :
182 4130 : _sc_is_initialized = true;
183 : }
184 :
185 10780 : void StaticCondensation::close()
186 : {
187 10780 : _communicator.max(_have_cached_values);
188 10780 : if (!_have_cached_values)
189 : {
190 2450 : bool closed = _reduced_sys_mat->closed();
191 : // closed is not collective
192 2450 : _communicator.min(closed);
193 2450 : if (!closed)
194 0 : _reduced_sys_mat->close();
195 70 : return;
196 : }
197 :
198 8806 : DenseMatrix<Number> shim;
199 238 : std::vector<dof_id_type> reduced_space_indices;
200 358659 : for (auto & [elem_id, matrix_data] : _elem_to_matrix_data)
201 : {
202 350329 : const auto & dof_data = libmesh_map_find(_reduced_dof_map._elem_to_dof_data, elem_id);
203 29585 : reduced_space_indices.clear();
204 :
205 : // The result matrix is either a Schur complement or it's simply the result of summing element
206 : // matrices of the uncondensed degrees of freedom
207 59170 : EigenMatrix result = matrix_data.Auu;
208 350329 : if (!_uncondensed_dofs_only)
209 : {
210 350329 : matrix_data.AccFactor = matrix_data.Acc.partialPivLu();
211 350329 : result -= matrix_data.Auc * matrix_data.AccFactor.solve(matrix_data.Acu);
212 : }
213 379914 : shim.resize(result.rows(), result.cols());
214 4741532 : for (const auto i : make_range(result.rows()))
215 61779900 : for (const auto j : make_range(result.cols()))
216 62281893 : shim(i, j) = result(i, j);
217 706994 : for (const auto & var_reduced_space_indices : dof_data.reduced_space_indices)
218 326504 : reduced_space_indices.insert(reduced_space_indices.end(),
219 : var_reduced_space_indices.begin(),
220 120644 : var_reduced_space_indices.end());
221 350329 : _reduced_sys_mat->add_matrix(shim, reduced_space_indices);
222 : }
223 :
224 8330 : _reduced_sys_mat->close();
225 :
226 8330 : _have_cached_values = false;
227 7854 : }
228 :
229 118 : bool StaticCondensation::closed() const
230 : {
231 118 : return _reduced_sys_mat->closed() && !_have_cached_values;
232 : }
233 :
234 12460 : void StaticCondensation::zero()
235 : {
236 12460 : _reduced_sys_mat->zero();
237 548718 : for (auto & [elem_id, matrix_data] : _elem_to_matrix_data)
238 : {
239 45328 : libmesh_ignore(elem_id);
240 490930 : matrix_data.Acc.setZero();
241 490930 : matrix_data.Acu.setZero();
242 490930 : matrix_data.Auc.setZero();
243 490930 : matrix_data.Auu.setZero();
244 : }
245 12460 : }
246 :
247 4130 : void StaticCondensation::setup() { libmesh_assert(this->closed()); }
248 :
249 0 : numeric_index_type StaticCondensation::m() const { return _full_dof_map.n_dofs(); }
250 :
251 0 : numeric_index_type StaticCondensation::row_start() const { return _full_dof_map.first_dof(); }
252 :
253 0 : numeric_index_type StaticCondensation::row_stop() const { return _full_dof_map.end_dof(); }
254 :
255 0 : void StaticCondensation::set(const numeric_index_type full_i,
256 : const numeric_index_type full_j,
257 : const Number val)
258 : {
259 0 : const auto reduced_i = _reduced_dof_map.get_reduced_from_global_constraint_dof(full_i);
260 0 : const auto reduced_j = _reduced_dof_map.get_reduced_from_global_constraint_dof(full_j);
261 0 : _reduced_sys_mat->set(reduced_i, reduced_j, val);
262 0 : }
263 :
264 350329 : void StaticCondensation::set_current_elem(const Elem & elem)
265 : {
266 29585 : libmesh_assert(!Threads::in_threads || libMesh::n_threads() == 1);
267 350329 : _current_elem_id = elem.id();
268 350329 : }
269 :
270 0 : void StaticCondensation::add(const numeric_index_type i,
271 : const numeric_index_type j,
272 : const Number value)
273 : {
274 0 : _size_one_mat(0, 0) = value;
275 0 : this->add_matrix(_size_one_mat, {i}, {j});
276 0 : }
277 :
278 439033 : void StaticCondensation::add_matrix(const DenseMatrix<Number> & dm,
279 : const std::vector<numeric_index_type> & rows,
280 : const std::vector<numeric_index_type> & cols)
281 : {
282 439033 : if (rows.empty() || cols.empty())
283 0 : return;
284 :
285 37649 : libmesh_assert(_current_elem_id != DofObject::invalid_id);
286 439033 : auto & matrix_data = libmesh_map_find(_elem_to_matrix_data, _current_elem_id);
287 439033 : const auto & dof_data = libmesh_map_find(_reduced_dof_map._elem_to_dof_data, _current_elem_id);
288 : EigenMatrix * mat;
289 :
290 218633634 : auto info_from_index = [&dof_data](const auto global_index) {
291 14899804 : auto index_it = dof_data.condensed_global_to_local_map.find(global_index);
292 14899804 : const bool index_is_condensed = index_it != dof_data.condensed_global_to_local_map.end();
293 173934222 : if (!index_is_condensed)
294 : {
295 12007480 : index_it = dof_data.uncondensed_global_to_local_map.find(global_index);
296 140505916 : if (index_it == dof_data.uncondensed_global_to_local_map.end())
297 0 : libmesh_error_msg("Failed to find the global index "
298 : << global_index
299 : << " in our current element's degree of freedom information. One way "
300 : "this can happen is when using a discontinuous Galerkin method, "
301 : "adding element matrices to both + and - sides of a face");
302 : }
303 : else
304 : // We found the dof in the condensed container. Let's assert that it's not also in the
305 : // uncondensed container
306 2892324 : libmesh_assert(dof_data.uncondensed_global_to_local_map.find(global_index) ==
307 : dof_data.uncondensed_global_to_local_map.end());
308 :
309 203733830 : return std::make_pair(index_is_condensed, index_it->second);
310 401384 : };
311 :
312 6133636 : for (const auto i : make_range(dm.m()))
313 92661714 : for (const auto j : make_range(dm.n()))
314 : {
315 86967111 : const auto global_i = rows[i];
316 86967111 : const auto global_j = cols[j];
317 86967111 : const auto [i_is_condensed, local_i] = info_from_index(global_i);
318 86967111 : const auto [j_is_condensed, local_j] = info_from_index(global_j);
319 86967111 : if (i_is_condensed)
320 : {
321 16764841 : if (j_is_condensed)
322 3806948 : mat = &matrix_data.Acc;
323 : else
324 12957893 : mat = &matrix_data.Acu;
325 : }
326 : else
327 : {
328 70202270 : if (j_is_condensed)
329 12856517 : mat = &matrix_data.Auc;
330 : else
331 57345753 : mat = &matrix_data.Auu;
332 : }
333 94417013 : (*mat)(local_i, local_j) += dm(i, j);
334 : }
335 :
336 439033 : _have_cached_values = true;
337 : }
338 :
339 347513 : void StaticCondensation::add_matrix(const DenseMatrix<Number> & dm,
340 : const std::vector<numeric_index_type> & dof_indices)
341 : {
342 347513 : this->add_matrix(dm, dof_indices, dof_indices);
343 347513 : }
344 :
345 0 : void StaticCondensation::add(const Number, const SparseMatrix<Number> &)
346 : {
347 0 : libmesh_not_implemented();
348 : }
349 :
350 0 : Number StaticCondensation::operator()(const numeric_index_type, const numeric_index_type) const
351 : {
352 0 : libmesh_not_implemented();
353 : }
354 :
355 0 : Real StaticCondensation::l1_norm() const { libmesh_not_implemented(); }
356 :
357 0 : Real StaticCondensation::linfty_norm() const { libmesh_not_implemented(); }
358 :
359 0 : void StaticCondensation::print_personal(std::ostream &) const { libmesh_not_implemented(); }
360 :
361 0 : void StaticCondensation::get_diagonal(NumericVector<Number> &) const { libmesh_not_implemented(); }
362 :
363 0 : void StaticCondensation::get_transpose(SparseMatrix<Number> &) const { libmesh_not_implemented(); }
364 :
365 0 : void StaticCondensation::get_row(numeric_index_type,
366 : std::vector<numeric_index_type> &,
367 : std::vector<Number> &) const
368 : {
369 0 : libmesh_not_implemented();
370 : }
371 :
372 1050987 : void StaticCondensation::set_local_vectors(const NumericVector<Number> & global_vector,
373 : const std::vector<dof_id_type> & elem_dof_indices,
374 : std::vector<Number> & elem_dof_values_vec,
375 : EigenVector & elem_dof_values)
376 : {
377 1050987 : global_vector.get(elem_dof_indices, elem_dof_values_vec);
378 177510 : elem_dof_values.resize(elem_dof_indices.size());
379 7419614 : for (const auto i : index_range(elem_dof_indices))
380 7454643 : elem_dof_values(i) = elem_dof_values_vec[i];
381 1050987 : }
382 :
383 8330 : void StaticCondensation::forward_elimination(const NumericVector<Number> & full_rhs)
384 : {
385 476 : std::vector<dof_id_type> elem_condensed_dofs;
386 476 : std::vector<Number> elem_condensed_rhs_vec;
387 476 : EigenVector elem_condensed_rhs, elem_uncondensed_rhs;
388 :
389 8330 : full_rhs.create_subvector(
390 8330 : *_reduced_rhs, _reduced_dof_map._local_uncondensed_dofs, /*all_global_entries=*/false);
391 :
392 476 : std::vector<dof_id_type> reduced_space_indices;
393 949069 : for (auto elem : _mesh.active_local_element_ptr_range())
394 : {
395 350329 : auto & matrix_data = libmesh_map_find(_elem_to_matrix_data, elem->id());
396 29585 : reduced_space_indices.clear();
397 350329 : const auto & dof_data = libmesh_map_find(_reduced_dof_map._elem_to_dof_data, elem->id());
398 706994 : for (const auto & var_reduced_space_indices : dof_data.reduced_space_indices)
399 326504 : reduced_space_indices.insert(reduced_space_indices.end(),
400 : var_reduced_space_indices.begin(),
401 120644 : var_reduced_space_indices.end());
402 350329 : elem_condensed_dofs.resize(dof_data.condensed_global_to_local_map.size());
403 1339041 : for (const auto & [global_dof, local_dof] : dof_data.condensed_global_to_local_map)
404 : {
405 85198 : libmesh_assert(local_dof < elem_condensed_dofs.size());
406 1073910 : elem_condensed_dofs[local_dof] = global_dof;
407 : }
408 :
409 350329 : set_local_vectors(full_rhs, elem_condensed_dofs, elem_condensed_rhs_vec, elem_condensed_rhs);
410 350329 : elem_uncondensed_rhs = -matrix_data.Auc * matrix_data.AccFactor.solve(elem_condensed_rhs);
411 :
412 29585 : libmesh_assert(cast_int<std::size_t>(elem_uncondensed_rhs.size()) ==
413 : reduced_space_indices.size());
414 379914 : _reduced_rhs->add_vector(elem_uncondensed_rhs.data(), reduced_space_indices);
415 7854 : }
416 8330 : _reduced_rhs->close();
417 8330 : }
418 :
419 8330 : void StaticCondensation::backwards_substitution(const NumericVector<Number> & full_rhs,
420 : NumericVector<Number> & full_sol)
421 : {
422 476 : std::vector<dof_id_type> elem_condensed_dofs, elem_uncondensed_dofs;
423 476 : std::vector<Number> elem_condensed_rhs_vec, elem_uncondensed_sol_vec;
424 476 : EigenVector elem_condensed_rhs, elem_uncondensed_sol, elem_condensed_sol;
425 :
426 949069 : for (auto elem : _mesh.active_local_element_ptr_range())
427 : {
428 350329 : auto & matrix_data = libmesh_map_find(_elem_to_matrix_data, elem->id());
429 350329 : const auto & dof_data = libmesh_map_find(_reduced_dof_map._elem_to_dof_data, elem->id());
430 350329 : elem_condensed_dofs.resize(dof_data.condensed_global_to_local_map.size());
431 350329 : elem_uncondensed_dofs.resize(dof_data.uncondensed_global_to_local_map.size());
432 1339041 : for (const auto & [global_dof, local_dof] : dof_data.condensed_global_to_local_map)
433 : {
434 85198 : libmesh_assert(local_dof < elem_condensed_dofs.size());
435 1073910 : elem_condensed_dofs[local_dof] = global_dof;
436 : }
437 4741532 : for (const auto & [global_dof, local_dof] : dof_data.uncondensed_global_to_local_map)
438 : {
439 372612 : libmesh_assert(local_dof < elem_uncondensed_dofs.size());
440 4763815 : elem_uncondensed_dofs[local_dof] = global_dof;
441 : }
442 :
443 350329 : set_local_vectors(full_rhs, elem_condensed_dofs, elem_condensed_rhs_vec, elem_condensed_rhs);
444 350329 : set_local_vectors(*_ghosted_full_sol,
445 : elem_uncondensed_dofs,
446 : elem_uncondensed_sol_vec,
447 : elem_uncondensed_sol);
448 :
449 : elem_condensed_sol =
450 379914 : matrix_data.AccFactor.solve(elem_condensed_rhs - matrix_data.Acu * elem_uncondensed_sol);
451 379914 : full_sol.insert(elem_condensed_sol.data(), elem_condensed_dofs);
452 7854 : }
453 :
454 8330 : full_sol.close();
455 8330 : }
456 :
457 8330 : void StaticCondensation::apply(const NumericVector<Number> & full_rhs,
458 : NumericVector<Number> & full_parallel_sol)
459 : {
460 8330 : forward_elimination(full_rhs);
461 : // Apparently PETSc will send us the yvec without zeroing it ahead of time. This can be a poor
462 : // initial guess for the Krylov solve as well as lead to bewildered users who expect their initial
463 : // residual norm to equal the norm of the RHS
464 8330 : full_parallel_sol.zero();
465 16422 : _reduced_sol = full_parallel_sol.get_subvector(_reduced_dof_map._local_uncondensed_dofs);
466 8568 : _reduced_solver->solve(*_reduced_sys_mat, *_reduced_sol, *_reduced_rhs, 1e-5, 300);
467 : // Must restore to the full solution because during backwards substitution we will need to be able
468 : // to read ghosted dofs and we don't support ghosting of subvectors
469 8568 : full_parallel_sol.restore_subvector(std::move(_reduced_sol),
470 8330 : _reduced_dof_map._local_uncondensed_dofs);
471 8330 : *_ghosted_full_sol = full_parallel_sol;
472 8330 : backwards_substitution(full_rhs, full_parallel_sol);
473 8330 : }
474 :
475 0 : SolverPackage StaticCondensation::solver_package() { return libMesh::default_solver_package(); }
476 :
477 280 : void StaticCondensation::dont_condense_vars(const std::unordered_set<unsigned int> & vars)
478 : {
479 280 : _reduced_dof_map.dont_condense_vars(vars);
480 280 : }
481 :
482 : }
483 : #else
484 :
485 : #include "libmesh/dof_map.h"
486 :
487 : namespace libMesh
488 : {
489 0 : StaticCondensation::StaticCondensation(const MeshBase &,
490 : const System &,
491 : const DofMap & full_dof_map,
492 0 : StaticCondensationDofMap &)
493 0 : : SparseMatrix<Number>(full_dof_map.comm())
494 : {
495 0 : libmesh_error_msg(
496 : "Static condensation requires configuring libMesh with PETSc and Eigen support");
497 : }
498 : }
499 :
500 : #endif // LIBMESH_HAVE_EIGEN && LIBMESH_HAVE_PETSC
|
