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 :
19 :
20 : // C++ includes
21 :
22 : // Local includes
23 : #include "libmesh/linear_implicit_system.h"
24 : #include "libmesh/linear_solver.h"
25 : #include "libmesh/equation_systems.h"
26 : #include "libmesh/libmesh_call_mpi.h"
27 : #include "libmesh/numeric_vector.h" // for parameter sensitivity calcs
28 : //#include "libmesh/parameter_vector.h"
29 : #include "libmesh/sparse_matrix.h" // for get_transpose
30 : #include "libmesh/system_subset.h"
31 : #include "libmesh/static_condensation.h"
32 : #include "libmesh/static_condensation_preconditioner.h"
33 :
34 : namespace libMesh
35 : {
36 :
37 : namespace
38 : {
39 :
40 : /**
41 : * Temporarily attaches a per-call solver configuration and restores the
42 : * previously attached configuration on every exit path. This prevents a
43 : * borrowed solver from retaining configuration intended for another system,
44 : * including when solver initialization or solve() throws.
45 : */
46 : class ScopedSolverConfiguration
47 : {
48 : public:
49 209635 : ScopedSolverConfiguration(LinearSolver<Number> & solver,
50 : SolverConfiguration * temporary_configuration)
51 209635 : : _solver(solver),
52 215779 : _previous_configuration(solver.solver_configuration()),
53 209635 : _restore(temporary_configuration != nullptr)
54 : {
55 209635 : if (_restore)
56 213 : _solver.set_solver_configuration(temporary_configuration);
57 209635 : }
58 :
59 203491 : ~ScopedSolverConfiguration()
60 6144 : {
61 209633 : if (_restore)
62 213 : _solver.set_solver_configuration(_previous_configuration);
63 203491 : }
64 :
65 : private:
66 : LinearSolver<Number> & _solver;
67 : SolverConfiguration * const _previous_configuration;
68 : const bool _restore;
69 : };
70 :
71 : /**
72 : * Applies a system's subset restriction to the active solver and removes it
73 : * on every exit path. Without this guard, an exception from solve() could
74 : * leave a borrowed solver restricted and affect a later solve on another
75 : * system.
76 : */
77 : class SubsetRestrictionGuard
78 : {
79 : public:
80 6144 : explicit SubsetRestrictionGuard(LinearSolver<Number> & solver)
81 203491 : : _solver(solver), _restricted(false)
82 : {
83 6144 : }
84 :
85 203491 : ~SubsetRestrictionGuard()
86 6144 : {
87 209633 : if (_restricted)
88 71 : _solver.restrict_solve_to(nullptr);
89 203491 : }
90 :
91 8 : void restrict(const std::vector<unsigned int> * dofs, const SubsetSolveMode subset_solve_mode)
92 : {
93 273 : _restricted = true;
94 281 : _solver.restrict_solve_to(dofs, subset_solve_mode);
95 273 : }
96 :
97 203422 : void clear()
98 : {
99 209564 : if (_restricted)
100 : {
101 210 : _solver.restrict_solve_to(nullptr);
102 204 : _restricted = false;
103 : }
104 203422 : }
105 :
106 : private:
107 : LinearSolver<Number> & _solver;
108 : bool _restricted;
109 : };
110 :
111 : } // anonymous namespace
112 :
113 11902 : LinearImplicitSystem::LinearImplicitSystem (EquationSystems & es,
114 : const std::string & name_in,
115 11902 : const unsigned int number_in) :
116 :
117 : Parent (es, name_in, number_in),
118 11238 : _n_linear_iterations (0),
119 11238 : _final_linear_residual (1.e20),
120 11238 : _shell_matrix(nullptr),
121 11238 : _subset(nullptr),
122 11902 : _subset_solve_mode(SUBSET_ZERO)
123 : {
124 : // linear_solver is now in the ImplicitSystem base class, but we are
125 : // going to keep using it basically the way we did before it was
126 : // moved.
127 23140 : linear_solver = LinearSolver<Number>::build(es.comm());
128 :
129 11902 : if (this->has_static_condensation())
130 140 : this->setup_static_condensation_preconditioner(*linear_solver);
131 11902 : }
132 :
133 :
134 :
135 21135 : LinearImplicitSystem::~LinearImplicitSystem () = default;
136 :
137 :
138 :
139 140 : void LinearImplicitSystem::create_static_condensation()
140 : {
141 140 : Parent::create_static_condensation();
142 140 : this->setup_static_condensation_preconditioner(*linear_solver);
143 140 : }
144 :
145 :
146 :
147 217 : void LinearImplicitSystem::clear ()
148 : {
149 : // clear the linear solver
150 209 : linear_solver->clear();
151 :
152 217 : this->restrict_solve_to(nullptr);
153 :
154 : // clear the parent data
155 217 : Parent::clear();
156 :
157 : // And restore any StaticCondensation to defaults
158 217 : if (this->has_static_condensation())
159 70 : this->setup_static_condensation_preconditioner(*linear_solver);
160 217 : }
161 :
162 :
163 :
164 11902 : void LinearImplicitSystem::init_data ()
165 : {
166 : // initialize parent data
167 11902 : Parent::init_data();
168 :
169 : // re-initialize the linear solver interface
170 11902 : linear_solver->clear();
171 11902 : }
172 :
173 :
174 :
175 18036 : void LinearImplicitSystem::reinit ()
176 : {
177 : // re-initialize the linear solver interface
178 18036 : linear_solver->clear();
179 :
180 : // initialize parent data
181 18036 : Parent::reinit();
182 18036 : }
183 :
184 :
185 :
186 569 : void LinearImplicitSystem::restrict_solve_to (const SystemSubset * subset,
187 : const SubsetSolveMode subset_solve_mode)
188 : {
189 569 : _subset = subset;
190 569 : _subset_solve_mode = subset_solve_mode;
191 :
192 18 : if (subset != nullptr)
193 8 : libmesh_assert_equal_to (&subset->get_system(), this);
194 569 : }
195 :
196 209351 : void LinearImplicitSystem::solve ()
197 : {
198 209351 : this->solve(LinearImplicitSystemSolveOptions{});
199 209351 : }
200 :
201 209635 : void LinearImplicitSystem::solve(const LinearImplicitSystemSolveOptions & options)
202 : {
203 209635 : LinearSolver<Number> * const solver = options.solver ? options.solver : linear_solver.get();
204 6144 : libmesh_assert(solver);
205 :
206 : #ifndef NDEBUG
207 : #ifdef LIBMESH_HAVE_MPI
208 : int communicator_comparison;
209 6144 : libmesh_call_mpi(
210 : MPI_Comm_compare(solver->comm().get(), this->comm().get(), &communicator_comparison));
211 6144 : libmesh_assert(communicator_comparison == MPI_IDENT || communicator_comparison == MPI_CONGRUENT);
212 : #else
213 : libmesh_assert_equal_to(solver->comm().get(), this->comm().get());
214 : #endif
215 : #endif
216 :
217 6144 : const bool assemble = options.assemble_before_solve.value_or(this->assemble_before_solve);
218 :
219 209635 : if (assemble)
220 : // Assemble the linear system
221 67280 : this->assemble();
222 :
223 215779 : ScopedSolverConfiguration scoped_configuration(*solver, options.solver_configuration);
224 :
225 : // If the linear solver hasn't been initialized, we do so here.
226 209635 : if (this->prefix_with_name())
227 2 : solver->init(this->prefix().c_str());
228 : else
229 209635 : solver->init();
230 :
231 209635 : solver->init_systems(*this);
232 :
233 12288 : SubsetRestrictionGuard subset_guard(*solver);
234 :
235 209635 : if (_subset != nullptr)
236 281 : subset_guard.restrict(&_subset->dof_ids(), _subset_solve_mode);
237 :
238 : // Solve the linear system. Several cases:
239 6144 : std::pair<unsigned int, Real> rval = std::make_pair(0,0.0);
240 209635 : SparseMatrix<Number> * const preconditioner = this->request_matrix("Preconditioner");
241 :
242 209635 : if (options.solver_configuration)
243 : {
244 213 : if (_shell_matrix)
245 : // 1.) Shell matrix with or without user-supplied preconditioner.
246 0 : rval = solver->solve(*_shell_matrix, preconditioner, *solution, *rhs);
247 : else
248 : // 2.) No shell matrix, with or without user-supplied preconditioner.
249 219 : rval = solver->solve(*matrix, preconditioner, *solution, *rhs);
250 : }
251 : else
252 : {
253 : // Get the user-specified linear solver tolerance and iteration limit.
254 209422 : const auto [maxits, tol] = this->get_linear_solve_parameters();
255 :
256 209422 : if (_shell_matrix)
257 : // 1.) Shell matrix with or without user-supplied preconditioner.
258 72 : rval = solver->solve(*_shell_matrix, preconditioner, *solution, *rhs, tol, maxits);
259 : else
260 : // 2.) No shell matrix, with or without user-supplied preconditioner.
261 215490 : rval = solver->solve(*matrix, preconditioner, *solution, *rhs, tol, maxits);
262 : }
263 :
264 203422 : subset_guard.clear();
265 :
266 : // Store the number of linear iterations required to
267 : // solve and the final residual.
268 209564 : _n_linear_iterations = rval.first;
269 209564 : _final_linear_residual = rval.second;
270 :
271 : // Update the system after the solve
272 209564 : this->update();
273 209698 : }
274 :
275 140 : void LinearImplicitSystem::attach_shell_matrix (ShellMatrix<Number> * shell_matrix)
276 : {
277 140 : _shell_matrix = shell_matrix;
278 140 : }
279 :
280 :
281 : /*
282 : void LinearImplicitSystem::sensitivity_solve (const ParameterVector & parameters)
283 : {
284 : if (this->assemble_before_solve)
285 : {
286 : // Assemble the linear system
287 : this->assemble ();
288 :
289 : // But now assemble right hand sides with the residual's
290 : // parameter derivatives
291 : this->assemble_residual_derivatives(parameters);
292 : }
293 :
294 : // Get a reference to the EquationSystems
295 : const EquationSystems & es =
296 : this->get_equation_systems();
297 :
298 : // Get the user-specified linear solver tolerance
299 : const Real tol =
300 : es.parameters.get<Real>("sensitivity solver tolerance");
301 :
302 : // Get the user-specified maximum # of linear solver iterations
303 : const unsigned int maxits =
304 : es.parameters.get<unsigned int>("sensitivity solver maximum iterations");
305 :
306 : // Our iteration counts and residuals will be sums of the individual
307 : // results
308 : _n_linear_iterations = 0;
309 : _final_linear_residual = 0.0;
310 : std::pair<unsigned int, Real> rval = std::make_pair(0,0.0);
311 :
312 : // Solve the linear system.
313 : SparseMatrix<Number> * pc = this->request_matrix("Preconditioner");
314 : for (std::size_t p=0; p != parameters.size(); ++p)
315 : {
316 : rval = linear_solver->solve (*matrix, pc,
317 : this->get_sensitivity_solution(p),
318 : this->get_sensitivity_rhs(p), tol, maxits);
319 : _n_linear_iterations += rval.first;
320 : _final_linear_residual += rval.second;
321 : }
322 :
323 : // Our matrix is the *negative* of the Jacobian for b-A*u, so our
324 : // solutions are all inverted
325 : for (std::size_t p=0; p != parameters.size(); ++p)
326 : {
327 : this->get_sensitivity_solution(p) *= -1.0;
328 : }
329 : }
330 :
331 :
332 :
333 : void LinearImplicitSystem::adjoint_solve (const QoISet & qoi_indices)
334 : {
335 : const unsigned int Nq = this->n_qois();
336 :
337 : // We currently don't support adjoint solves of shell matrices
338 : // FIXME - we should let shell matrices support
339 : // vector_transpose_mult so that we can use them here.
340 : if (_shell_matrix!=nullptr)
341 : libmesh_not_implemented();
342 :
343 : if (this->assemble_before_solve)
344 : {
345 : // Assemble the linear system
346 : this->assemble ();
347 :
348 : // And take the adjoint
349 : matrix->get_transpose(*matrix);
350 :
351 : // Including of any separate preconditioner
352 : SparseMatrix<Number> * pc = this->request_matrix("Preconditioner");
353 : if (pc)
354 : pc->get_transpose(*pc);
355 :
356 : // But now replace the right hand sides with the quantity of
357 : // interest functional's derivative
358 : this->assemble_qoi_derivative(qoi_indices);
359 : }
360 :
361 : // Get a reference to the EquationSystems
362 : const EquationSystems & es =
363 : this->get_equation_systems();
364 :
365 : // Get the user-specified linear solver tolerance
366 : const Real tol =
367 : es.parameters.get<Real>("adjoint solver tolerance");
368 :
369 : // Get the user-specified maximum # of linear solver iterations
370 : const unsigned int maxits =
371 : es.parameters.get<unsigned int>("adjoint solver maximum iterations");
372 :
373 : // Our iteration counts and residuals will be sums of the individual
374 : // results
375 : _n_linear_iterations = 0;
376 : _final_linear_residual = 0.0;
377 : std::pair<unsigned int, Real> rval = std::make_pair(0,0.0);
378 :
379 : // Solve the linear system.
380 : SparseMatrix<Number> * pc = this->request_matrix("Preconditioner");
381 : for (unsigned int i=0; i != Nq; ++i)
382 : if (qoi_indices.has_index(i))
383 : {
384 : rval = linear_solver->solve (*matrix, pc,
385 : this->add_adjoint_solution(i),
386 : this->get_adjoint_rhs(i), tol, maxits);
387 : _n_linear_iterations += rval.first;
388 : _final_linear_residual += rval.second;
389 : }
390 : }
391 :
392 :
393 :
394 : void LinearImplicitSystem::forward_qoi_parameter_sensitivity
395 : (const QoISet & qoi_indices,
396 : const ParameterVector & parameters,
397 : SensitivityData & sensitivities)
398 : {
399 : const unsigned int Np = parameters.size();
400 : const unsigned int Nq = this->n_qois();
401 :
402 : // An introduction to the problem:
403 : //
404 : // A(p)*u(p) = b(p), where x is determined implicitly.
405 : // Residual R(u(p),p) := b(p) - A(p)*u(p)
406 : // partial R / partial u = -A
407 : //
408 : // This implies that:
409 : // d/dp(R) = 0
410 : // (partial b / partial p) -
411 : // (partial A / partial p) * u -
412 : // A * (partial u / partial p) = 0
413 : // A * (partial u / partial p) = (partial R / partial p)
414 : // = (partial b / partial p) - (partial A / partial p) * u
415 :
416 : // We first solve for (partial u / partial p) for each parameter:
417 : // -A * (partial u / partial p) = - (partial R / partial p)
418 :
419 : this->sensitivity_solve(parameters);
420 :
421 : // Get ready to fill in sensitivities:
422 : sensitivities.allocate_data(qoi_indices, *this, parameters);
423 :
424 : // We use the identity:
425 : // dq/dp = (partial q / partial p) + (partial q / partial u) *
426 : // (partial u / partial p)
427 :
428 : // We get (partial q / partial u) from the user
429 : this->assemble_qoi_derivative(qoi_indices);
430 :
431 : for (unsigned int j=0; j != Np; ++j)
432 : {
433 : // We currently get partial derivatives via central differencing
434 : Number delta_p = 1e-6;
435 :
436 : // (partial q / partial p) ~= (q(p+dp)-q(p-dp))/(2*dp)
437 :
438 : Number old_parameter = *parameters[j];
439 :
440 : *parameters[j] = old_parameter - delta_p;
441 : this->assemble_qoi(qoi_indices);
442 : std::vector<Number> qoi_minus = this->qoi;
443 :
444 : *parameters[j] = old_parameter + delta_p;
445 : this->assemble_qoi(qoi_indices);
446 : std::vector<Number> & qoi_plus = this->qoi;
447 : std::vector<Number> partialq_partialp(Nq, 0);
448 : for (unsigned int i=0; i != Nq; ++i)
449 : if (qoi_indices.has_index(i))
450 : partialq_partialp[i] = (qoi_plus[i] - qoi_minus[i]) / (2.*delta_p);
451 :
452 : for (unsigned int i=0; i != Nq; ++i)
453 : if (qoi_indices.has_index(i))
454 : sensitivities[i][j] = partialq_partialp[i] +
455 : this->get_adjoint_rhs(i).dot(this->get_sensitivity_solution(i));
456 : }
457 :
458 : // All parameters have been reset.
459 : // Don't leave the qoi or system changed - principle of least
460 : // surprise.
461 : this->assemble();
462 : this->rhs->close();
463 : this->matrix->close();
464 : this->assemble_qoi(qoi_indices);
465 : }
466 : */
467 :
468 :
469 :
470 289244 : LinearSolver<Number> * LinearImplicitSystem::get_linear_solver() const
471 : {
472 289244 : return linear_solver.get();
473 : }
474 :
475 :
476 :
477 0 : void LinearImplicitSystem::assembly(bool,
478 : bool,
479 : bool,
480 : bool)
481 : {
482 : // Residual R(u(p),p) := A(p)*u(p) - b(p)
483 : // partial R / partial u = A
484 :
485 0 : this->assemble();
486 0 : this->rhs->close();
487 0 : this->matrix->close();
488 :
489 0 : *(this->rhs) *= -1.0;
490 0 : this->rhs->add_vector(*this->solution, *this->matrix);
491 0 : }
492 :
493 : } // namespace libMesh
|
