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 : void
134 0 : CondensedEigenSystem::clear()
135 : {
136 0 : Parent::clear();
137 0 : _condensed_matrix_A = nullptr;
138 0 : _condensed_matrix_B = nullptr;
139 0 : _condensed_precond_matrix = nullptr;
140 0 : }
141 :
142 : void
143 208 : CondensedEigenSystem::add_matrices()
144 : {
145 208 : Parent::add_matrices();
146 :
147 208 : if (!this->use_shell_matrices())
148 : {
149 142 : if (!_condensed_matrix_A)
150 280 : _condensed_matrix_A = SparseMatrix<Number>::build(this->comm());
151 142 : if (!_condensed_matrix_B)
152 280 : _condensed_matrix_B = SparseMatrix<Number>::build(this->comm());
153 :
154 : // When not using shell matrices we use A for P as well so we don't need to add P
155 : }
156 : // we *are* using shell matrices but we might not be using a shell matrix for P
157 66 : else if (!this->use_shell_precond_matrix())
158 : {
159 66 : if (!_condensed_precond_matrix)
160 132 : _condensed_precond_matrix = SparseMatrix<Number>::build(this->comm());
161 : }
162 208 : }
163 :
164 219 : void CondensedEigenSystem::solve()
165 : {
166 4 : LOG_SCOPE("solve()", "CondensedEigenSystem");
167 :
168 : // If we haven't initialized any condensed dofs,
169 : // just use the default eigen_system
170 219 : if (!_have_condensed_dofs)
171 : {
172 0 : Parent::solve();
173 0 : return;
174 : }
175 :
176 : // check that necessary parameters have been set
177 4 : libmesh_assert(
178 : this->get_equation_systems().parameters.have_parameter<unsigned int>("eigenpairs"));
179 4 : libmesh_assert(
180 : this->get_equation_systems().parameters.have_parameter<unsigned int>("basis vectors"));
181 :
182 219 : if (this->assemble_before_solve)
183 : {
184 : // Assemble the linear system
185 140 : this->assemble();
186 :
187 : // And close the assembled matrices; using a non-closed matrix
188 : // with create_submatrix() is deprecated.
189 140 : if (matrix_A)
190 140 : matrix_A->close();
191 140 : if (generalized() && matrix_B)
192 140 : matrix_B->close();
193 140 : if (precond_matrix)
194 0 : precond_matrix->close();
195 : }
196 :
197 : // If we reach here, then there should be some non-condensed dofs
198 4 : libmesh_assert(!local_non_condensed_dofs_vector.empty());
199 :
200 219 : if (_create_submatrices_in_solve)
201 : {
202 153 : if (matrix_A)
203 153 : matrix_A->create_submatrix(
204 153 : *_condensed_matrix_A, local_non_condensed_dofs_vector, local_non_condensed_dofs_vector);
205 153 : if (generalized() && matrix_B)
206 153 : matrix_B->create_submatrix(
207 153 : *_condensed_matrix_B, local_non_condensed_dofs_vector, local_non_condensed_dofs_vector);
208 153 : if (precond_matrix)
209 0 : precond_matrix->create_submatrix(*_condensed_precond_matrix,
210 0 : local_non_condensed_dofs_vector,
211 0 : local_non_condensed_dofs_vector);
212 : }
213 :
214 223 : solve_helper(
215 : _condensed_matrix_A.get(), _condensed_matrix_B.get(), _condensed_precond_matrix.get());
216 : }
217 :
218 779 : std::pair<Real, Real> CondensedEigenSystem::get_eigenpair(dof_id_type i)
219 : {
220 40 : LOG_SCOPE("get_eigenpair()", "CondensedEigenSystem");
221 :
222 : // If we haven't initialized any condensed dofs,
223 : // just use the default eigen_system
224 779 : if (!_have_condensed_dofs)
225 0 : return Parent::get_eigenpair(i);
226 :
227 : // If we reach here, then there should be some non-condensed dofs
228 20 : libmesh_assert(!local_non_condensed_dofs_vector.empty());
229 :
230 : // This function assumes that condensed_solve has just been called.
231 : // If this is not the case, then we will trip an asset in get_eigenpair
232 779 : std::unique_ptr<NumericVector<Number>> temp = NumericVector<Number>::build(this->comm());
233 : const dof_id_type n_local =
234 40 : cast_int<dof_id_type>(local_non_condensed_dofs_vector.size());
235 779 : dof_id_type n = n_local;
236 779 : this->comm().sum(n);
237 :
238 779 : temp->init (n, n_local, false, PARALLEL);
239 :
240 799 : std::pair<Real, Real> eval = eigen_solver->get_eigenpair (i, *temp);
241 :
242 : // Now map temp to solution. Loop over local entries of local_non_condensed_dofs_vector
243 20 : libmesh_assert(this->comm().verify(this->solution->closed()));
244 779 : if (!this->solution->closed())
245 0 : this->solution->close();
246 779 : this->solution->zero();
247 129238 : for (auto j : make_range(n_local))
248 : {
249 138509 : const dof_id_type index = local_non_condensed_dofs_vector[j];
250 138509 : solution->set(index,(*temp)(temp->first_local_index()+j));
251 : }
252 :
253 : // Enforcing constraints requires creating a ghosted version of the solution, which requires the
254 : // solution be assembled
255 779 : solution->close();
256 779 : get_dof_map().enforce_constraints_exactly(*this);
257 779 : this->update();
258 :
259 779 : return eval;
260 739 : }
261 :
262 :
263 :
264 0 : SparseMatrix<Number> & CondensedEigenSystem::get_condensed_matrix_A()
265 : {
266 0 : libmesh_assert(_condensed_matrix_A);
267 0 : return *_condensed_matrix_A;
268 : }
269 :
270 0 : SparseMatrix<Number> & CondensedEigenSystem::get_condensed_matrix_B()
271 : {
272 0 : libmesh_assert(_condensed_matrix_B);
273 0 : return *_condensed_matrix_B;
274 : }
275 :
276 132 : SparseMatrix<Number> & CondensedEigenSystem::get_condensed_precond_matrix()
277 : {
278 0 : libmesh_assert(_condensed_precond_matrix);
279 132 : return *_condensed_precond_matrix;
280 : }
281 :
282 : void
283 3762 : CondensedEigenSystem::copy_sub_to_super(const NumericVector<Number> & sub,
284 : NumericVector<Number> & super)
285 : {
286 0 : libmesh_assert_equal_to(sub.local_size(), local_non_condensed_dofs_vector.size());
287 0 : libmesh_assert_equal_to(sub.local_size() + this->get_dof_map().n_local_constrained_dofs(),
288 : super.local_size());
289 3762 : auto super_sub_view = super.get_subvector(local_non_condensed_dofs_vector);
290 3762 : *super_sub_view = sub;
291 3762 : super.restore_subvector(std::move(super_sub_view), local_non_condensed_dofs_vector);
292 3762 : }
293 :
294 : void
295 3498 : CondensedEigenSystem::copy_super_to_sub(NumericVector<Number> & super,
296 : NumericVector<Number> & sub)
297 : {
298 0 : libmesh_assert_equal_to(sub.local_size(), local_non_condensed_dofs_vector.size());
299 0 : libmesh_assert_equal_to(sub.local_size() + this->get_dof_map().n_local_constrained_dofs(),
300 : super.local_size());
301 3498 : auto super_sub_view = super.get_subvector(local_non_condensed_dofs_vector);
302 3498 : sub = *super_sub_view;
303 3498 : super.restore_subvector(std::move(super_sub_view), local_non_condensed_dofs_vector);
304 3498 : }
305 :
306 : void
307 264 : CondensedEigenSystem::copy_super_to_sub(const SparseMatrix<Number> & super,
308 : SparseMatrix<Number> & sub)
309 : {
310 0 : parallel_object_only();
311 :
312 0 : libmesh_assert_equal_to(sub.local_m(), local_non_condensed_dofs_vector.size());
313 0 : libmesh_assert_equal_to(sub.local_m() + this->get_dof_map().n_local_constrained_dofs(),
314 : super.local_m());
315 264 : auto super_sub_view = SparseMatrix<Number>::build(super.comm());
316 264 : super.create_submatrix(
317 264 : *super_sub_view, local_non_condensed_dofs_vector, local_non_condensed_dofs_vector);
318 264 : sub = *super_sub_view;
319 264 : }
320 :
321 : } // namespace libMesh
322 :
323 :
324 : #endif // LIBMESH_HAVE_SLEPC
|
