Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 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 : #ifndef LIBMESH_STATIC_CONDENSATION_H
19 : #define LIBMESH_STATIC_CONDENSATION_H
20 :
21 : #include "libmesh/libmesh_config.h"
22 :
23 : // subvectors currently only work with petsc and we rely on Eigen for local LU factorizations
24 : #if defined(LIBMESH_HAVE_EIGEN) && defined(LIBMESH_HAVE_PETSC)
25 :
26 : #include "libmesh/petsc_matrix_shell_matrix.h"
27 : #include "libmesh/id_types.h"
28 : #include "libmesh/libmesh_common.h"
29 : #include "libmesh/dense_matrix.h"
30 : #include "libmesh/variable.h"
31 : #include "libmesh/sparsity_pattern.h"
32 :
33 : #include <unordered_map>
34 : #include <memory>
35 : #include <vector>
36 :
37 : // Warnings about stack protection
38 : #include "libmesh/ignore_warnings.h"
39 : #include <Eigen/Dense>
40 : #include "libmesh/restore_warnings.h"
41 :
42 : namespace libMesh
43 : {
44 : class MeshBase;
45 : class System;
46 : class DofMap;
47 : class Elem;
48 : template <typename>
49 : class LinearSolver;
50 : template <typename>
51 : class NumericVector;
52 : template <typename>
53 : class Preconditioner;
54 : class StaticCondensationPreconditioner;
55 : class StaticCondensationDofMap;
56 :
57 : typedef Eigen::Matrix<Number, Eigen::Dynamic, Eigen::Dynamic> EigenMatrix;
58 : typedef Eigen::Matrix<Number, Eigen::Dynamic, 1> EigenVector;
59 :
60 50 : class StaticCondensation : public PetscMatrixShellMatrix<Number>
61 : {
62 : public:
63 : StaticCondensation(const MeshBase & mesh,
64 : System & system,
65 : const DofMap & full_dof_map,
66 : StaticCondensationDofMap & reduced_dof_map);
67 : virtual ~StaticCondensation();
68 :
69 : //
70 : // SparseMatrix overrides
71 : //
72 :
73 : virtual SparseMatrix<Number> & operator=(const SparseMatrix<Number> &) override;
74 :
75 : virtual SolverPackage solver_package() override;
76 :
77 : virtual void init(const numeric_index_type m,
78 : const numeric_index_type n,
79 : const numeric_index_type m_l,
80 : const numeric_index_type n_l,
81 : const numeric_index_type nnz = 30,
82 : const numeric_index_type noz = 10,
83 : const numeric_index_type blocksize = 1) override;
84 :
85 : virtual void init(const ParallelType type) override;
86 :
87 : virtual bool initialized() const override;
88 :
89 : virtual void clear() noexcept override;
90 :
91 : virtual void zero() override;
92 :
93 : virtual std::unique_ptr<SparseMatrix<Number>> zero_clone() const override;
94 :
95 : virtual std::unique_ptr<SparseMatrix<Number>> clone() const override;
96 :
97 : virtual void close() override;
98 :
99 : virtual numeric_index_type m() const override;
100 :
101 0 : virtual numeric_index_type n() const override { return this->m(); }
102 :
103 : virtual numeric_index_type row_start() const override;
104 :
105 : virtual numeric_index_type row_stop() const override;
106 :
107 0 : virtual numeric_index_type col_start() const override { return this->row_start(); }
108 :
109 0 : virtual numeric_index_type col_stop() const override { return this->row_stop(); }
110 :
111 : virtual void
112 : set(const numeric_index_type i, const numeric_index_type j, const Number value) override;
113 :
114 : virtual void
115 : add(const numeric_index_type i, const numeric_index_type j, const Number value) override;
116 :
117 : virtual void add_matrix(const DenseMatrix<Number> & dm,
118 : const std::vector<numeric_index_type> & rows,
119 : const std::vector<numeric_index_type> & cols) override;
120 :
121 : virtual void add_matrix(const DenseMatrix<Number> & dm,
122 : const std::vector<numeric_index_type> & dof_indices) override;
123 :
124 : virtual void add(const Number a, const SparseMatrix<Number> & X) override;
125 :
126 : virtual Number operator()(const numeric_index_type i, const numeric_index_type j) const override;
127 :
128 : virtual Real l1_norm() const override;
129 :
130 : virtual Real linfty_norm() const override;
131 :
132 : virtual bool closed() const override;
133 :
134 : virtual void print_personal(std::ostream & os = libMesh::out) const override;
135 :
136 : virtual void get_diagonal(NumericVector<Number> & dest) const override;
137 :
138 : virtual void get_transpose(SparseMatrix<Number> & dest) const override;
139 :
140 : virtual void get_row(numeric_index_type i,
141 : std::vector<numeric_index_type> & indices,
142 : std::vector<Number> & values) const override;
143 :
144 : //
145 : // Unique methods
146 : //
147 :
148 : /**
149 : * Size the element matrices
150 : */
151 : void init();
152 :
153 : /**
154 : * A no-op to be consistent with shimming from the \p StaticCondenstionPreconditioner. "setup"
155 : * should take place at the end of matrix assembly, e.g. during a call to \p close()
156 : */
157 : void setup();
158 :
159 : /**
160 : * Perform our three stages, \p forward_elimination(), a \p solve() on the condensed system, and
161 : * finally a \p backwards_substitution()
162 : */
163 : void apply(const NumericVector<Number> & full_rhs, NumericVector<Number> & full_sol);
164 :
165 : const SparseMatrix<Number> & get_condensed_mat() const;
166 :
167 : /**
168 : * Set the current element. This enables fast lookups of local indices from global indices
169 : */
170 : void set_current_elem(const Elem & elem);
171 :
172 : /**
173 : * Get the preconditioning wrapper
174 : */
175 10 : StaticCondensationPreconditioner & get_preconditioner() { return *_scp; }
176 :
177 : /**
178 : * @returns The reduced system linear solver
179 : */
180 : LinearSolver<Number> & reduced_system_solver();
181 :
182 350 : virtual bool require_sparsity_pattern() const override { return false; }
183 :
184 : /**
185 : * Sets whether this matrix represents uncondensed dofs only. In that case when building the Schur
186 : * complement we won't attempt to invert zero element matrices corresponding to the condensed dofs
187 : */
188 : void uncondensed_dofs_only() { _uncondensed_dofs_only = true; }
189 :
190 : /**
191 : * Add \p vars to the list of variables not to condense. This can be useful when some variable's
192 : * equation is discretized with a DG method or if including the variable in the condensed block
193 : * diagonal would result in it being singular
194 : */
195 : void dont_condense_vars(const std::unordered_set<unsigned int> & vars);
196 :
197 : private:
198 : /**
199 : * Retrieves the degree of freedom values from \p global_vector corresponding to \p
200 : * elem_dof_indices, filling both \p elem_dof_values_vec and \p elem_dof_values
201 : */
202 : static void set_local_vectors(const NumericVector<Number> & global_vector,
203 : const std::vector<dof_id_type> & elem_dof_indices,
204 : std::vector<Number> & elem_dof_values_vec,
205 : EigenVector & elem_dof_values);
206 :
207 : /**
208 : * Takes an incoming "full" RHS from the solver, where full means the union of uncondensed and
209 : * condennsed dofs, and condenses it down into the condensed \p _reduced_rhs data member using
210 : * element Schur complements
211 : */
212 : void forward_elimination(const NumericVector<Number> & full_rhs);
213 :
214 : /**
215 : * After performing the solve with the \p _reduced_rhs and Schur complement matrix (\p
216 : * _reduced_sys_mat) to determine the \p _reduced_sol, we use the uncondensed \p full_rhs along
217 : * with the \p _reduced_sol to back substitute into the \p full_sol, which is the final output
218 : * data of the static condensation preconditioner application
219 : */
220 : void backwards_substitution(const NumericVector<Number> & full_rhs,
221 : NumericVector<Number> & full_sol);
222 :
223 : /**
224 : * Data stored on a per-element basis used to compute element Schur complements and their
225 : * applications to vectors
226 : */
227 : struct MatrixData
228 : {
229 : /// condensed-condensed matrix entries
230 : EigenMatrix Acc;
231 : /// condensed-uncondensed matrix entries
232 : EigenMatrix Acu;
233 : /// uncondensed-condensed matrix entries
234 : EigenMatrix Auc;
235 : /// uncondensed-uncondensed matrix entries
236 : EigenMatrix Auu;
237 :
238 : // Acc LU decompositions
239 : typename std::remove_const<decltype(Acc.partialPivLu())>::type AccFactor;
240 : };
241 :
242 : /// A map from element ID to Schur complement data
243 : std::unordered_map<dof_id_type, MatrixData> _elem_to_matrix_data;
244 :
245 : const MeshBase & _mesh;
246 : System & _system;
247 : const DofMap & _full_dof_map;
248 : StaticCondensationDofMap & _reduced_dof_map;
249 :
250 : /// global sparse matrix for the uncondensed degrees of freedom
251 : std::unique_ptr<SparseMatrix<Number>> _reduced_sys_mat;
252 : /// solution for the uncondensed degrees of freedom
253 : std::unique_ptr<NumericVector<Number>> _reduced_sol;
254 : /// RHS corresponding to the uncondensed degrees of freedom. This is constructed by applying the
255 : /// element Schur complements to an incoming RHS which contains both condensed and uncondensed
256 : /// degrees of freedom
257 : std::unique_ptr<NumericVector<Number>> _reduced_rhs;
258 : /// The solver for the uncondensed degrees of freedom
259 : std::unique_ptr<LinearSolver<Number>> _reduced_solver;
260 : /// This is a ghosted representation of the full (uncondensed + condensed) solution. Note that
261 : // this is, in general, *not equal* to the system solution, e.g. this may correspond to the
262 : // solution for the Newton *update*
263 : std::unique_ptr<NumericVector<Number>> _ghosted_full_sol;
264 :
265 : /// The current element ID. This is one half of a key, along with the global index, that maps to a
266 : /// local index
267 : dof_id_type _current_elem_id;
268 :
269 : /// Helper data member for adding individual matrix elements
270 : DenseMatrix<Number> _size_one_mat;
271 :
272 : /// Preconditioner object which will call back to us for the preconditioning action
273 : std::unique_ptr<StaticCondensationPreconditioner> _scp;
274 :
275 : /// Whether our object has been initialized
276 : bool _sc_is_initialized;
277 :
278 : /// Whether we have cached values via add_XXX()
279 : bool _have_cached_values;
280 :
281 : /// The parallel type to use for the reduced matrix
282 : ParallelType _parallel_type;
283 :
284 : /// whether this matrix represents uncondensed dofs only. In that case when building the Schur
285 : /// complement we won't attempt to invert zero element matrices corresponding to the condensed
286 : /// dofs
287 : bool _uncondensed_dofs_only;
288 : };
289 :
290 : inline const SparseMatrix<Number> & StaticCondensation::get_condensed_mat() const
291 : {
292 : libmesh_assert(_reduced_sys_mat);
293 : return *_reduced_sys_mat;
294 : }
295 :
296 : inline LinearSolver<Number> & StaticCondensation::reduced_system_solver()
297 : {
298 : libmesh_assert_msg(_reduced_solver, "Reduced system solver not built yet");
299 : return *_reduced_solver;
300 : }
301 :
302 : }
303 :
304 : #else
305 :
306 : #include "libmesh/sparse_matrix.h"
307 : #include <unordered_set>
308 :
309 : namespace libMesh
310 : {
311 : class MeshBase;
312 : class System;
313 : class DofMap;
314 : class StaticCondensationPreconditioner;
315 : class StaticCondensationDofMap;
316 :
317 : class StaticCondensation : public SparseMatrix<Number>
318 : {
319 : public:
320 : StaticCondensation(const MeshBase &,
321 : const System &,
322 : const DofMap & full_dof_map,
323 : StaticCondensationDofMap & reduced_dof_map);
324 :
325 : const std::unordered_set<unsigned int> & uncondensed_vars() const { libmesh_not_implemented(); }
326 0 : StaticCondensationPreconditioner & get_preconditioner() { libmesh_not_implemented(); }
327 0 : virtual SparseMatrix<Number> & operator=(const SparseMatrix<Number> &) override
328 : {
329 0 : libmesh_not_implemented();
330 : }
331 0 : virtual SolverPackage solver_package() override { libmesh_not_implemented(); }
332 0 : virtual void init(const numeric_index_type,
333 : const numeric_index_type,
334 : const numeric_index_type,
335 : const numeric_index_type,
336 : const numeric_index_type = 30,
337 : const numeric_index_type = 10,
338 : const numeric_index_type = 1) override
339 : {
340 0 : libmesh_not_implemented();
341 : }
342 0 : virtual void init(ParallelType) override { libmesh_not_implemented(); }
343 0 : virtual void clear() override { libmesh_not_implemented(); }
344 0 : virtual void zero() override { libmesh_not_implemented(); }
345 0 : virtual std::unique_ptr<SparseMatrix<Number>> zero_clone() const override
346 : {
347 0 : libmesh_not_implemented();
348 : }
349 0 : virtual std::unique_ptr<SparseMatrix<Number>> clone() const override
350 : {
351 0 : libmesh_not_implemented();
352 : }
353 0 : virtual void close() override { libmesh_not_implemented(); }
354 0 : virtual numeric_index_type m() const override { libmesh_not_implemented(); }
355 0 : virtual numeric_index_type n() const override { libmesh_not_implemented(); }
356 0 : virtual numeric_index_type row_start() const override { libmesh_not_implemented(); }
357 0 : virtual numeric_index_type row_stop() const override { libmesh_not_implemented(); }
358 0 : virtual numeric_index_type col_start() const override { libmesh_not_implemented(); }
359 0 : virtual numeric_index_type col_stop() const override { libmesh_not_implemented(); }
360 0 : virtual void set(const numeric_index_type, const numeric_index_type, const Number) override
361 : {
362 0 : libmesh_not_implemented();
363 : }
364 0 : virtual void add(const numeric_index_type, const numeric_index_type, const Number) override
365 : {
366 0 : libmesh_not_implemented();
367 : }
368 0 : virtual void add_matrix(const DenseMatrix<Number> &,
369 : const std::vector<numeric_index_type> &,
370 : const std::vector<numeric_index_type> &) override
371 : {
372 0 : libmesh_not_implemented();
373 : }
374 0 : virtual void add_matrix(const DenseMatrix<Number> &,
375 : const std::vector<numeric_index_type> &) override
376 : {
377 0 : libmesh_not_implemented();
378 : }
379 0 : virtual void add(const Number, const SparseMatrix<Number> &) override
380 : {
381 0 : libmesh_not_implemented();
382 : }
383 0 : virtual Number operator()(const numeric_index_type, const numeric_index_type) const override
384 : {
385 0 : libmesh_not_implemented();
386 : }
387 0 : virtual Real l1_norm() const override { libmesh_not_implemented(); }
388 0 : virtual Real linfty_norm() const override { libmesh_not_implemented(); }
389 0 : virtual bool closed() const override { libmesh_not_implemented(); }
390 0 : virtual void print_personal(std::ostream & = libMesh::out) const override
391 : {
392 0 : libmesh_not_implemented();
393 : }
394 0 : virtual void get_diagonal(NumericVector<Number> &) const override { libmesh_not_implemented(); }
395 0 : virtual void get_transpose(SparseMatrix<Number> &) const override { libmesh_not_implemented(); }
396 0 : virtual void get_row(numeric_index_type,
397 : std::vector<numeric_index_type> &,
398 : std::vector<Number> &) const override
399 : {
400 0 : libmesh_not_implemented();
401 : }
402 0 : void init() { libmesh_not_implemented(); }
403 0 : void setup() { libmesh_not_implemented(); }
404 0 : void apply(const NumericVector<Number> &, NumericVector<Number> &) { libmesh_not_implemented(); }
405 : };
406 : }
407 :
408 : #endif // LIBMESH_HAVE_EIGEN && LIBMESH_HAVE_PETSC
409 : #endif // LIBMESH_STATIC_CONDENSATION_H
|
