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