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 60 : StaticCondensation::StaticCondensation(const MeshBase & mesh,
40 : System & system,
41 : const DofMap & full_dof_map,
42 60 : StaticCondensationDofMap & reduced_dof_map)
43 : : PetscMatrixShellMatrix<Number>(full_dof_map.comm()),
44 20 : _mesh(mesh),
45 20 : _system(system),
46 20 : _reduced_dof_map(reduced_dof_map),
47 20 : _current_elem_id(DofObject::invalid_id),
48 20 : _sc_is_initialized(false),
49 20 : _have_cached_values(false),
50 20 : _parallel_type(INVALID_PARALLELIZATION),
51 60 : _uncondensed_dofs_only(false)
52 : {
53 40 : _size_one_mat.resize(1, 1);
54 80 : _scp = std::make_unique<StaticCondensationPreconditioner>(*this);
55 60 : }
56 :
57 40 : StaticCondensation::~StaticCondensation() = default;
58 :
59 0 : SparseMatrix<Number> & StaticCondensation::operator=(const SparseMatrix<Number> &)
60 : {
61 0 : libmesh_not_implemented();
62 : }
63 :
64 0 : std::unique_ptr<SparseMatrix<Number>> StaticCondensation::zero_clone() const
65 : {
66 0 : libmesh_not_implemented();
67 : }
68 :
69 0 : std::unique_ptr<SparseMatrix<Number>> StaticCondensation::clone() const
70 : {
71 0 : libmesh_not_implemented();
72 : }
73 :
74 360 : void StaticCondensation::clear() noexcept
75 : {
76 360 : PetscMatrixShellMatrix<Number>::clear();
77 :
78 120 : _elem_to_matrix_data.clear();
79 120 : _reduced_sys_mat.reset();
80 120 : _reduced_sol.reset();
81 120 : _reduced_rhs.reset();
82 120 : _reduced_solver.reset();
83 360 : _current_elem_id = DofObject::invalid_id;
84 360 : _have_cached_values = false;
85 360 : _sc_is_initialized = false;
86 360 : }
87 :
88 0 : void StaticCondensation::init(const numeric_index_type m,
89 : const numeric_index_type n,
90 : const numeric_index_type m_l,
91 : const numeric_index_type n_l,
92 : const numeric_index_type nnz,
93 : const numeric_index_type noz,
94 : const numeric_index_type blocksize)
95 : {
96 0 : if (!this->initialized())
97 : {
98 0 : PetscMatrixShellMatrix<Number>::init(m, n, m_l, n_l, nnz, noz, blocksize);
99 0 : _parallel_type = ((m == m_l) && (this->n_processors() > 1)) ? SERIAL : PARALLEL;
100 0 : this->init();
101 : }
102 0 : }
103 :
104 360 : void StaticCondensation::init(const ParallelType type)
105 : {
106 360 : if (!this->initialized())
107 : {
108 : // The operator acts on the parent space of the static condensation dof map
109 360 : const auto n = _reduced_dof_map.parent_n_dofs();
110 360 : const auto n_local = _reduced_dof_map.parent_n_local_dofs();
111 120 : libmesh_assert(n == n_local || type != SERIAL);
112 360 : PetscMatrixShellMatrix<Number>::init(n, n, n_local, n_local);
113 360 : _parallel_type = type;
114 360 : this->init();
115 : }
116 360 : }
117 :
118 3182 : bool StaticCondensation::initialized() const
119 : {
120 3182 : return PetscMatrixShellMatrix<Number>::initialized() && _sc_is_initialized;
121 : }
122 :
123 1368 : void StaticCondensation::init()
124 : {
125 1368 : if (this->initialized())
126 336 : return;
127 :
128 120 : libmesh_assert(_reduced_dof_map.initialized());
129 :
130 : // This API is public, so it can be called without going through the other init overloads which
131 : // would give us an indication of what kind of parallel type the user wants. If we've gotten no
132 : // indication so far, we will default to parallel
133 360 : if (_parallel_type == INVALID_PARALLELIZATION)
134 0 : _parallel_type = PARALLEL;
135 120 : libmesh_assert(_parallel_type != SERIAL);
136 :
137 47583 : for (const auto & [elem_id, dof_data] : _reduced_dof_map._elem_to_dof_data)
138 : {
139 15741 : auto & matrix_data = _elem_to_matrix_data[elem_id];
140 :
141 15741 : const auto condensed_dof_size = dof_data.condensed_global_to_local_map.size();
142 15741 : const auto uncondensed_dof_size = dof_data.uncondensed_global_to_local_map.size();
143 :
144 47223 : matrix_data.Acc.setZero(condensed_dof_size, condensed_dof_size);
145 47223 : matrix_data.Acu.setZero(condensed_dof_size, uncondensed_dof_size);
146 47223 : matrix_data.Auc.setZero(uncondensed_dof_size, condensed_dof_size);
147 47223 : matrix_data.Auu.setZero(uncondensed_dof_size, uncondensed_dof_size);
148 : }
149 :
150 : //
151 : // Build the reduced system data
152 : //
153 360 : const auto n = _reduced_dof_map.n_dofs();
154 : const auto n_local =
155 360 : (_parallel_type == SERIAL) ? _reduced_dof_map.n_dofs() : _reduced_dof_map.n_local_dofs();
156 600 : _reduced_solver = LinearSolver<Number>::build(this->comm());
157 480 : _reduced_solver->init((_system.name() + "_condensed_").c_str());
158 600 : _reduced_rhs = NumericVector<Number>::build(this->comm());
159 : // Init the RHS vector so we can conveniently get processor row offsets
160 360 : _reduced_rhs->init(n, n_local);
161 :
162 : // Initialize the reduced system matrix
163 360 : _sp = &_reduced_dof_map.reduced_sparsity();
164 600 : _reduced_sys_mat = SparseMatrix<Number>::build(this->comm());
165 360 : if (auto * const petsc_mat = dynamic_cast<PetscMatrix<Number> *>(_reduced_sys_mat.get()))
166 : {
167 : // Optimization for PETSc. This is critical for problems in which there are SCALAR dofs that
168 : // introduce dense rows to avoid allocating a dense matrix
169 360 : petsc_mat->init(
170 360 : n, n, n_local, n_local, _reduced_dof_map.reduced_n_nz(), _reduced_dof_map.reduced_n_oz());
171 : }
172 : else
173 : {
174 0 : const auto & nnz = _sp->get_n_nz();
175 0 : const auto & noz = _sp->get_n_oz();
176 0 : const auto nz = nnz.empty() ? dof_id_type(0) : *std::max_element(nnz.begin(), nnz.end());
177 0 : const auto oz = noz.empty() ? dof_id_type(0) : *std::max_element(noz.begin(), noz.end());
178 0 : _reduced_sys_mat->init(n, n, n_local, n_local, nz, oz);
179 : }
180 :
181 : // Build the ghosted parent-space solution vector. Note that this is, in general, *not equal* to
182 : // the system solution, e.g. this may correspond to the solution for the Newton *update*
183 600 : _ghosted_parent_sol = NumericVector<Number>::build(this->comm());
184 480 : _ghosted_parent_sol->init(_reduced_dof_map.parent_n_dofs(),
185 360 : _reduced_dof_map.parent_n_local_dofs(),
186 360 : _reduced_dof_map.parent_send_list(),
187 : /*fast=*/false,
188 240 : GHOSTED);
189 :
190 360 : _sc_is_initialized = true;
191 : }
192 :
193 1056 : void StaticCondensation::close()
194 : {
195 1056 : _communicator.max(_have_cached_values);
196 1056 : if (!_have_cached_values)
197 : {
198 294 : bool closed = _reduced_sys_mat->closed();
199 : // closed is not collective
200 294 : _communicator.min(closed);
201 294 : if (!closed)
202 0 : _reduced_sys_mat->close();
203 98 : return;
204 : }
205 :
206 1270 : DenseMatrix<Number> shim;
207 254 : std::vector<dof_id_type> reduced_space_indices;
208 89847 : for (auto & [elem_id, matrix_data] : _elem_to_matrix_data)
209 : {
210 29695 : reduced_space_indices.clear();
211 :
212 : // The result matrix is either a Schur complement or it's simply the result of summing element
213 : // matrices of the uncondensed degrees of freedom
214 59390 : EigenMatrix result = matrix_data.Auu;
215 89085 : if (!_uncondensed_dofs_only)
216 : {
217 89085 : matrix_data.AccFactor = matrix_data.Acc.partialPivLu();
218 89085 : result -= matrix_data.Auc * matrix_data.AccFactor.solve(matrix_data.Acu);
219 : }
220 118780 : shim.resize(result.rows(), result.cols());
221 1210461 : for (const auto i : make_range(result.rows()))
222 15839424 : for (const auto j : make_range(result.cols()))
223 19624064 : shim(i, j) = result(i, j);
224 179898 : for (const auto & var_reduced_space_indices : _reduced_dof_map.elem_dof_indices(elem_id))
225 60542 : reduced_space_indices.insert(reduced_space_indices.end(),
226 : var_reduced_space_indices.begin(),
227 121084 : var_reduced_space_indices.end());
228 89085 : _reduced_sys_mat->add_matrix(shim, reduced_space_indices);
229 : }
230 :
231 762 : _reduced_sys_mat->close();
232 :
233 762 : _have_cached_values = false;
234 254 : }
235 :
236 120 : bool StaticCondensation::closed() const
237 : {
238 120 : return _reduced_sys_mat->closed() && !_have_cached_values;
239 : }
240 :
241 1122 : void StaticCondensation::zero()
242 : {
243 1122 : _reduced_sys_mat->zero();
244 137430 : for (auto & [elem_id, matrix_data] : _elem_to_matrix_data)
245 : {
246 45436 : libmesh_ignore(elem_id);
247 90872 : matrix_data.Acc.setZero();
248 90872 : matrix_data.Acu.setZero();
249 90872 : matrix_data.Auc.setZero();
250 90872 : matrix_data.Auu.setZero();
251 : }
252 1122 : }
253 :
254 360 : void StaticCondensation::setup() { libmesh_assert(this->closed()); }
255 :
256 0 : numeric_index_type StaticCondensation::m() const { return _reduced_dof_map.parent_n_dofs(); }
257 :
258 0 : numeric_index_type StaticCondensation::row_start() const
259 : {
260 0 : return _reduced_dof_map.parent_first_dof();
261 : }
262 :
263 0 : numeric_index_type StaticCondensation::row_stop() const
264 : {
265 0 : return _reduced_dof_map.parent_end_dof();
266 : }
267 :
268 0 : void StaticCondensation::set(const numeric_index_type full_i,
269 : const numeric_index_type full_j,
270 : const Number val)
271 : {
272 0 : const auto reduced_i = _reduced_dof_map.find_reduced_dof(full_i);
273 0 : const auto reduced_j = _reduced_dof_map.find_reduced_dof(full_j);
274 0 : libmesh_assert_msg((reduced_i != DofObject::invalid_id && reduced_j != DofObject::invalid_id) ||
275 : _reduced_dof_map.parent(),
276 : "Only uncondensed entries may be set in a statically condensed matrix");
277 0 : if (reduced_i != DofObject::invalid_id && reduced_j != DofObject::invalid_id)
278 0 : _reduced_sys_mat->set(reduced_i, reduced_j, val);
279 0 : }
280 :
281 89085 : void StaticCondensation::set_current_elem(const Elem & elem)
282 : {
283 29695 : libmesh_assert(!Threads::in_threads || libMesh::n_threads() == 1);
284 89085 : _current_elem_id = elem.id();
285 89085 : }
286 :
287 0 : void StaticCondensation::add(const numeric_index_type i,
288 : const numeric_index_type j,
289 : const Number value)
290 : {
291 0 : _size_one_mat(0, 0) = value;
292 0 : this->add_matrix(_size_one_mat, {i}, {j});
293 0 : }
294 :
295 113277 : void StaticCondensation::add_matrix(const DenseMatrix<Number> & dm,
296 : const std::vector<numeric_index_type> & rows,
297 : const std::vector<numeric_index_type> & cols)
298 : {
299 113277 : if (rows.empty() || cols.empty())
300 0 : return;
301 :
302 37759 : libmesh_assert(_current_elem_id != DofObject::invalid_id);
303 113277 : auto & matrix_data = libmesh_map_find(_elem_to_matrix_data, _current_elem_id);
304 113277 : const auto & dof_data = libmesh_map_find(_reduced_dof_map._elem_to_dof_data, _current_elem_id);
305 : EigenMatrix * mat;
306 :
307 : // Entries in rows or columns the parent space eliminates are dropped; the element matrix already
308 : // carries the parent reduction, e.g. C^T K C for eliminated constraints
309 113277 : const auto * const parent = _reduced_dof_map.parent();
310 30046360 : auto in_parent = [parent](const auto global_index) {
311 37446992 : return !parent || parent->find_reduced_dof(global_index) != DofObject::invalid_id;
312 75518 : };
313 :
314 89732640 : auto info_from_index = [&dof_data](const auto global_index) {
315 14955440 : auto index_it = dof_data.condensed_global_to_local_map.find(global_index);
316 14955440 : const bool index_is_condensed = index_it != dof_data.condensed_global_to_local_map.end();
317 44866320 : if (!index_is_condensed)
318 : {
319 12045728 : index_it = dof_data.uncondensed_global_to_local_map.find(global_index);
320 36137184 : if (index_it == dof_data.uncondensed_global_to_local_map.end())
321 0 : libmesh_error_msg("Failed to find the global index "
322 : << global_index
323 : << " in our current element's degree of freedom information. One way "
324 : "this can happen is when using a discontinuous Galerkin method, "
325 : "adding element matrices to both + and - sides of a face");
326 : }
327 : else
328 : // We found the dof in the condensed container. Let's assert that it's not also in the
329 : // uncondensed container
330 2909712 : libmesh_assert(dof_data.uncondensed_global_to_local_map.find(global_index) ==
331 : dof_data.uncondensed_global_to_local_map.end());
332 :
333 74777200 : return std::make_pair(index_is_condensed, index_it->second);
334 75518 : };
335 :
336 1579089 : for (const auto i : make_range(dm.m()))
337 23944500 : for (const auto j : make_range(dm.n()))
338 : {
339 22478688 : const auto global_i = rows[i];
340 29971584 : const auto global_j = cols[j];
341 14985792 : if (!in_parent(global_i) || !in_parent(global_j))
342 45528 : continue;
343 22433160 : const auto [i_is_condensed, local_i] = info_from_index(global_i);
344 22433160 : const auto [j_is_condensed, local_j] = info_from_index(global_j);
345 22433160 : if (i_is_condensed)
346 : {
347 4378392 : if (j_is_condensed)
348 1002312 : mat = &matrix_data.Acc;
349 : else
350 3376080 : mat = &matrix_data.Acu;
351 : }
352 : else
353 : {
354 18054768 : if (j_is_condensed)
355 3348432 : mat = &matrix_data.Auc;
356 : else
357 14706336 : mat = &matrix_data.Auu;
358 : }
359 29910880 : (*mat)(local_i, local_j) += dm(i, j);
360 : }
361 :
362 113277 : _have_cached_values = true;
363 : }
364 :
365 88317 : void StaticCondensation::add_matrix(const DenseMatrix<Number> & dm,
366 : const std::vector<numeric_index_type> & dof_indices)
367 : {
368 88317 : this->add_matrix(dm, dof_indices, dof_indices);
369 88317 : }
370 :
371 0 : void StaticCondensation::add(const Number, const SparseMatrix<Number> &)
372 : {
373 0 : libmesh_not_implemented();
374 : }
375 :
376 0 : Number StaticCondensation::operator()(const numeric_index_type, const numeric_index_type) const
377 : {
378 0 : libmesh_not_implemented();
379 : }
380 :
381 0 : Real StaticCondensation::l1_norm() const { libmesh_not_implemented(); }
382 :
383 0 : Real StaticCondensation::linfty_norm() const { libmesh_not_implemented(); }
384 :
385 0 : void StaticCondensation::print_personal(std::ostream &) const { libmesh_not_implemented(); }
386 :
387 0 : void StaticCondensation::get_diagonal(NumericVector<Number> &) const { libmesh_not_implemented(); }
388 :
389 0 : void StaticCondensation::get_transpose(SparseMatrix<Number> &) const { libmesh_not_implemented(); }
390 :
391 0 : void StaticCondensation::get_row(numeric_index_type,
392 : std::vector<numeric_index_type> &,
393 : std::vector<Number> &) const
394 : {
395 0 : libmesh_not_implemented();
396 : }
397 :
398 267084 : void StaticCondensation::set_local_vectors(const NumericVector<Number> & global_vector,
399 : const std::vector<dof_id_type> & elem_dof_indices,
400 : std::vector<Number> & elem_dof_values_vec,
401 : EigenVector & elem_dof_values)
402 : {
403 267084 : global_vector.get(elem_dof_indices, elem_dof_values_vec);
404 178056 : elem_dof_values.resize(elem_dof_indices.size());
405 1901952 : for (const auto i : index_range(elem_dof_indices))
406 2724780 : elem_dof_values(i) = elem_dof_values_vec[i];
407 267084 : }
408 :
409 756 : void StaticCondensation::forward_elimination(const NumericVector<Number> & parent_rhs)
410 : {
411 504 : std::vector<Number> elem_condensed_rhs_vec;
412 504 : EigenVector elem_condensed_rhs, elem_uncondensed_rhs;
413 :
414 756 : parent_rhs.create_subvector(
415 756 : *_reduced_rhs, _reduced_dof_map.local_retained_parent_dofs(), /*all_global_entries=*/false);
416 :
417 504 : std::vector<dof_id_type> reduced_space_indices;
418 149640 : for (auto elem : _mesh.active_local_element_ptr_range())
419 : {
420 89028 : auto & matrix_data = libmesh_map_find(_elem_to_matrix_data, elem->id());
421 29676 : reduced_space_indices.clear();
422 89028 : const auto & dof_data = libmesh_map_find(_reduced_dof_map._elem_to_dof_data, elem->id());
423 179784 : for (const auto & var_reduced_space_indices : _reduced_dof_map.elem_dof_indices(elem->id()))
424 60504 : reduced_space_indices.insert(reduced_space_indices.end(),
425 : var_reduced_space_indices.begin(),
426 121008 : var_reduced_space_indices.end());
427 :
428 89028 : set_local_vectors(
429 89028 : parent_rhs, dof_data.condensed_parent_dofs, elem_condensed_rhs_vec, elem_condensed_rhs);
430 89028 : elem_uncondensed_rhs = -matrix_data.Auc * matrix_data.AccFactor.solve(elem_condensed_rhs);
431 :
432 29676 : libmesh_assert(cast_int<std::size_t>(elem_uncondensed_rhs.size()) ==
433 : reduced_space_indices.size());
434 118704 : _reduced_rhs->add_vector(elem_uncondensed_rhs.data(), reduced_space_indices);
435 252 : }
436 756 : _reduced_rhs->close();
437 756 : }
438 :
439 756 : void StaticCondensation::backwards_substitution(const NumericVector<Number> & parent_rhs,
440 : NumericVector<Number> & parent_sol)
441 : {
442 504 : std::vector<Number> elem_condensed_rhs_vec, elem_uncondensed_sol_vec;
443 504 : EigenVector elem_condensed_rhs, elem_uncondensed_sol, elem_condensed_sol;
444 :
445 149640 : for (auto elem : _mesh.active_local_element_ptr_range())
446 : {
447 89028 : auto & matrix_data = libmesh_map_find(_elem_to_matrix_data, elem->id());
448 89028 : const auto & dof_data = libmesh_map_find(_reduced_dof_map._elem_to_dof_data, elem->id());
449 :
450 89028 : set_local_vectors(
451 89028 : parent_rhs, dof_data.condensed_parent_dofs, elem_condensed_rhs_vec, elem_condensed_rhs);
452 89028 : set_local_vectors(*_ghosted_parent_sol,
453 89028 : dof_data.uncondensed_parent_dofs,
454 : elem_uncondensed_sol_vec,
455 : elem_uncondensed_sol);
456 :
457 : elem_condensed_sol =
458 118704 : matrix_data.AccFactor.solve(elem_condensed_rhs - matrix_data.Acu * elem_uncondensed_sol);
459 118704 : parent_sol.insert(elem_condensed_sol.data(), dof_data.condensed_parent_dofs);
460 252 : }
461 :
462 756 : parent_sol.close();
463 756 : }
464 :
465 756 : void StaticCondensation::apply(const NumericVector<Number> & parent_rhs,
466 : NumericVector<Number> & parent_sol)
467 : {
468 756 : forward_elimination(parent_rhs);
469 : // Apparently PETSc will send us the yvec without zeroing it ahead of time. This can be a poor
470 : // initial guess for the Krylov solve as well as lead to bewildered users who expect their initial
471 : // residual norm to equal the norm of the RHS
472 756 : parent_sol.zero();
473 1512 : _reduced_sol = parent_sol.get_subvector(_reduced_dof_map.local_retained_parent_dofs());
474 1008 : _reduced_solver->solve(*_reduced_sys_mat, *_reduced_sol, *_reduced_rhs, 1e-5, 300);
475 : // Must restore to the parent-space solution because during backwards substitution we will need
476 : // to be able to read ghosted dofs and we don't support ghosting of subvectors
477 1008 : parent_sol.restore_subvector(std::move(_reduced_sol),
478 756 : _reduced_dof_map.local_retained_parent_dofs());
479 756 : *_ghosted_parent_sol = parent_sol;
480 756 : backwards_substitution(parent_rhs, parent_sol);
481 756 : }
482 :
483 0 : SolverPackage StaticCondensation::solver_package() { return libMesh::default_solver_package(); }
484 :
485 24 : void StaticCondensation::dont_condense_vars(const std::unordered_set<unsigned int> & vars)
486 : {
487 24 : _reduced_dof_map.dont_condense_vars(vars);
488 24 : }
489 :
490 : }
491 : #else
492 :
493 : #include "libmesh/dof_map.h"
494 :
495 : namespace libMesh
496 : {
497 0 : StaticCondensation::StaticCondensation(const MeshBase &,
498 : const System &,
499 : const DofMap & full_dof_map,
500 0 : StaticCondensationDofMap &)
501 0 : : SparseMatrix<Number>(full_dof_map.comm())
502 : {
503 0 : libmesh_error_msg(
504 : "Static condensation requires configuring libMesh with PETSc and Eigen support");
505 : }
506 : }
507 :
508 : #endif // LIBMESH_HAVE_EIGEN && LIBMESH_HAVE_PETSC
|
