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 : // Local includes
21 : #include "libmesh/dof_map.h"
22 : #include "libmesh/equation_systems.h"
23 : #include "libmesh/implicit_system.h"
24 : #include "libmesh/int_range.h"
25 : #include "libmesh/libmesh_logging.h"
26 : #include "libmesh/linear_solver.h"
27 : #include "libmesh/mesh_base.h"
28 : #include "libmesh/numeric_vector.h"
29 : #include "libmesh/parameters.h"
30 : #include "libmesh/parameter_vector.h"
31 : #include "libmesh/qoi_set.h"
32 : #include "libmesh/sensitivity_data.h"
33 : #include "libmesh/sparse_matrix.h"
34 : #include "libmesh/diagonal_matrix.h"
35 : #include "libmesh/utility.h"
36 : #include "libmesh/static_condensation.h"
37 : #include "libmesh/static_condensation_preconditioner.h"
38 : #include "libmesh/nonlinear_solver.h"
39 :
40 : namespace libMesh
41 : {
42 :
43 18022 : ImplicitSystem::ImplicitSystem (EquationSystems & es,
44 : const std::string & name_in,
45 18022 : const unsigned int number_in) :
46 :
47 : Parent (es, name_in, number_in),
48 16974 : matrix (nullptr),
49 18022 : zero_out_matrix_and_rhs(true)
50 : {
51 18022 : if (this->has_static_condensation())
52 280 : this->create_static_condensation_system_matrix();
53 18022 : }
54 :
55 16974 : ImplicitSystem::~ImplicitSystem () = default;
56 :
57 350 : void ImplicitSystem::create_static_condensation_system_matrix()
58 : {
59 350 : auto sc_system_matrix = std::make_unique<StaticCondensation>(this->get_mesh(), *this, this->get_dof_map(), this->get_dof_map().get_static_condensation());
60 350 : _sc_system_matrix = sc_system_matrix.get();
61 350 : matrix = &(this->add_matrix ("System Matrix", std::move(sc_system_matrix)));
62 350 : }
63 :
64 70 : void ImplicitSystem::create_static_condensation ()
65 : {
66 70 : Parent::create_static_condensation();
67 70 : this->create_static_condensation_system_matrix();
68 70 : }
69 :
70 : template <typename T>
71 350 : void ImplicitSystem::setup_static_condensation_preconditioner(T & solver)
72 : {
73 10 : libmesh_assert(_sc_system_matrix);
74 360 : solver.attach_preconditioner(&_sc_system_matrix->get_preconditioner());
75 350 : }
76 :
77 :
78 427 : void ImplicitSystem::clear ()
79 : {
80 : // clear the parent data
81 427 : Parent::clear();
82 :
83 : // Restore us to a "basic" state
84 427 : matrix = nullptr;
85 427 : _sc_system_matrix = nullptr;
86 427 : }
87 :
88 :
89 :
90 38939 : void ImplicitSystem::assemble ()
91 : {
92 1344 : libmesh_assert(matrix);
93 1344 : libmesh_assert (matrix->initialized());
94 1344 : libmesh_assert(rhs);
95 1344 : libmesh_assert (rhs->initialized());
96 :
97 38939 : if (zero_out_matrix_and_rhs)
98 : {
99 38939 : matrix->zero ();
100 38939 : rhs->zero ();
101 : }
102 :
103 : // Call the base class assemble function
104 38939 : Parent::assemble ();
105 38939 : }
106 :
107 :
108 :
109 18022 : void ImplicitSystem::add_matrices ()
110 : {
111 18022 : Parent::add_matrices();
112 :
113 : // Possible that we cleared the _matrices but
114 : // forgot to update the matrix pointer?
115 18022 : if (this->n_matrices() == 0)
116 17389 : matrix = nullptr;
117 :
118 : // Only need to add the matrix if it isn't there
119 : // already!
120 18022 : if (matrix == nullptr)
121 17672 : matrix = &(this->add_matrix ("System Matrix"));
122 :
123 524 : libmesh_assert(matrix);
124 18022 : }
125 :
126 :
127 :
128 141120 : void ImplicitSystem::disable_cache () {
129 141120 : this->assemble_before_solve = true;
130 141120 : this->get_linear_solver()->reuse_preconditioner(false);
131 141120 : }
132 :
133 :
134 :
135 : std::pair<unsigned int, Real>
136 71 : ImplicitSystem::sensitivity_solve (const ParameterVector & parameters_vec)
137 : {
138 : // Log how long the linear solve takes.
139 2 : LOG_SCOPE("sensitivity_solve()", "ImplicitSystem");
140 :
141 : // The forward system should now already be solved.
142 : // Now assemble the corresponding sensitivity system.
143 :
144 71 : if (this->assemble_before_solve)
145 : {
146 : // Build the Jacobian
147 71 : this->assembly(false, true);
148 71 : this->matrix->close();
149 :
150 : // Reset and build the RHS from the residual derivatives
151 71 : this->assemble_residual_derivatives(parameters_vec);
152 : }
153 :
154 : // The sensitivity problem is linear
155 71 : LinearSolver<Number> * solver = this->get_linear_solver();
156 :
157 : // Our iteration counts and residuals will be sums of the individual
158 : // results
159 : std::pair<unsigned int, Real> solver_params =
160 71 : this->get_linear_solve_parameters();
161 2 : std::pair<unsigned int, Real> totalrval = std::make_pair(0,0.0);
162 :
163 : // Solve the linear system.
164 71 : SparseMatrix<Number> * pc = this->request_matrix("Preconditioner");
165 213 : for (auto p : make_range(parameters_vec.size()))
166 : {
167 : std::pair<unsigned int, Real> rval =
168 142 : solver->solve (*matrix, pc,
169 : this->add_sensitivity_solution(p),
170 : this->get_sensitivity_rhs(p),
171 138 : double(solver_params.second),
172 8 : solver_params.first);
173 :
174 142 : totalrval.first += rval.first;
175 142 : totalrval.second += rval.second;
176 : }
177 :
178 : // The linear solver may not have fit our constraints exactly
179 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
180 213 : for (auto p : make_range(parameters_vec.size()))
181 4 : this->get_dof_map().enforce_constraints_exactly
182 142 : (*this, &this->get_sensitivity_solution(p),
183 : /* homogeneous = */ true);
184 : #endif
185 :
186 73 : return totalrval;
187 : }
188 :
189 :
190 :
191 : std::pair<unsigned int, Real>
192 14583 : ImplicitSystem::adjoint_solve (const QoISet & qoi_indices)
193 : {
194 : // Log how long the linear solve takes.
195 462 : LOG_SCOPE("adjoint_solve()", "ImplicitSystem");
196 :
197 14583 : if (this->assemble_before_solve)
198 : // Assemble the linear system
199 14583 : this->assembly (/* get_residual = */ false,
200 924 : /* get_jacobian = */ true);
201 :
202 : // The adjoint problem is linear
203 14583 : LinearSolver<Number> * solver = this->get_linear_solver();
204 :
205 : // Reset and build the RHS from the QOI derivative
206 15507 : this->assemble_qoi_derivative(qoi_indices,
207 : /* include_liftfunc = */ false,
208 924 : /* apply_constraints = */ true);
209 :
210 : // Our iteration counts and residuals will be sums of the individual
211 : // results
212 : std::pair<unsigned int, Real> solver_params =
213 14583 : this->get_linear_solve_parameters();
214 462 : std::pair<unsigned int, Real> totalrval = std::make_pair(0,0.0);
215 :
216 39329 : for (auto i : make_range(this->n_qois()))
217 24746 : if (qoi_indices.has_index(i))
218 : {
219 : const std::pair<unsigned int, Real> rval =
220 24746 : solver->adjoint_solve (*matrix, this->add_adjoint_solution(i),
221 : this->get_adjoint_rhs(i),
222 23948 : double(solver_params.second),
223 2394 : solver_params.first);
224 :
225 24746 : totalrval.first += rval.first;
226 24746 : totalrval.second += rval.second;
227 : }
228 :
229 : // The linear solver may not have fit our constraints exactly
230 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
231 39329 : for (auto i : make_range(this->n_qois()))
232 24746 : if (qoi_indices.has_index(i))
233 798 : this->get_dof_map().enforce_adjoint_constraints_exactly
234 24746 : (this->get_adjoint_solution(i), i);
235 : #endif
236 :
237 15045 : return totalrval;
238 : }
239 :
240 :
241 :
242 : std::pair<unsigned int, Real>
243 0 : ImplicitSystem::weighted_sensitivity_adjoint_solve (const ParameterVector & parameters_in,
244 : const ParameterVector & weights,
245 : const QoISet & qoi_indices)
246 : {
247 : // Log how long the linear solve takes.
248 0 : LOG_SCOPE("weighted_sensitivity_adjoint_solve()", "ImplicitSystem");
249 :
250 : // We currently get partial derivatives via central differencing
251 0 : const Real delta_p = TOLERANCE;
252 :
253 0 : ParameterVector & parameters_vec =
254 : const_cast<ParameterVector &>(parameters_in);
255 :
256 : // The forward system should now already be solved.
257 : // The adjoint system should now already be solved.
258 : // Now we're assembling a weighted sum of adjoint-adjoint systems:
259 : //
260 : // dR/du (u, sum_l(w_l*z^l)) = sum_l(w_l*(Q''_ul - R''_ul (u, z)))
261 :
262 : // FIXME: The derivation here does not yet take adjoint boundary
263 : // conditions into account.
264 : #ifdef LIBMESH_ENABLE_DIRICHLET
265 0 : for (auto i : make_range(this->n_qois()))
266 0 : if (qoi_indices.has_index(i))
267 0 : libmesh_assert(!this->get_dof_map().has_adjoint_dirichlet_boundaries(i));
268 : #endif
269 :
270 : // We'll assemble the rhs first, because the R'' term will require
271 : // perturbing the jacobian
272 :
273 : // We'll use temporary rhs vectors, because we haven't (yet) found
274 : // any good reasons why users might want to save these:
275 :
276 0 : std::vector<std::unique_ptr<NumericVector<Number>>> temprhs(this->n_qois());
277 0 : for (auto i : make_range(this->n_qois()))
278 0 : if (qoi_indices.has_index(i))
279 0 : temprhs[i] = this->rhs->zero_clone();
280 :
281 : // We approximate the _l partial derivatives via a central
282 : // differencing perturbation in the w_l direction:
283 : //
284 : // sum_l(w_l*v_l) ~= (v(p + dp*w_l*e_l) - v(p - dp*w_l*e_l))/(2*dp)
285 :
286 : // PETSc doesn't implement SGEMX, so neither does NumericVector,
287 : // so we want to avoid calculating f -= R'*z. We'll thus evaluate
288 : // the above equation by first adding -v(p+dp...), then multiplying
289 : // the intermediate result vectors by -1, then adding -v(p-dp...),
290 : // then finally dividing by 2*dp.
291 :
292 0 : ParameterVector oldparameters, parameterperturbation;
293 0 : parameters_vec.deep_copy(oldparameters);
294 0 : weights.deep_copy(parameterperturbation);
295 0 : parameterperturbation *= delta_p;
296 0 : parameters_vec += parameterperturbation;
297 :
298 0 : this->assembly(false, true);
299 0 : this->matrix->close();
300 :
301 : // Take the discrete adjoint, so that we can calculate R_u(u,z) with
302 : // a matrix-vector product of R_u and z.
303 0 : matrix->get_transpose(*matrix);
304 :
305 0 : this->assemble_qoi_derivative(qoi_indices,
306 : /* include_liftfunc = */ false,
307 0 : /* apply_constraints = */ true);
308 0 : for (auto i : make_range(this->n_qois()))
309 0 : if (qoi_indices.has_index(i))
310 : {
311 0 : this->get_adjoint_rhs(i).close();
312 0 : *(temprhs[i]) -= this->get_adjoint_rhs(i);
313 0 : this->matrix->vector_mult_add(*(temprhs[i]), this->get_adjoint_solution(i));
314 0 : *(temprhs[i]) *= -1.0;
315 : }
316 :
317 0 : oldparameters.value_copy(parameters_vec);
318 0 : parameterperturbation *= -1.0;
319 0 : parameters_vec += parameterperturbation;
320 :
321 0 : this->assembly(false, true);
322 0 : this->matrix->close();
323 0 : matrix->get_transpose(*matrix);
324 :
325 0 : this->assemble_qoi_derivative(qoi_indices,
326 : /* include_liftfunc = */ false,
327 0 : /* apply_constraints = */ true);
328 0 : for (auto i : make_range(this->n_qois()))
329 0 : if (qoi_indices.has_index(i))
330 : {
331 0 : this->get_adjoint_rhs(i).close();
332 0 : *(temprhs[i]) -= this->get_adjoint_rhs(i);
333 0 : this->matrix->vector_mult_add(*(temprhs[i]), this->get_adjoint_solution(i));
334 0 : *(temprhs[i]) /= (2.0*delta_p);
335 : }
336 :
337 : // Finally, assemble the jacobian at the non-perturbed parameter
338 : // values. Ignore assemble_before_solve; if we had a good
339 : // non-perturbed matrix before we've already overwritten it.
340 0 : oldparameters.value_copy(parameters_vec);
341 :
342 : // if (this->assemble_before_solve)
343 : {
344 : // Build the Jacobian
345 0 : this->assembly(false, true);
346 0 : this->matrix->close();
347 :
348 : // Take the discrete adjoint
349 0 : matrix->get_transpose(*matrix);
350 : }
351 :
352 : // The weighted adjoint-adjoint problem is linear
353 0 : LinearSolver<Number> * solver = this->get_linear_solver();
354 :
355 : // Our iteration counts and residuals will be sums of the individual
356 : // results
357 : std::pair<unsigned int, Real> solver_params =
358 0 : this->get_linear_solve_parameters();
359 0 : std::pair<unsigned int, Real> totalrval = std::make_pair(0,0.0);
360 :
361 0 : for (auto i : make_range(this->n_qois()))
362 0 : if (qoi_indices.has_index(i))
363 : {
364 : const std::pair<unsigned int, Real> rval =
365 0 : solver->solve (*matrix, this->add_weighted_sensitivity_adjoint_solution(i),
366 0 : *(temprhs[i]),
367 0 : double(solver_params.second),
368 0 : solver_params.first);
369 :
370 0 : totalrval.first += rval.first;
371 0 : totalrval.second += rval.second;
372 : }
373 :
374 : // The linear solver may not have fit our constraints exactly
375 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
376 0 : for (auto i : make_range(this->n_qois()))
377 0 : if (qoi_indices.has_index(i))
378 0 : this->get_dof_map().enforce_constraints_exactly
379 0 : (*this, &this->get_weighted_sensitivity_adjoint_solution(i),
380 : /* homogeneous = */ true);
381 : #endif
382 :
383 0 : return totalrval;
384 0 : }
385 :
386 :
387 :
388 : std::pair<unsigned int, Real>
389 0 : ImplicitSystem::weighted_sensitivity_solve (const ParameterVector & parameters_in,
390 : const ParameterVector & weights)
391 : {
392 : // Log how long the linear solve takes.
393 0 : LOG_SCOPE("weighted_sensitivity_solve()", "ImplicitSystem");
394 :
395 : // We currently get partial derivatives via central differencing
396 0 : const Real delta_p = TOLERANCE;
397 :
398 0 : ParameterVector & parameters_vec =
399 : const_cast<ParameterVector &>(parameters_in);
400 :
401 : // The forward system should now already be solved.
402 :
403 : // Now we're assembling a weighted sum of sensitivity systems:
404 : //
405 : // dR/du (u, v)(sum(w_l*u'_l)) = -sum_l(w_l*R'_l (u, v)) forall v
406 :
407 : // We'll assemble the rhs first, because the R' term will require
408 : // perturbing the system, and some applications may not be able to
409 : // assemble a perturbed residual without simultaneously constructing
410 : // a perturbed jacobian.
411 :
412 : // We approximate the _l partial derivatives via a central
413 : // differencing perturbation in the w_l direction:
414 : //
415 : // sum_l(w_l*v_l) ~= (v(p + dp*w_l*e_l) - v(p - dp*w_l*e_l))/(2*dp)
416 :
417 0 : ParameterVector oldparameters, parameterperturbation;
418 0 : parameters_vec.deep_copy(oldparameters);
419 0 : weights.deep_copy(parameterperturbation);
420 0 : parameterperturbation *= delta_p;
421 0 : parameters_vec += parameterperturbation;
422 :
423 0 : this->assembly(true, false, true);
424 0 : this->rhs->close();
425 :
426 0 : std::unique_ptr<NumericVector<Number>> temprhs = this->rhs->clone();
427 :
428 0 : oldparameters.value_copy(parameters_vec);
429 0 : parameterperturbation *= -1.0;
430 0 : parameters_vec += parameterperturbation;
431 :
432 0 : this->assembly(true, false, true);
433 0 : this->rhs->close();
434 :
435 0 : *temprhs -= *(this->rhs);
436 0 : *temprhs /= (2.0*delta_p);
437 :
438 : // Finally, assemble the jacobian at the non-perturbed parameter
439 : // values
440 0 : oldparameters.value_copy(parameters_vec);
441 :
442 : // Build the Jacobian
443 0 : this->assembly(false, true);
444 0 : this->matrix->close();
445 :
446 : // The weighted sensitivity problem is linear
447 0 : LinearSolver<Number> * solver = this->get_linear_solver();
448 :
449 : std::pair<unsigned int, Real> solver_params =
450 0 : this->get_linear_solve_parameters();
451 :
452 : const std::pair<unsigned int, Real> rval =
453 0 : solver->solve (*matrix, this->add_weighted_sensitivity_solution(),
454 0 : *temprhs,
455 0 : double(solver_params.second),
456 0 : solver_params.first);
457 :
458 : // The linear solver may not have fit our constraints exactly
459 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
460 0 : this->get_dof_map().enforce_constraints_exactly
461 0 : (*this, &this->get_weighted_sensitivity_solution(),
462 : /* homogeneous = */ true);
463 : #endif
464 :
465 0 : return rval;
466 0 : }
467 :
468 :
469 :
470 7075 : void ImplicitSystem::assemble_residual_derivatives(const ParameterVector & parameters_in)
471 : {
472 202 : ParameterVector & parameters_vec =
473 : const_cast<ParameterVector &>(parameters_in);
474 :
475 : const unsigned int Np = cast_int<unsigned int>
476 202 : (parameters_vec.size());
477 :
478 14505 : for (unsigned int p=0; p != Np; ++p)
479 : {
480 7430 : NumericVector<Number> & sensitivity_rhs = this->add_sensitivity_rhs(p);
481 :
482 : // Approximate -(partial R / partial p) by
483 : // (R(p-dp) - R(p+dp)) / (2*dp)
484 :
485 7430 : Number old_parameter = *parameters_vec[p];
486 :
487 : const Real delta_p =
488 7642 : TOLERANCE * std::max(std::abs(old_parameter), 1e-3);
489 :
490 7218 : *parameters_vec[p] -= delta_p;
491 :
492 : // this->assembly(true, false, true);
493 7430 : this->assembly(true, false, false);
494 7430 : this->rhs->close();
495 7430 : sensitivity_rhs = *this->rhs;
496 :
497 7430 : *parameters_vec[p] = old_parameter + delta_p;
498 :
499 : // this->assembly(true, false, true);
500 7430 : this->assembly(true, false, false);
501 7430 : this->rhs->close();
502 :
503 7430 : sensitivity_rhs -= *this->rhs;
504 7430 : sensitivity_rhs /= (2*delta_p);
505 7430 : sensitivity_rhs.close();
506 :
507 212 : *parameters_vec[p] = old_parameter;
508 : }
509 7075 : }
510 :
511 :
512 :
513 7004 : void ImplicitSystem::adjoint_qoi_parameter_sensitivity (const QoISet & qoi_indices,
514 : const ParameterVector & parameters_in,
515 : SensitivityData & sensitivities)
516 : {
517 200 : ParameterVector & parameters_vec =
518 : const_cast<ParameterVector &>(parameters_in);
519 :
520 : const unsigned int Np = cast_int<unsigned int>
521 200 : (parameters_vec.size());
522 200 : const unsigned int Nq = this->n_qois();
523 :
524 : // An introduction to the problem:
525 : //
526 : // Residual R(u(p),p) = 0
527 : // partial R / partial u = J = system matrix
528 : //
529 : // This implies that:
530 : // d/dp(R) = 0
531 : // (partial R / partial p) +
532 : // (partial R / partial u) * (partial u / partial p) = 0
533 :
534 : // We first do an adjoint solve:
535 : // J^T * z = (partial q / partial u)
536 : // if we haven't already or dont have an initial condition for the adjoint
537 7004 : if (!this->is_adjoint_already_solved())
538 : {
539 284 : this->adjoint_solve(qoi_indices);
540 : }
541 :
542 7004 : this->assemble_residual_derivatives(parameters_in);
543 :
544 : // Get ready to fill in sensitivities:
545 7004 : sensitivities.allocate_data(qoi_indices, *this, parameters_vec);
546 :
547 : // We use the identities:
548 : // dq/dp = (partial q / partial p) + (partial q / partial u) *
549 : // (partial u / partial p)
550 : // dq/dp = (partial q / partial p) + (J^T * z) *
551 : // (partial u / partial p)
552 : // dq/dp = (partial q / partial p) + z * J *
553 : // (partial u / partial p)
554 :
555 : // Leading to our final formula:
556 : // dq/dp = (partial q / partial p) - z * (partial R / partial p)
557 :
558 : // In the case of adjoints with heterogenous Dirichlet boundary
559 : // function phi, where
560 : // q := S(u) - R(u,phi)
561 : // the final formula works out to:
562 : // dq/dp = (partial S / partial p) - z * (partial R / partial p)
563 : // Because we currently have no direct access to
564 : // (partial S / partial p), we use the identity
565 : // (partial S / partial p) = (partial q / partial p) +
566 : // phi * (partial R / partial p)
567 : // to derive an equivalent equation:
568 : // dq/dp = (partial q / partial p) - (z-phi) * (partial R / partial p)
569 :
570 : // Since z-phi degrees of freedom are zero for constrained indices,
571 : // we can use the same constrained -(partial R / partial p) that we
572 : // use for forward sensitivity solves, taking into account the
573 : // differing sign convention.
574 : //
575 : // Since that vector is constrained, its constrained indices are
576 : // zero, so its product with phi is zero, so we can neglect the
577 : // evaluation of phi terms.
578 :
579 14292 : for (unsigned int j=0; j != Np; ++j)
580 : {
581 : // We currently get partial derivatives via central differencing
582 :
583 : // (partial q / partial p) ~= (q(p+dp)-q(p-dp))/(2*dp)
584 : // (partial R / partial p) ~= (rhs(p+dp) - rhs(p-dp))/(2*dp)
585 :
586 7288 : Number old_parameter = *parameters_vec[j];
587 :
588 : const Real delta_p =
589 7496 : TOLERANCE * std::max(std::abs(old_parameter), 1e-3);
590 :
591 7288 : *parameters_vec[j] = old_parameter - delta_p;
592 7288 : this->assemble_qoi(qoi_indices);
593 7496 : const std::vector<Number> qoi_minus = this->get_qoi_values();
594 :
595 7288 : NumericVector<Number> & neg_partialR_partialp = this->get_sensitivity_rhs(j);
596 :
597 7288 : *parameters_vec[j] = old_parameter + delta_p;
598 7288 : this->assemble_qoi(qoi_indices);
599 7496 : const std::vector<Number> qoi_plus = this->get_qoi_values();
600 :
601 7496 : std::vector<Number> partialq_partialp(Nq, 0);
602 14576 : for (unsigned int i=0; i != Nq; ++i)
603 7288 : if (qoi_indices.has_index(i))
604 7912 : partialq_partialp[i] = (qoi_plus[i] - qoi_minus[i]) / (2.*delta_p);
605 :
606 : // Don't leave the parameter changed
607 208 : *parameters_vec[j] = old_parameter;
608 :
609 14576 : for (unsigned int i=0; i != Nq; ++i)
610 7288 : if (qoi_indices.has_index(i))
611 7912 : sensitivities[i][j] = partialq_partialp[i] +
612 7288 : neg_partialR_partialp.dot(this->get_adjoint_solution(i));
613 : }
614 :
615 : // All parameters_vec have been reset.
616 : // Reset the original qoi.
617 :
618 7004 : this->assemble_qoi(qoi_indices);
619 7004 : }
620 :
621 :
622 :
623 71 : void ImplicitSystem::forward_qoi_parameter_sensitivity (const QoISet & qoi_indices,
624 : const ParameterVector & parameters_in,
625 : SensitivityData & sensitivities)
626 : {
627 2 : ParameterVector & parameters_vec =
628 : const_cast<ParameterVector &>(parameters_in);
629 :
630 : const unsigned int Np = cast_int<unsigned int>
631 2 : (parameters_vec.size());
632 2 : const unsigned int Nq = this->n_qois();
633 :
634 : // An introduction to the problem:
635 : //
636 : // Residual R(u(p),p) = 0
637 : // partial R / partial u = J = system matrix
638 : //
639 : // This implies that:
640 : // d/dp(R) = 0
641 : // (partial R / partial p) +
642 : // (partial R / partial u) * (partial u / partial p) = 0
643 :
644 : // We first solve for (partial u / partial p) for each parameter:
645 : // J * (partial u / partial p) = - (partial R / partial p)
646 :
647 71 : this->sensitivity_solve(parameters_vec);
648 :
649 : // Get ready to fill in sensitivities:
650 71 : sensitivities.allocate_data(qoi_indices, *this, parameters_vec);
651 :
652 : // We use the identity:
653 : // dq/dp = (partial q / partial p) + (partial q / partial u) *
654 : // (partial u / partial p)
655 :
656 : // We get (partial q / partial u) from the user
657 75 : this->assemble_qoi_derivative(qoi_indices,
658 : /* include_liftfunc = */ true,
659 4 : /* apply_constraints = */ false);
660 :
661 : // We don't need these to be closed() in this function, but libMesh
662 : // standard practice is to have them closed() by the time the
663 : // function exits
664 142 : for (auto i : make_range(this->n_qois()))
665 71 : if (qoi_indices.has_index(i))
666 71 : this->get_adjoint_rhs(i).close();
667 :
668 213 : for (unsigned int j=0; j != Np; ++j)
669 : {
670 : // We currently get partial derivatives via central differencing
671 :
672 : // (partial q / partial p) ~= (q(p+dp)-q(p-dp))/(2*dp)
673 :
674 142 : Number old_parameter = *parameters_vec[j];
675 :
676 : const Real delta_p =
677 146 : TOLERANCE * std::max(std::abs(old_parameter), 1e-3);
678 :
679 142 : *parameters_vec[j] = old_parameter - delta_p;
680 142 : this->assemble_qoi(qoi_indices);
681 146 : const std::vector<Number> qoi_minus = this->get_qoi_values();
682 :
683 142 : *parameters_vec[j] = old_parameter + delta_p;
684 142 : this->assemble_qoi(qoi_indices);
685 146 : const std::vector<Number> qoi_plus = this->get_qoi_values();
686 :
687 146 : std::vector<Number> partialq_partialp(Nq, 0);
688 284 : for (unsigned int i=0; i != Nq; ++i)
689 142 : if (qoi_indices.has_index(i))
690 154 : partialq_partialp[i] = (qoi_plus[i] - qoi_minus[i]) / (2.*delta_p);
691 :
692 : // Don't leave the parameter changed
693 4 : *parameters_vec[j] = old_parameter;
694 :
695 284 : for (unsigned int i=0; i != Nq; ++i)
696 142 : if (qoi_indices.has_index(i))
697 154 : sensitivities[i][j] = partialq_partialp[i] +
698 142 : this->get_adjoint_rhs(i).dot(this->get_sensitivity_solution(j));
699 : }
700 :
701 : // All parameters_vec have been reset.
702 : // We didn't cache the original rhs or matrix for memory reasons,
703 : // but we can restore them to a state consistent solution -
704 : // principle of least surprise.
705 71 : this->assembly(true, true);
706 71 : this->rhs->close();
707 71 : this->matrix->close();
708 71 : this->assemble_qoi(qoi_indices);
709 71 : }
710 :
711 :
712 :
713 0 : void ImplicitSystem::qoi_parameter_hessian_vector_product (const QoISet & qoi_indices,
714 : const ParameterVector & parameters_in,
715 : const ParameterVector & vector,
716 : SensitivityData & sensitivities)
717 : {
718 : // We currently get partial derivatives via finite differencing
719 0 : const Real delta_p = TOLERANCE;
720 :
721 0 : ParameterVector & parameters_vec =
722 : const_cast<ParameterVector &>(parameters_in);
723 :
724 : // We'll use a single temporary vector for matrix-vector-vector products
725 0 : std::unique_ptr<NumericVector<Number>> tempvec = this->solution->zero_clone();
726 :
727 : const unsigned int Np = cast_int<unsigned int>
728 0 : (parameters_vec.size());
729 0 : const unsigned int Nq = this->n_qois();
730 :
731 : // For each quantity of interest q, the parameter sensitivity
732 : // Hessian is defined as q''_{kl} = {d^2 q}/{d p_k d p_l}.
733 : // Given a vector of parameter perturbation weights w_l, this
734 : // function evaluates the hessian-vector product sum_l(q''_{kl}*w_l)
735 : //
736 : // We calculate it from values and partial derivatives of the
737 : // quantity of interest function Q, solution u, adjoint solution z,
738 : // parameter sensitivity adjoint solutions z^l, and residual R, as:
739 : //
740 : // sum_l(q''_{kl}*w_l) =
741 : // sum_l(w_l * Q''_{kl}) + Q''_{uk}(u)*(sum_l(w_l u'_l)) -
742 : // R'_k(u, sum_l(w_l*z^l)) - R'_{uk}(u,z)*(sum_l(w_l u'_l) -
743 : // sum_l(w_l*R''_{kl}(u,z))
744 : //
745 : // See the adjoints model document for more details.
746 :
747 : // We first do an adjoint solve to get z for each quantity of
748 : // interest
749 : // if we haven't already or dont have an initial condition for the adjoint
750 0 : if (!this->is_adjoint_already_solved())
751 : {
752 0 : this->adjoint_solve(qoi_indices);
753 : }
754 :
755 : // Get ready to fill in sensitivities:
756 0 : sensitivities.allocate_data(qoi_indices, *this, parameters_vec);
757 :
758 : // We can't solve for all the solution sensitivities u'_l or for all
759 : // of the parameter sensitivity adjoint solutions z^l without
760 : // requiring O(Nq*Np) linear solves. So we'll solve directly for their
761 : // weighted sum - this is just O(Nq) solves.
762 :
763 : // First solve for sum_l(w_l u'_l).
764 0 : this->weighted_sensitivity_solve(parameters_vec, vector);
765 :
766 : // Then solve for sum_l(w_l z^l).
767 0 : this->weighted_sensitivity_adjoint_solve(parameters_vec, vector, qoi_indices);
768 :
769 0 : for (unsigned int k=0; k != Np; ++k)
770 : {
771 : // We approximate sum_l(w_l * Q''_{kl}) with a central
772 : // differencing perturbation:
773 : // sum_l(w_l * Q''_{kl}) ~=
774 : // (Q(p + dp*w_l*e_l + dp*e_k) - Q(p - dp*w_l*e_l + dp*e_k) -
775 : // Q(p + dp*w_l*e_l - dp*e_k) + Q(p - dp*w_l*e_l - dp*e_k))/(4*dp^2)
776 :
777 : // The sum(w_l*R''_kl) term requires the same sort of perturbation,
778 : // and so we subtract it in at the same time:
779 : // sum_l(w_l * R''_{kl}) ~=
780 : // (R(p + dp*w_l*e_l + dp*e_k) - R(p - dp*w_l*e_l + dp*e_k) -
781 : // R(p + dp*w_l*e_l - dp*e_k) + R(p - dp*w_l*e_l - dp*e_k))/(4*dp^2)
782 :
783 0 : ParameterVector oldparameters, parameterperturbation;
784 0 : parameters_vec.deep_copy(oldparameters);
785 0 : vector.deep_copy(parameterperturbation);
786 0 : parameterperturbation *= delta_p;
787 0 : parameters_vec += parameterperturbation;
788 :
789 0 : Number old_parameter = *parameters_vec[k];
790 :
791 0 : *parameters_vec[k] = old_parameter + delta_p;
792 0 : this->assemble_qoi(qoi_indices);
793 0 : this->assembly(true, false, true);
794 0 : this->rhs->close();
795 0 : std::vector<Number> partial2q_term = this->get_qoi_values();
796 0 : std::vector<Number> partial2R_term(this->n_qois());
797 0 : for (unsigned int i=0; i != Nq; ++i)
798 0 : if (qoi_indices.has_index(i))
799 0 : partial2R_term[i] = this->rhs->dot(this->get_adjoint_solution(i));
800 :
801 0 : *parameters_vec[k] = old_parameter - delta_p;
802 0 : this->assemble_qoi(qoi_indices);
803 0 : this->assembly(true, false, true);
804 0 : this->rhs->close();
805 0 : for (unsigned int i=0; i != Nq; ++i)
806 0 : if (qoi_indices.has_index(i))
807 : {
808 0 : partial2q_term[i] -= this->get_qoi_value(i);
809 0 : partial2R_term[i] -= this->rhs->dot(this->get_adjoint_solution(i));
810 : }
811 :
812 0 : oldparameters.value_copy(parameters_vec);
813 0 : parameterperturbation *= -1.0;
814 0 : parameters_vec += parameterperturbation;
815 :
816 : // Re-center old_parameter, which may be affected by vector
817 0 : old_parameter = *parameters_vec[k];
818 :
819 0 : *parameters_vec[k] = old_parameter + delta_p;
820 0 : this->assemble_qoi(qoi_indices);
821 0 : this->assembly(true, false, true);
822 0 : this->rhs->close();
823 0 : for (unsigned int i=0; i != Nq; ++i)
824 0 : if (qoi_indices.has_index(i))
825 : {
826 0 : partial2q_term[i] -= this->get_qoi_value(i);
827 0 : partial2R_term[i] -= this->rhs->dot(this->get_adjoint_solution(i));
828 : }
829 :
830 0 : *parameters_vec[k] = old_parameter - delta_p;
831 0 : this->assemble_qoi(qoi_indices);
832 0 : this->assembly(true, false, true);
833 0 : this->rhs->close();
834 0 : for (unsigned int i=0; i != Nq; ++i)
835 0 : if (qoi_indices.has_index(i))
836 : {
837 0 : partial2q_term[i] += this->get_qoi_value(i);
838 0 : partial2R_term[i] += this->rhs->dot(this->get_adjoint_solution(i));
839 : }
840 :
841 0 : for (unsigned int i=0; i != Nq; ++i)
842 0 : if (qoi_indices.has_index(i))
843 : {
844 0 : partial2q_term[i] /= (4. * delta_p * delta_p);
845 0 : partial2R_term[i] /= (4. * delta_p * delta_p);
846 : }
847 :
848 0 : for (unsigned int i=0; i != Nq; ++i)
849 0 : if (qoi_indices.has_index(i))
850 0 : sensitivities[i][k] = partial2q_term[i] - partial2R_term[i];
851 :
852 : // We get (partial q / partial u), R, and
853 : // (partial R / partial u) from the user, but centrally
854 : // difference to get q_uk, R_k, and R_uk terms:
855 : // (partial R / partial k)
856 : // R_k*sum(w_l*z^l) = (R(p+dp*e_k)*sum(w_l*z^l) - R(p-dp*e_k)*sum(w_l*z^l))/(2*dp)
857 : // (partial^2 q / partial u partial k)
858 : // q_uk = (q_u(p+dp*e_k) - q_u(p-dp*e_k))/(2*dp)
859 : // (partial^2 R / partial u partial k)
860 : // R_uk*z*sum(w_l*u'_l) = (R_u(p+dp*e_k)*z*sum(w_l*u'_l) - R_u(p-dp*e_k)*z*sum(w_l*u'_l))/(2*dp)
861 :
862 : // To avoid creating Nq temporary vectors for q_uk or R_uk, we add
863 : // subterms to the sensitivities output one by one.
864 : //
865 : // FIXME: this is probably a bad order of operations for
866 : // controlling floating point error.
867 :
868 0 : *parameters_vec[k] = old_parameter + delta_p;
869 0 : this->assembly(true, true);
870 0 : this->rhs->close();
871 0 : this->matrix->close();
872 0 : this->assemble_qoi_derivative(qoi_indices,
873 : /* include_liftfunc = */ true,
874 0 : /* apply_constraints = */ false);
875 :
876 0 : this->matrix->vector_mult(*tempvec, this->get_weighted_sensitivity_solution());
877 :
878 0 : for (unsigned int i=0; i != Nq; ++i)
879 0 : if (qoi_indices.has_index(i))
880 : {
881 0 : this->get_adjoint_rhs(i).close();
882 0 : sensitivities[i][k] += (this->get_adjoint_rhs(i).dot(this->get_weighted_sensitivity_solution()) -
883 0 : this->rhs->dot(this->get_weighted_sensitivity_adjoint_solution(i)) -
884 0 : this->get_adjoint_solution(i).dot(*tempvec)) / (2.*delta_p);
885 : }
886 :
887 0 : *parameters_vec[k] = old_parameter - delta_p;
888 0 : this->assembly(true, true);
889 0 : this->rhs->close();
890 0 : this->matrix->close();
891 0 : this->assemble_qoi_derivative(qoi_indices,
892 : /* include_liftfunc = */ true,
893 0 : /* apply_constraints = */ false);
894 :
895 0 : this->matrix->vector_mult(*tempvec, this->get_weighted_sensitivity_solution());
896 :
897 0 : for (unsigned int i=0; i != Nq; ++i)
898 0 : if (qoi_indices.has_index(i))
899 : {
900 0 : this->get_adjoint_rhs(i).close();
901 0 : sensitivities[i][k] += (-this->get_adjoint_rhs(i).dot(this->get_weighted_sensitivity_solution()) +
902 0 : this->rhs->dot(this->get_weighted_sensitivity_adjoint_solution(i)) +
903 0 : this->get_adjoint_solution(i).dot(*tempvec)) / (2.*delta_p);
904 : }
905 : }
906 :
907 : // All parameters have been reset.
908 : // Don't leave the qoi or system changed - principle of least
909 : // surprise.
910 0 : this->assembly(true, true);
911 0 : this->rhs->close();
912 0 : this->matrix->close();
913 0 : this->assemble_qoi(qoi_indices);
914 0 : }
915 :
916 :
917 :
918 0 : void ImplicitSystem::qoi_parameter_hessian (const QoISet & qoi_indices,
919 : const ParameterVector & parameters_in,
920 : SensitivityData & sensitivities)
921 : {
922 : // We currently get partial derivatives via finite differencing
923 0 : const Real delta_p = TOLERANCE;
924 :
925 0 : ParameterVector & parameters_vec =
926 : const_cast<ParameterVector &>(parameters_in);
927 :
928 : // We'll use one temporary vector for matrix-vector-vector products
929 0 : std::unique_ptr<NumericVector<Number>> tempvec = this->solution->zero_clone();
930 :
931 : // And another temporary vector to hold a copy of the true solution
932 : // so we can safely perturb this->solution.
933 0 : std::unique_ptr<NumericVector<Number>> oldsolution = this->solution->clone();
934 :
935 : const unsigned int Np = cast_int<unsigned int>
936 0 : (parameters_vec.size());
937 0 : const unsigned int Nq = this->n_qois();
938 :
939 : // For each quantity of interest q, the parameter sensitivity
940 : // Hessian is defined as q''_{kl} = {d^2 q}/{d p_k d p_l}.
941 : //
942 : // We calculate it from values and partial derivatives of the
943 : // quantity of interest function Q, solution u, adjoint solution z,
944 : // and residual R, as:
945 : //
946 : // q''_{kl} =
947 : // Q''_{kl} + Q''_{uk}(u)*u'_l + Q''_{ul}(u) * u'_k +
948 : // Q''_{uu}(u)*u'_k*u'_l -
949 : // R''_{kl}(u,z) -
950 : // R''_{uk}(u,z)*u'_l - R''_{ul}(u,z)*u'_k -
951 : // R''_{uu}(u,z)*u'_k*u'_l
952 : //
953 : // See the adjoints model document for more details.
954 :
955 : // We first do an adjoint solve to get z for each quantity of
956 : // interest
957 : // if we haven't already or dont have an initial condition for the adjoint
958 0 : if (!this->is_adjoint_already_solved())
959 : {
960 0 : this->adjoint_solve(qoi_indices);
961 : }
962 :
963 : // And a sensitivity solve to get u_k for each parameter
964 0 : this->sensitivity_solve(parameters_vec);
965 :
966 : // Get ready to fill in second derivatives:
967 0 : sensitivities.allocate_hessian_data(qoi_indices, *this, parameters_vec);
968 :
969 0 : for (unsigned int k=0; k != Np; ++k)
970 : {
971 0 : Number old_parameterk = *parameters_vec[k];
972 :
973 : // The Hessian is symmetric, so we just calculate the lower
974 : // triangle and the diagonal, and we get the upper triangle from
975 : // the transpose of the lower
976 :
977 0 : for (unsigned int l=0; l != k+1; ++l)
978 : {
979 : // The second partial derivatives with respect to parameters_vec
980 : // are all calculated via a central finite difference
981 : // stencil:
982 : // F''_{kl} ~= (F(p+dp*e_k+dp*e_l) - F(p+dp*e_k-dp*e_l) -
983 : // F(p-dp*e_k+dp*e_l) + F(p-dp*e_k-dp*e_l))/(4*dp^2)
984 : // We will add Q''_{kl}(u) and subtract R''_{kl}(u,z) at the
985 : // same time.
986 : //
987 : // We have to be careful with the perturbations to handle
988 : // the k=l case
989 :
990 0 : Number old_parameterl = *parameters_vec[l];
991 :
992 0 : *parameters_vec[k] += delta_p;
993 0 : *parameters_vec[l] += delta_p;
994 0 : this->assemble_qoi(qoi_indices);
995 0 : this->assembly(true, false, true);
996 0 : this->rhs->close();
997 0 : std::vector<Number> partial2q_term = this->get_qoi_values();
998 0 : std::vector<Number> partial2R_term(this->n_qois());
999 0 : for (unsigned int i=0; i != Nq; ++i)
1000 0 : if (qoi_indices.has_index(i))
1001 0 : partial2R_term[i] = this->rhs->dot(this->get_adjoint_solution(i));
1002 :
1003 0 : *parameters_vec[l] -= 2.*delta_p;
1004 0 : this->assemble_qoi(qoi_indices);
1005 0 : this->assembly(true, false, true);
1006 0 : this->rhs->close();
1007 0 : for (unsigned int i=0; i != Nq; ++i)
1008 0 : if (qoi_indices.has_index(i))
1009 : {
1010 0 : partial2q_term[i] -= this->get_qoi_value(i);
1011 0 : partial2R_term[i] -= this->rhs->dot(this->get_adjoint_solution(i));
1012 : }
1013 :
1014 0 : *parameters_vec[k] -= 2.*delta_p;
1015 0 : this->assemble_qoi(qoi_indices);
1016 0 : this->assembly(true, false, true);
1017 0 : this->rhs->close();
1018 0 : for (unsigned int i=0; i != Nq; ++i)
1019 0 : if (qoi_indices.has_index(i))
1020 : {
1021 0 : partial2q_term[i] += this->get_qoi_value(i);
1022 0 : partial2R_term[i] += this->rhs->dot(this->get_adjoint_solution(i));
1023 : }
1024 :
1025 0 : *parameters_vec[l] += 2.*delta_p;
1026 0 : this->assemble_qoi(qoi_indices);
1027 0 : this->assembly(true, false, true);
1028 0 : this->rhs->close();
1029 0 : for (unsigned int i=0; i != Nq; ++i)
1030 0 : if (qoi_indices.has_index(i))
1031 : {
1032 0 : partial2q_term[i] -= this->get_qoi_value(i);
1033 0 : partial2R_term[i] -= this->rhs->dot(this->get_adjoint_solution(i));
1034 0 : partial2q_term[i] /= (4. * delta_p * delta_p);
1035 0 : partial2R_term[i] /= (4. * delta_p * delta_p);
1036 : }
1037 :
1038 0 : for (unsigned int i=0; i != Nq; ++i)
1039 0 : if (qoi_indices.has_index(i))
1040 : {
1041 0 : Number current_terms = partial2q_term[i] - partial2R_term[i];
1042 0 : sensitivities.second_derivative(i,k,l) += current_terms;
1043 0 : if (k != l)
1044 0 : sensitivities.second_derivative(i,l,k) += current_terms;
1045 : }
1046 :
1047 : // Don't leave the parameters_vec perturbed
1048 0 : *parameters_vec[l] = old_parameterl;
1049 0 : *parameters_vec[k] = old_parameterk;
1050 : }
1051 :
1052 : // We get (partial q / partial u) and
1053 : // (partial R / partial u) from the user, but centrally
1054 : // difference to get q_uk and R_uk terms:
1055 : // (partial^2 q / partial u partial k)
1056 : // q_uk*u'_l = (q_u(p+dp*e_k)*u'_l - q_u(p-dp*e_k)*u'_l)/(2*dp)
1057 : // R_uk*z*u'_l = (R_u(p+dp*e_k)*z*u'_l - R_u(p-dp*e_k)*z*u'_l)/(2*dp)
1058 : //
1059 : // To avoid creating Nq temporary vectors, we add these
1060 : // subterms to the sensitivities output one by one.
1061 : //
1062 : // FIXME: this is probably a bad order of operations for
1063 : // controlling floating point error.
1064 :
1065 0 : *parameters_vec[k] = old_parameterk + delta_p;
1066 0 : this->assembly(false, true);
1067 0 : this->matrix->close();
1068 0 : this->assemble_qoi_derivative(qoi_indices,
1069 : /* include_liftfunc = */ true,
1070 0 : /* apply_constraints = */ false);
1071 :
1072 0 : for (unsigned int l=0; l != Np; ++l)
1073 : {
1074 0 : this->matrix->vector_mult(*tempvec, this->get_sensitivity_solution(l));
1075 0 : for (unsigned int i=0; i != Nq; ++i)
1076 0 : if (qoi_indices.has_index(i))
1077 : {
1078 0 : this->get_adjoint_rhs(i).close();
1079 : Number current_terms =
1080 0 : (this->get_adjoint_rhs(i).dot(this->get_sensitivity_solution(l)) -
1081 0 : tempvec->dot(this->get_adjoint_solution(i))) / (2.*delta_p);
1082 0 : sensitivities.second_derivative(i,k,l) += current_terms;
1083 :
1084 : // We use the _uk terms twice; symmetry lets us reuse
1085 : // these calculations for the _ul terms.
1086 :
1087 0 : sensitivities.second_derivative(i,l,k) += current_terms;
1088 : }
1089 : }
1090 :
1091 0 : *parameters_vec[k] = old_parameterk - delta_p;
1092 0 : this->assembly(false, true);
1093 0 : this->matrix->close();
1094 0 : this->assemble_qoi_derivative(qoi_indices,
1095 : /* include_liftfunc = */ true,
1096 0 : /* apply_constraints = */ false);
1097 :
1098 0 : for (unsigned int l=0; l != Np; ++l)
1099 : {
1100 0 : this->matrix->vector_mult(*tempvec, this->get_sensitivity_solution(l));
1101 0 : for (unsigned int i=0; i != Nq; ++i)
1102 0 : if (qoi_indices.has_index(i))
1103 : {
1104 0 : this->get_adjoint_rhs(i).close();
1105 : Number current_terms =
1106 0 : (-this->get_adjoint_rhs(i).dot(this->get_sensitivity_solution(l)) +
1107 0 : tempvec->dot(this->get_adjoint_solution(i))) / (2.*delta_p);
1108 0 : sensitivities.second_derivative(i,k,l) += current_terms;
1109 :
1110 : // We use the _uk terms twice; symmetry lets us reuse
1111 : // these calculations for the _ul terms.
1112 :
1113 0 : sensitivities.second_derivative(i,l,k) += current_terms;
1114 : }
1115 : }
1116 :
1117 : // Don't leave the parameter perturbed
1118 0 : *parameters_vec[k] = old_parameterk;
1119 :
1120 : // Our last remaining terms are -R_uu(u,z)*u_k*u_l and
1121 : // Q_uu(u)*u_k*u_l
1122 : //
1123 : // We take directional central finite differences of R_u and Q_u
1124 : // to approximate these terms, e.g.:
1125 : //
1126 : // Q_uu(u)*u_k ~= (Q_u(u+dp*u_k) - Q_u(u-dp*u_k))/(2*dp)
1127 :
1128 0 : *this->solution = this->get_sensitivity_solution(k);
1129 0 : *this->solution *= delta_p;
1130 0 : *this->solution += *oldsolution;
1131 :
1132 : // We've modified solution, so we need to update before calling
1133 : // assembly since assembly may only use current_local_solution
1134 0 : this->update();
1135 0 : this->assembly(false, true);
1136 0 : this->matrix->close();
1137 0 : this->assemble_qoi_derivative(qoi_indices,
1138 : /* include_liftfunc = */ true,
1139 0 : /* apply_constraints = */ false);
1140 :
1141 : // The Hessian is symmetric, so we just calculate the lower
1142 : // triangle and the diagonal, and we get the upper triangle from
1143 : // the transpose of the lower
1144 : //
1145 : // Note that, because we took the directional finite difference
1146 : // with respect to k and not l, we've added an O(delta_p^2)
1147 : // error to any permutational symmetry in the Hessian...
1148 0 : for (unsigned int l=0; l != k+1; ++l)
1149 : {
1150 0 : this->matrix->vector_mult(*tempvec, this->get_sensitivity_solution(l));
1151 0 : for (unsigned int i=0; i != Nq; ++i)
1152 0 : if (qoi_indices.has_index(i))
1153 : {
1154 0 : this->get_adjoint_rhs(i).close();
1155 : Number current_terms =
1156 0 : (this->get_adjoint_rhs(i).dot(this->get_sensitivity_solution(l)) -
1157 0 : tempvec->dot(this->get_adjoint_solution(i))) / (2.*delta_p);
1158 0 : sensitivities.second_derivative(i,k,l) += current_terms;
1159 0 : if (k != l)
1160 0 : sensitivities.second_derivative(i,l,k) += current_terms;
1161 : }
1162 : }
1163 :
1164 0 : *this->solution = this->get_sensitivity_solution(k);
1165 0 : *this->solution *= -delta_p;
1166 0 : *this->solution += *oldsolution;
1167 :
1168 : // We've modified solution, so we need to update before calling
1169 : // assembly since assembly may only use current_local_solution
1170 0 : this->update();
1171 0 : this->assembly(false, true);
1172 0 : this->matrix->close();
1173 0 : this->assemble_qoi_derivative(qoi_indices,
1174 : /* include_liftfunc = */ true,
1175 0 : /* apply_constraints = */ false);
1176 :
1177 0 : for (unsigned int l=0; l != k+1; ++l)
1178 : {
1179 0 : this->matrix->vector_mult(*tempvec, this->get_sensitivity_solution(l));
1180 0 : for (unsigned int i=0; i != Nq; ++i)
1181 0 : if (qoi_indices.has_index(i))
1182 : {
1183 0 : this->get_adjoint_rhs(i).close();
1184 : Number current_terms =
1185 0 : (-this->get_adjoint_rhs(i).dot(this->get_sensitivity_solution(l)) +
1186 0 : tempvec->dot(this->get_adjoint_solution(i))) / (2.*delta_p);
1187 0 : sensitivities.second_derivative(i,k,l) += current_terms;
1188 0 : if (k != l)
1189 0 : sensitivities.second_derivative(i,l,k) += current_terms;
1190 : }
1191 : }
1192 :
1193 : // Don't leave the solution perturbed
1194 0 : *this->solution = *oldsolution;
1195 : }
1196 :
1197 : // All parameters_vec have been reset.
1198 : // Don't leave the qoi or system changed - principle of least
1199 : // surprise.
1200 : // We've modified solution, so we need to update before calling
1201 : // assembly since assembly may only use current_local_solution
1202 0 : this->update();
1203 0 : this->assembly(true, true);
1204 0 : this->rhs->close();
1205 0 : this->matrix->close();
1206 0 : this->assemble_qoi(qoi_indices);
1207 0 : }
1208 :
1209 :
1210 :
1211 0 : LinearSolver<Number> * ImplicitSystem::get_linear_solver() const
1212 : {
1213 : // Note: we always start "from scratch" to mimic the original
1214 : // behavior of this function. The goal is not to reuse the
1215 : // LinearSolver object, but to manage its lifetime in a more
1216 : // consistent manner.
1217 0 : linear_solver.reset();
1218 :
1219 0 : linear_solver = LinearSolver<Number>::build(this->comm());
1220 :
1221 0 : if (this->prefix_with_name())
1222 0 : linear_solver->init(this->prefix().c_str());
1223 : else
1224 0 : linear_solver->init();
1225 :
1226 0 : linear_solver->init_names(*this);
1227 :
1228 0 : return linear_solver.get();
1229 : }
1230 :
1231 :
1232 :
1233 231571 : std::pair<unsigned int, Real> ImplicitSystem::get_linear_solve_parameters() const
1234 : {
1235 218007 : return std::make_pair(
1236 231571 : parameters.have_parameter<unsigned int>("linear solver maximum iterations")
1237 71 : ? parameters.get<unsigned int>("linear solver maximum iterations")
1238 231500 : : this->get_equation_systems().parameters.get<unsigned int>(
1239 13560 : "linear solver maximum iterations"),
1240 231571 : parameters.have_parameter<Real>("linear solver tolerance")
1241 71 : ? parameters.get<Real>("linear solver tolerance")
1242 476639 : : this->get_equation_systems().parameters.get<Real>("linear solver tolerance"));
1243 : }
1244 :
1245 :
1246 :
1247 : #ifdef LIBMESH_ENABLE_DEPRECATED
1248 0 : void ImplicitSystem::release_linear_solver(LinearSolver<Number> *) const
1249 : {
1250 : // This function was originally paired with get_linear_solver()
1251 : // calls when that returned a dumb pointer which needed to be
1252 : // cleaned up. Since get_linear_solver() now just returns a pointer
1253 : // to a LinearSolver object managed by this class, this function no
1254 : // longer needs to do any cleanup.
1255 : libmesh_deprecated();
1256 0 : }
1257 : #endif // LIBMESH_ENABLE_DEPRECATED
1258 :
1259 :
1260 :
1261 0 : const SparseMatrix<Number> & ImplicitSystem::get_system_matrix() const
1262 : {
1263 0 : libmesh_assert(matrix);
1264 0 : libmesh_assert_equal_to(&get_matrix("System Matrix"), matrix);
1265 0 : return *matrix;
1266 : }
1267 :
1268 :
1269 :
1270 6752969 : SparseMatrix<Number> & ImplicitSystem::get_system_matrix()
1271 : {
1272 595194 : libmesh_assert(matrix);
1273 595194 : libmesh_assert_equal_to(&get_matrix("System Matrix"), matrix);
1274 6752969 : return *matrix;
1275 : }
1276 :
1277 : template void ImplicitSystem::setup_static_condensation_preconditioner(LinearSolver<Number> & solver);
1278 : template void ImplicitSystem::setup_static_condensation_preconditioner(NonlinearSolver<Number> & solver);
1279 : } // namespace libMesh
|
