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 : #include "libmesh/libmesh_config.h"
19 :
20 : // Currently, the EigenSystem should only be available
21 : // if SLEPc support is enabled.
22 : #if defined(LIBMESH_HAVE_SLEPC)
23 :
24 : #include "libmesh/condensed_eigen_system.h"
25 :
26 : #include "libmesh/dof_map.h"
27 : #include "libmesh/equation_systems.h"
28 : #include "libmesh/int_range.h"
29 : #include "libmesh/libmesh_logging.h"
30 : #include "libmesh/numeric_vector.h"
31 : #include "libmesh/parallel.h"
32 :
33 : namespace libMesh
34 : {
35 :
36 208 : CondensedEigenSystem::CondensedEigenSystem (EquationSystems & es,
37 : const std::string & name_in,
38 208 : const unsigned int number_in)
39 : : Parent(es, name_in, number_in),
40 200 : _condensed_dofs_initialized(false),
41 200 : _have_condensed_dofs(false),
42 208 : _create_submatrices_in_solve(true)
43 : {
44 208 : }
45 :
46 598 : CondensedEigenSystem::~CondensedEigenSystem() = default;
47 :
48 : void
49 207 : CondensedEigenSystem::initialize_condensed_dofs(const std::set<dof_id_type> & global_dirichlet_dofs_set)
50 : {
51 4 : const DofMap & dof_map = this->get_dof_map();
52 207 : if (global_dirichlet_dofs_set.empty() && !dof_map.n_constrained_dofs())
53 : {
54 0 : _have_condensed_dofs = false;
55 0 : _condensed_dofs_initialized = true;
56 0 : return;
57 : }
58 :
59 : // First, put all unconstrained local dofs into non_dirichlet_dofs_set
60 4 : std::set<dof_id_type> local_non_condensed_dofs_set;
61 36328 : for (auto i : make_range(dof_map.first_dof(), dof_map.end_dof()))
62 : #if LIBMESH_ENABLE_CONSTRAINTS
63 5550 : if (!dof_map.is_constrained_dof(i))
64 : #endif
65 30235 : local_non_condensed_dofs_set.insert(i);
66 :
67 : // Now erase the condensed dofs
68 233 : for (const auto dof : global_dirichlet_dofs_set)
69 26 : if ((dof_map.first_dof() <= dof) && (dof < dof_map.end_dof()))
70 0 : local_non_condensed_dofs_set.erase(dof);
71 :
72 : // Finally, move local_non_condensed_dofs_set over to a vector for convenience in solve()
73 207 : this->local_non_condensed_dofs_vector.clear();
74 :
75 32426 : for (const auto dof : local_non_condensed_dofs_set)
76 32219 : this->local_non_condensed_dofs_vector.push_back(dof);
77 :
78 207 : _have_condensed_dofs = true;
79 207 : _condensed_dofs_initialized = true;
80 : }
81 :
82 : void
83 0 : CondensedEigenSystem::reinit()
84 : {
85 0 : Parent::reinit();
86 0 : _condensed_dofs_initialized = false;
87 0 : }
88 :
89 : void
90 66 : CondensedEigenSystem::initialize_condensed_matrices()
91 : {
92 66 : if (!_condensed_dofs_initialized)
93 132 : this->initialize_condensed_dofs();
94 :
95 66 : if (!_have_condensed_dofs)
96 0 : return;
97 :
98 0 : const auto m = cast_int<numeric_index_type>(this->local_non_condensed_dofs_vector.size());
99 66 : auto M = m;
100 66 : _communicator.sum(M);
101 66 : if (this->has_condensed_matrix_A())
102 : {
103 0 : this->get_condensed_matrix_A().init(M, M, m, m);
104 : // A bit ludicrously MatCopy requires the matrix being copied to to be assembled
105 0 : this->get_condensed_matrix_A().close();
106 : }
107 66 : if (this->has_condensed_matrix_B())
108 : {
109 0 : this->get_condensed_matrix_B().init(M, M, m, m);
110 0 : this->get_condensed_matrix_B().close();
111 : }
112 66 : if (this->has_condensed_precond_matrix())
113 : {
114 66 : this->get_condensed_precond_matrix().init(M, M, m, m);
115 66 : this->get_condensed_precond_matrix().close();
116 : }
117 : }
118 :
119 0 : dof_id_type CondensedEigenSystem::n_global_non_condensed_dofs() const
120 : {
121 0 : if (!_have_condensed_dofs)
122 0 : return this->n_dofs();
123 : else
124 : {
125 : dof_id_type n_global_non_condensed_dofs =
126 0 : cast_int<dof_id_type>(local_non_condensed_dofs_vector.size());
127 0 : this->comm().sum(n_global_non_condensed_dofs);
128 :
129 0 : return n_global_non_condensed_dofs;
130 : }
131 : }
132 :
133 : #ifdef LIBMESH_ENABLE_DEPRECATED
134 : void
135 208 : CondensedEigenSystem::set_raw_pointers()
136 : {
137 208 : condensed_matrix_A = _condensed_matrix_A.get();
138 208 : condensed_matrix_B = _condensed_matrix_B.get();
139 208 : }
140 : #endif
141 :
142 : void
143 0 : CondensedEigenSystem::clear()
144 : {
145 0 : Parent::clear();
146 0 : _condensed_matrix_A = nullptr;
147 0 : _condensed_matrix_B = nullptr;
148 0 : _condensed_precond_matrix = nullptr;
149 : #ifdef LIBMESH_ENABLE_DEPRECATED
150 0 : set_raw_pointers();
151 : #endif
152 0 : }
153 :
154 : void
155 208 : CondensedEigenSystem::add_matrices()
156 : {
157 208 : Parent::add_matrices();
158 :
159 208 : if (!this->use_shell_matrices())
160 : {
161 142 : if (!_condensed_matrix_A)
162 280 : _condensed_matrix_A = SparseMatrix<Number>::build(this->comm());
163 142 : if (!_condensed_matrix_B)
164 280 : _condensed_matrix_B = SparseMatrix<Number>::build(this->comm());
165 :
166 : // When not using shell matrices we use A for P as well so we don't need to add P
167 : }
168 : // we *are* using shell matrices but we might not be using a shell matrix for P
169 66 : else if (!this->use_shell_precond_matrix())
170 : {
171 66 : if (!_condensed_precond_matrix)
172 132 : _condensed_precond_matrix = SparseMatrix<Number>::build(this->comm());
173 : }
174 :
175 : #ifdef LIBMESH_ENABLE_DEPRECATED
176 208 : set_raw_pointers();
177 : #endif
178 208 : }
179 :
180 219 : void CondensedEigenSystem::solve()
181 : {
182 4 : LOG_SCOPE("solve()", "CondensedEigenSystem");
183 :
184 : // If we haven't initialized any condensed dofs,
185 : // just use the default eigen_system
186 219 : if (!_have_condensed_dofs)
187 : {
188 0 : Parent::solve();
189 0 : return;
190 : }
191 :
192 : // check that necessary parameters have been set
193 4 : libmesh_assert(
194 : this->get_equation_systems().parameters.have_parameter<unsigned int>("eigenpairs"));
195 4 : libmesh_assert(
196 : this->get_equation_systems().parameters.have_parameter<unsigned int>("basis vectors"));
197 :
198 219 : if (this->assemble_before_solve)
199 : {
200 : // Assemble the linear system
201 140 : this->assemble();
202 :
203 : // And close the assembled matrices; using a non-closed matrix
204 : // with create_submatrix() is deprecated.
205 140 : if (matrix_A)
206 140 : matrix_A->close();
207 140 : if (generalized() && matrix_B)
208 140 : matrix_B->close();
209 140 : if (precond_matrix)
210 0 : precond_matrix->close();
211 : }
212 :
213 : // If we reach here, then there should be some non-condensed dofs
214 4 : libmesh_assert(!local_non_condensed_dofs_vector.empty());
215 :
216 219 : if (_create_submatrices_in_solve)
217 : {
218 153 : if (matrix_A)
219 153 : matrix_A->create_submatrix(
220 153 : *_condensed_matrix_A, local_non_condensed_dofs_vector, local_non_condensed_dofs_vector);
221 153 : if (generalized() && matrix_B)
222 153 : matrix_B->create_submatrix(
223 153 : *_condensed_matrix_B, local_non_condensed_dofs_vector, local_non_condensed_dofs_vector);
224 153 : if (precond_matrix)
225 0 : precond_matrix->create_submatrix(*_condensed_precond_matrix,
226 0 : local_non_condensed_dofs_vector,
227 0 : local_non_condensed_dofs_vector);
228 : }
229 :
230 223 : solve_helper(
231 : _condensed_matrix_A.get(), _condensed_matrix_B.get(), _condensed_precond_matrix.get());
232 : }
233 :
234 779 : std::pair<Real, Real> CondensedEigenSystem::get_eigenpair(dof_id_type i)
235 : {
236 40 : LOG_SCOPE("get_eigenpair()", "CondensedEigenSystem");
237 :
238 : // If we haven't initialized any condensed dofs,
239 : // just use the default eigen_system
240 779 : if (!_have_condensed_dofs)
241 0 : return Parent::get_eigenpair(i);
242 :
243 : // If we reach here, then there should be some non-condensed dofs
244 20 : libmesh_assert(!local_non_condensed_dofs_vector.empty());
245 :
246 : // This function assumes that condensed_solve has just been called.
247 : // If this is not the case, then we will trip an asset in get_eigenpair
248 779 : std::unique_ptr<NumericVector<Number>> temp = NumericVector<Number>::build(this->comm());
249 : const dof_id_type n_local =
250 40 : cast_int<dof_id_type>(local_non_condensed_dofs_vector.size());
251 779 : dof_id_type n = n_local;
252 779 : this->comm().sum(n);
253 :
254 779 : temp->init (n, n_local, false, PARALLEL);
255 :
256 799 : std::pair<Real, Real> eval = eigen_solver->get_eigenpair (i, *temp);
257 :
258 : // Now map temp to solution. Loop over local entries of local_non_condensed_dofs_vector
259 20 : libmesh_assert(this->comm().verify(this->solution->closed()));
260 779 : if (!this->solution->closed())
261 0 : this->solution->close();
262 779 : this->solution->zero();
263 129238 : for (auto j : make_range(n_local))
264 : {
265 138509 : const dof_id_type index = local_non_condensed_dofs_vector[j];
266 138509 : solution->set(index,(*temp)(temp->first_local_index()+j));
267 : }
268 :
269 : // Enforcing constraints requires creating a ghosted version of the solution, which requires the
270 : // solution be assembled
271 779 : solution->close();
272 779 : get_dof_map().enforce_constraints_exactly(*this);
273 779 : this->update();
274 :
275 779 : return eval;
276 739 : }
277 :
278 :
279 :
280 0 : SparseMatrix<Number> & CondensedEigenSystem::get_condensed_matrix_A()
281 : {
282 0 : libmesh_assert(_condensed_matrix_A);
283 0 : return *_condensed_matrix_A;
284 : }
285 :
286 0 : SparseMatrix<Number> & CondensedEigenSystem::get_condensed_matrix_B()
287 : {
288 0 : libmesh_assert(_condensed_matrix_B);
289 0 : return *_condensed_matrix_B;
290 : }
291 :
292 132 : SparseMatrix<Number> & CondensedEigenSystem::get_condensed_precond_matrix()
293 : {
294 0 : libmesh_assert(_condensed_precond_matrix);
295 132 : return *_condensed_precond_matrix;
296 : }
297 :
298 : void
299 3762 : CondensedEigenSystem::copy_sub_to_super(const NumericVector<Number> & sub,
300 : NumericVector<Number> & super)
301 : {
302 0 : libmesh_assert_equal_to(sub.local_size(), local_non_condensed_dofs_vector.size());
303 0 : libmesh_assert_equal_to(sub.local_size() + this->get_dof_map().n_local_constrained_dofs(),
304 : super.local_size());
305 3762 : auto super_sub_view = super.get_subvector(local_non_condensed_dofs_vector);
306 3762 : *super_sub_view = sub;
307 3762 : super.restore_subvector(std::move(super_sub_view), local_non_condensed_dofs_vector);
308 3762 : }
309 :
310 : void
311 3498 : CondensedEigenSystem::copy_super_to_sub(NumericVector<Number> & super,
312 : NumericVector<Number> & sub)
313 : {
314 0 : libmesh_assert_equal_to(sub.local_size(), local_non_condensed_dofs_vector.size());
315 0 : libmesh_assert_equal_to(sub.local_size() + this->get_dof_map().n_local_constrained_dofs(),
316 : super.local_size());
317 3498 : auto super_sub_view = super.get_subvector(local_non_condensed_dofs_vector);
318 3498 : sub = *super_sub_view;
319 3498 : super.restore_subvector(std::move(super_sub_view), local_non_condensed_dofs_vector);
320 3498 : }
321 :
322 : void
323 264 : CondensedEigenSystem::copy_super_to_sub(const SparseMatrix<Number> & super,
324 : SparseMatrix<Number> & sub)
325 : {
326 0 : parallel_object_only();
327 :
328 0 : libmesh_assert_equal_to(sub.local_m(), local_non_condensed_dofs_vector.size());
329 0 : libmesh_assert_equal_to(sub.local_m() + this->get_dof_map().n_local_constrained_dofs(),
330 : super.local_m());
331 264 : auto super_sub_view = SparseMatrix<Number>::build(super.comm());
332 264 : super.create_submatrix(
333 264 : *super_sub_view, local_non_condensed_dofs_vector, local_non_condensed_dofs_vector);
334 264 : sub = *super_sub_view;
335 264 : }
336 :
337 : } // namespace libMesh
338 :
339 :
340 : #endif // LIBMESH_HAVE_SLEPC
|
