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/numeric_vector.h" // for parameter sensitivity calcs
27 : //#include "libmesh/parameter_vector.h"
28 : #include "libmesh/sparse_matrix.h" // for get_transpose
29 : #include "libmesh/system_subset.h"
30 : #include "libmesh/static_condensation.h"
31 : #include "libmesh/static_condensation_preconditioner.h"
32 :
33 : namespace libMesh
34 : {
35 :
36 :
37 11689 : LinearImplicitSystem::LinearImplicitSystem (EquationSystems & es,
38 : const std::string & name_in,
39 11689 : const unsigned int number_in) :
40 :
41 : Parent (es, name_in, number_in),
42 11037 : _n_linear_iterations (0),
43 11037 : _final_linear_residual (1.e20),
44 11037 : _shell_matrix(nullptr),
45 11037 : _subset(nullptr),
46 11689 : _subset_solve_mode(SUBSET_ZERO)
47 : {
48 : // linear_solver is now in the ImplicitSystem base class, but we are
49 : // going to keep using it basically the way we did before it was
50 : // moved.
51 22726 : linear_solver = LinearSolver<Number>::build(es.comm());
52 :
53 11689 : if (this->has_static_condensation())
54 140 : this->setup_static_condensation_preconditioner(*linear_solver);
55 11689 : }
56 :
57 :
58 :
59 20733 : LinearImplicitSystem::~LinearImplicitSystem () = default;
60 :
61 :
62 :
63 140 : void LinearImplicitSystem::create_static_condensation()
64 : {
65 140 : Parent::create_static_condensation();
66 140 : this->setup_static_condensation_preconditioner(*linear_solver);
67 140 : }
68 :
69 :
70 :
71 217 : void LinearImplicitSystem::clear ()
72 : {
73 : // clear the linear solver
74 217 : linear_solver->clear();
75 :
76 217 : this->restrict_solve_to(nullptr);
77 :
78 : // clear the parent data
79 217 : Parent::clear();
80 :
81 : // And restore any StaticCondensation to defaults
82 217 : if (this->has_static_condensation())
83 70 : this->setup_static_condensation_preconditioner(*linear_solver);
84 217 : }
85 :
86 :
87 :
88 11689 : void LinearImplicitSystem::init_data ()
89 : {
90 : // initialize parent data
91 11689 : Parent::init_data();
92 :
93 : // re-initialize the linear solver interface
94 11689 : linear_solver->clear();
95 11689 : }
96 :
97 :
98 :
99 18036 : void LinearImplicitSystem::reinit ()
100 : {
101 : // re-initialize the linear solver interface
102 18036 : linear_solver->clear();
103 :
104 : // initialize parent data
105 18036 : Parent::reinit();
106 18036 : }
107 :
108 :
109 :
110 427 : void LinearImplicitSystem::restrict_solve_to (const SystemSubset * subset,
111 : const SubsetSolveMode subset_solve_mode)
112 : {
113 427 : _subset = subset;
114 427 : _subset_solve_mode = subset_solve_mode;
115 :
116 14 : if (subset != nullptr)
117 6 : libmesh_assert_equal_to (&subset->get_system(), this);
118 427 : }
119 :
120 :
121 :
122 209280 : void LinearImplicitSystem::solve ()
123 : {
124 209280 : if (this->assemble_before_solve)
125 : // Assemble the linear system
126 67209 : this->assemble ();
127 :
128 : // If the linear solver hasn't been initialized, we do so here.
129 209280 : if (this->prefix_with_name())
130 0 : linear_solver->init(this->prefix().c_str());
131 : else
132 209280 : linear_solver->init();
133 :
134 209280 : linear_solver->init_systems(*this);
135 :
136 : // Get the user-specified linear solver tolerance
137 209280 : const auto [maxits, tol] = this->get_linear_solve_parameters();
138 :
139 209280 : if (_subset != nullptr)
140 210 : linear_solver->restrict_solve_to(&_subset->dof_ids(),_subset_solve_mode);
141 :
142 : // Solve the linear system. Several cases:
143 6134 : std::pair<unsigned int, Real> rval = std::make_pair(0,0.0);
144 209280 : if (_shell_matrix)
145 : // 1.) Shell matrix with or without user-supplied preconditioner.
146 70 : rval = linear_solver->solve(*_shell_matrix, this->request_matrix("Preconditioner"), *solution, *rhs, tol, maxits);
147 : else
148 : // 2.) No shell matrix, with or without user-supplied preconditioner
149 209210 : rval = linear_solver->solve (*matrix, this->request_matrix("Preconditioner"), *solution, *rhs, tol, maxits);
150 :
151 209280 : if (_subset != nullptr)
152 210 : linear_solver->restrict_solve_to(nullptr);
153 :
154 : // Store the number of linear iterations required to
155 : // solve and the final residual.
156 209280 : _n_linear_iterations = rval.first;
157 209280 : _final_linear_residual = rval.second;
158 :
159 : // Update the system after the solve
160 209280 : this->update();
161 209280 : }
162 :
163 :
164 :
165 140 : void LinearImplicitSystem::attach_shell_matrix (ShellMatrix<Number> * shell_matrix)
166 : {
167 140 : _shell_matrix = shell_matrix;
168 140 : }
169 :
170 :
171 : /*
172 : void LinearImplicitSystem::sensitivity_solve (const ParameterVector & parameters)
173 : {
174 : if (this->assemble_before_solve)
175 : {
176 : // Assemble the linear system
177 : this->assemble ();
178 :
179 : // But now assemble right hand sides with the residual's
180 : // parameter derivatives
181 : this->assemble_residual_derivatives(parameters);
182 : }
183 :
184 : // Get a reference to the EquationSystems
185 : const EquationSystems & es =
186 : this->get_equation_systems();
187 :
188 : // Get the user-specified linear solver tolerance
189 : const Real tol =
190 : es.parameters.get<Real>("sensitivity solver tolerance");
191 :
192 : // Get the user-specified maximum # of linear solver iterations
193 : const unsigned int maxits =
194 : es.parameters.get<unsigned int>("sensitivity solver maximum iterations");
195 :
196 : // Our iteration counts and residuals will be sums of the individual
197 : // results
198 : _n_linear_iterations = 0;
199 : _final_linear_residual = 0.0;
200 : std::pair<unsigned int, Real> rval = std::make_pair(0,0.0);
201 :
202 : // Solve the linear system.
203 : SparseMatrix<Number> * pc = this->request_matrix("Preconditioner");
204 : for (std::size_t p=0; p != parameters.size(); ++p)
205 : {
206 : rval = linear_solver->solve (*matrix, pc,
207 : this->get_sensitivity_solution(p),
208 : this->get_sensitivity_rhs(p), tol, maxits);
209 : _n_linear_iterations += rval.first;
210 : _final_linear_residual += rval.second;
211 : }
212 :
213 : // Our matrix is the *negative* of the Jacobian for b-A*u, so our
214 : // solutions are all inverted
215 : for (std::size_t p=0; p != parameters.size(); ++p)
216 : {
217 : this->get_sensitivity_solution(p) *= -1.0;
218 : }
219 : }
220 :
221 :
222 :
223 : void LinearImplicitSystem::adjoint_solve (const QoISet & qoi_indices)
224 : {
225 : const unsigned int Nq = this->n_qois();
226 :
227 : // We currently don't support adjoint solves of shell matrices
228 : // FIXME - we should let shell matrices support
229 : // vector_transpose_mult so that we can use them here.
230 : if (_shell_matrix!=nullptr)
231 : libmesh_not_implemented();
232 :
233 : if (this->assemble_before_solve)
234 : {
235 : // Assemble the linear system
236 : this->assemble ();
237 :
238 : // And take the adjoint
239 : matrix->get_transpose(*matrix);
240 :
241 : // Including of any separate preconditioner
242 : SparseMatrix<Number> * pc = this->request_matrix("Preconditioner");
243 : if (pc)
244 : pc->get_transpose(*pc);
245 :
246 : // But now replace the right hand sides with the quantity of
247 : // interest functional's derivative
248 : this->assemble_qoi_derivative(qoi_indices);
249 : }
250 :
251 : // Get a reference to the EquationSystems
252 : const EquationSystems & es =
253 : this->get_equation_systems();
254 :
255 : // Get the user-specified linear solver tolerance
256 : const Real tol =
257 : es.parameters.get<Real>("adjoint solver tolerance");
258 :
259 : // Get the user-specified maximum # of linear solver iterations
260 : const unsigned int maxits =
261 : es.parameters.get<unsigned int>("adjoint solver maximum iterations");
262 :
263 : // Our iteration counts and residuals will be sums of the individual
264 : // results
265 : _n_linear_iterations = 0;
266 : _final_linear_residual = 0.0;
267 : std::pair<unsigned int, Real> rval = std::make_pair(0,0.0);
268 :
269 : // Solve the linear system.
270 : SparseMatrix<Number> * pc = this->request_matrix("Preconditioner");
271 : for (unsigned int i=0; i != Nq; ++i)
272 : if (qoi_indices.has_index(i))
273 : {
274 : rval = linear_solver->solve (*matrix, pc,
275 : this->add_adjoint_solution(i),
276 : this->get_adjoint_rhs(i), tol, maxits);
277 : _n_linear_iterations += rval.first;
278 : _final_linear_residual += rval.second;
279 : }
280 : }
281 :
282 :
283 :
284 : void LinearImplicitSystem::forward_qoi_parameter_sensitivity
285 : (const QoISet & qoi_indices,
286 : const ParameterVector & parameters,
287 : SensitivityData & sensitivities)
288 : {
289 : const unsigned int Np = parameters.size();
290 : const unsigned int Nq = this->n_qois();
291 :
292 : // An introduction to the problem:
293 : //
294 : // A(p)*u(p) = b(p), where x is determined implicitly.
295 : // Residual R(u(p),p) := b(p) - A(p)*u(p)
296 : // partial R / partial u = -A
297 : //
298 : // This implies that:
299 : // d/dp(R) = 0
300 : // (partial b / partial p) -
301 : // (partial A / partial p) * u -
302 : // A * (partial u / partial p) = 0
303 : // A * (partial u / partial p) = (partial R / partial p)
304 : // = (partial b / partial p) - (partial A / partial p) * u
305 :
306 : // We first solve for (partial u / partial p) for each parameter:
307 : // -A * (partial u / partial p) = - (partial R / partial p)
308 :
309 : this->sensitivity_solve(parameters);
310 :
311 : // Get ready to fill in sensitivities:
312 : sensitivities.allocate_data(qoi_indices, *this, parameters);
313 :
314 : // We use the identity:
315 : // dq/dp = (partial q / partial p) + (partial q / partial u) *
316 : // (partial u / partial p)
317 :
318 : // We get (partial q / partial u) from the user
319 : this->assemble_qoi_derivative(qoi_indices);
320 :
321 : for (unsigned int j=0; j != Np; ++j)
322 : {
323 : // We currently get partial derivatives via central differencing
324 : Number delta_p = 1e-6;
325 :
326 : // (partial q / partial p) ~= (q(p+dp)-q(p-dp))/(2*dp)
327 :
328 : Number old_parameter = *parameters[j];
329 :
330 : *parameters[j] = old_parameter - delta_p;
331 : this->assemble_qoi(qoi_indices);
332 : std::vector<Number> qoi_minus = this->qoi;
333 :
334 : *parameters[j] = old_parameter + delta_p;
335 : this->assemble_qoi(qoi_indices);
336 : std::vector<Number> & qoi_plus = this->qoi;
337 : std::vector<Number> partialq_partialp(Nq, 0);
338 : for (unsigned int i=0; i != Nq; ++i)
339 : if (qoi_indices.has_index(i))
340 : partialq_partialp[i] = (qoi_plus[i] - qoi_minus[i]) / (2.*delta_p);
341 :
342 : for (unsigned int i=0; i != Nq; ++i)
343 : if (qoi_indices.has_index(i))
344 : sensitivities[i][j] = partialq_partialp[i] +
345 : this->get_adjoint_rhs(i).dot(this->get_sensitivity_solution(i));
346 : }
347 :
348 : // All parameters have been reset.
349 : // Don't leave the qoi or system changed - principle of least
350 : // surprise.
351 : this->assemble();
352 : this->rhs->close();
353 : this->matrix->close();
354 : this->assemble_qoi(qoi_indices);
355 : }
356 : */
357 :
358 :
359 :
360 289244 : LinearSolver<Number> * LinearImplicitSystem::get_linear_solver() const
361 : {
362 289244 : return linear_solver.get();
363 : }
364 :
365 :
366 :
367 0 : void LinearImplicitSystem::assembly(bool,
368 : bool,
369 : bool,
370 : bool)
371 : {
372 : // Residual R(u(p),p) := A(p)*u(p) - b(p)
373 : // partial R / partial u = A
374 :
375 0 : this->assemble();
376 0 : this->rhs->close();
377 0 : this->matrix->close();
378 :
379 0 : *(this->rhs) *= -1.0;
380 0 : this->rhs->add_vector(*this->solution, *this->matrix);
381 0 : }
382 :
383 : } // namespace libMesh
|
