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