libMesh
implicit_system.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 // C++ includes
21 
22 // Local includes
23 #include "libmesh/dof_map.h"
24 #include "libmesh/equation_systems.h"
25 #include "libmesh/implicit_system.h"
26 #include "libmesh/libmesh_logging.h"
27 #include "libmesh/linear_solver.h"
28 #include "libmesh/mesh_base.h"
29 #include "libmesh/numeric_vector.h"
30 #include "libmesh/parameters.h"
31 #include "libmesh/parameter_vector.h"
32 #include "libmesh/qoi_set.h"
33 #include "libmesh/sensitivity_data.h"
34 #include "libmesh/sparse_matrix.h"
35 
36 namespace libMesh
37 {
38 
39 // ------------------------------------------------------------
40 // ImplicitSystem implementation
42  const std::string & name_in,
43  const unsigned int number_in) :
44 
45  Parent (es, name_in, number_in),
46  matrix (libmesh_nullptr),
47  zero_out_matrix_and_rhs(true),
48  _can_add_matrices (true)
49 {
50  // Add the system matrix.
51  this->add_system_matrix ();
52 }
53 
54 
55 
57 {
58  // Clear data
59  this->clear();
60 
61  remove_matrix("System Matrix");
62 }
63 
64 
65 
67 {
68  // clear the parent data
69  Parent::clear();
70 
71  // clear any user-added matrices
72  {
73  for (matrices_iterator pos = _matrices.begin();
74  pos != _matrices.end(); ++pos)
75  {
76  pos->second->clear ();
77  delete pos->second;
78  pos->second = libmesh_nullptr;
79  }
80 
81  _matrices.clear();
82  _can_add_matrices = true;
83  }
84 
85  // Restore us to a "basic" state
86  this->add_system_matrix ();
87 }
88 
89 
90 
92 {
93  // initialize parent data
95 
96  // Clear any existing matrices
97  for (matrices_iterator pos = _matrices.begin();
98  pos != _matrices.end(); ++pos)
99  pos->second->clear();
100 
101  // Initialize the matrices for the system
102  this->init_matrices ();
103 }
104 
105 
106 
108 {
110 
111  // Check for quick return in case the system matrix
112  // (and by extension all the matrices) has already
113  // been initialized
114  if (matrix->initialized())
115  return;
116 
117  // Get a reference to the DofMap
118  DofMap & dof_map = this->get_dof_map();
119 
120  // no chance to add other matrices
121  _can_add_matrices = false;
122 
123  // Tell the matrices about the dof map, and vice versa
124  for (matrices_iterator pos = _matrices.begin();
125  pos != _matrices.end(); ++pos)
126  {
127  SparseMatrix<Number> & m = *(pos->second);
129 
130  // We want to allow repeated init() on systems, but we don't
131  // want to attach the same matrix to the DofMap twice
132  if (!dof_map.is_attached(m))
133  dof_map.attach_matrix (m);
134  }
135 
136  // Compute the sparsity pattern for the current
137  // mesh and DOF distribution. This also updates
138  // additional matrices, \p DofMap now knows them
139  dof_map.compute_sparsity (this->get_mesh());
140 
141  // Initialize matrices
142  for (matrices_iterator pos = _matrices.begin();
143  pos != _matrices.end(); ++pos)
144  pos->second->init ();
145 
146  // Set the additional matrices to 0.
147  for (matrices_iterator pos = _matrices.begin();
148  pos != _matrices.end(); ++pos)
149  pos->second->zero ();
150 }
151 
152 
153 
155 {
156  // initialize parent data
157  Parent::reinit();
158 
159  // Get a reference to the DofMap
160  DofMap & dof_map = this->get_dof_map();
161 
162  // Clear the matrices
163  for (matrices_iterator pos = _matrices.begin();
164  pos != _matrices.end(); ++pos)
165  {
166  pos->second->clear();
167  pos->second->attach_dof_map (dof_map);
168  }
169 
170  // Clear the sparsity pattern
171  this->get_dof_map().clear_sparsity();
172 
173  // Compute the sparsity pattern for the current
174  // mesh and DOF distribution. This also updates
175  // additional matrices, \p DofMap now knows them
176  dof_map.compute_sparsity (this->get_mesh());
177 
178  // Initialize matrices
179  for (matrices_iterator pos = _matrices.begin();
180  pos != _matrices.end(); ++pos)
181  pos->second->init ();
182 
183  // Set the additional matrices to 0.
184  for (matrices_iterator pos = _matrices.begin();
185  pos != _matrices.end(); ++pos)
186  pos->second->zero ();
187 }
188 
189 
190 
192 {
197 
199  {
200  matrix->zero ();
201  rhs->zero ();
202  }
203 
204  // Call the base class assemble function
205  Parent::assemble ();
206 }
207 
208 
209 
210 SparseMatrix<Number> & ImplicitSystem::add_matrix (const std::string & mat_name)
211 {
212  // only add matrices before initializing...
213  if (!_can_add_matrices)
214  libmesh_error_msg("ERROR: Too late. Cannot add matrices to the system after initialization"
215  << "\n any more. You should have done this earlier.");
216 
217  // Return the matrix if it is already there.
218  if (this->have_matrix(mat_name))
219  return *(_matrices[mat_name]);
220 
221  // Otherwise build the matrix and return it.
222  SparseMatrix<Number> * buf = SparseMatrix<Number>::build(this->comm()).release();
223  _matrices.insert (std::make_pair (mat_name, buf));
224 
225  return *buf;
226 }
227 
228 
229 void ImplicitSystem::remove_matrix (const std::string & mat_name)
230 {
231  matrices_iterator pos = _matrices.find (mat_name);
232 
233  //Return if the matrix does not exist
234  if (pos == _matrices.end())
235  return;
236 
237  delete pos->second;
238 
239  _matrices.erase(pos);
240 }
241 
242 
243 
244 const SparseMatrix<Number> * ImplicitSystem::request_matrix (const std::string & mat_name) const
245 {
246  // Make sure the matrix exists
247  const_matrices_iterator pos = _matrices.find (mat_name);
248 
249  if (pos == _matrices.end())
250  return libmesh_nullptr;
251 
252  return pos->second;
253 }
254 
255 
256 
257 SparseMatrix<Number> * ImplicitSystem::request_matrix (const std::string & mat_name)
258 {
259  // Make sure the matrix exists
260  matrices_iterator pos = _matrices.find (mat_name);
261 
262  if (pos == _matrices.end())
263  return libmesh_nullptr;
264 
265  return pos->second;
266 }
267 
268 
269 
270 const SparseMatrix<Number> & ImplicitSystem::get_matrix (const std::string & mat_name) const
271 {
272  // Make sure the matrix exists
273  const_matrices_iterator pos = _matrices.find (mat_name);
274 
275  if (pos == _matrices.end())
276  libmesh_error_msg("ERROR: matrix " << mat_name << " does not exist in this system!");
277 
278  return *(pos->second);
279 }
280 
281 
282 
283 SparseMatrix<Number> & ImplicitSystem::get_matrix (const std::string & mat_name)
284 {
285  // Make sure the matrix exists
286  matrices_iterator pos = _matrices.find (mat_name);
287 
288  if (pos == _matrices.end())
289  libmesh_error_msg("ERROR: matrix " << mat_name << " does not exist in this system!");
290 
291  return *(pos->second);
292 }
293 
294 
295 
297 {
298  // Possible that we cleared the _matrices but
299  // forgot to NULL-out the matrix?
300  if (_matrices.empty()) matrix = libmesh_nullptr;
301 
302 
303  // Only need to add the matrix if it isn't there
304  // already!
305  if (matrix == libmesh_nullptr)
306  matrix = &(this->add_matrix ("System Matrix"));
307 
309 }
310 
311 
312 
314  this->assemble_before_solve = true;
315  this->get_linear_solver()->reuse_preconditioner(false);
316 }
317 
318 
319 
320 std::pair<unsigned int, Real>
322 {
323  // Log how long the linear solve takes.
324  LOG_SCOPE("sensitivity_solve()", "ImplicitSystem");
325 
326  // The forward system should now already be solved.
327  // Now assemble the corresponding sensitivity system.
328 
329  if (this->assemble_before_solve)
330  {
331  // Build the Jacobian
332  this->assembly(false, true);
333  this->matrix->close();
334 
335  // Reset and build the RHS from the residual derivatives
336  this->assemble_residual_derivatives(parameters);
337  }
338 
339  // The sensitivity problem is linear
340  LinearSolver<Number> * linear_solver = this->get_linear_solver();
341 
342  // Our iteration counts and residuals will be sums of the individual
343  // results
344  std::pair<unsigned int, Real> solver_params =
346  std::pair<unsigned int, Real> totalrval = std::make_pair(0,0.0);
347 
348  // Solve the linear system.
349  SparseMatrix<Number> * pc = this->request_matrix("Preconditioner");
350  for (std::size_t p=0; p != parameters.size(); ++p)
351  {
352  std::pair<unsigned int, Real> rval =
353  linear_solver->solve (*matrix, pc,
354  this->add_sensitivity_solution(p),
355  this->get_sensitivity_rhs(p),
356  solver_params.second,
357  solver_params.first);
358 
359  totalrval.first += rval.first;
360  totalrval.second += rval.second;
361  }
362 
363  // The linear solver may not have fit our constraints exactly
364 #ifdef LIBMESH_ENABLE_CONSTRAINTS
365  for (std::size_t p=0; p != parameters.size(); ++p)
367  (*this, &this->get_sensitivity_solution(p),
368  /* homogeneous = */ true);
369 #endif
370 
371  this->release_linear_solver(linear_solver);
372 
373  return totalrval;
374 }
375 
376 
377 
378 std::pair<unsigned int, Real>
380 {
381  // Log how long the linear solve takes.
382  LOG_SCOPE("adjoint_solve()", "ImplicitSystem");
383 
384  if (this->assemble_before_solve)
385  // Assemble the linear system
386  this->assembly (/* get_residual = */ false,
387  /* get_jacobian = */ true);
388 
389  // The adjoint problem is linear
390  LinearSolver<Number> * linear_solver = this->get_linear_solver();
391 
392  // Reset and build the RHS from the QOI derivative
393  this->assemble_qoi_derivative(qoi_indices,
394  /* include_liftfunc = */ false,
395  /* apply_constraints = */ true);
396 
397  // Our iteration counts and residuals will be sums of the individual
398  // results
399  std::pair<unsigned int, Real> solver_params =
401  std::pair<unsigned int, Real> totalrval = std::make_pair(0,0.0);
402 
403  for (std::size_t i=0; i != this->qoi.size(); ++i)
404  if (qoi_indices.has_index(i))
405  {
406  const std::pair<unsigned int, Real> rval =
407  linear_solver->adjoint_solve (*matrix, this->add_adjoint_solution(i),
408  this->get_adjoint_rhs(i),
409  solver_params.second,
410  solver_params.first);
411 
412  totalrval.first += rval.first;
413  totalrval.second += rval.second;
414  }
415 
416  this->release_linear_solver(linear_solver);
417 
418  // The linear solver may not have fit our constraints exactly
419 #ifdef LIBMESH_ENABLE_CONSTRAINTS
420  for (std::size_t i=0; i != this->qoi.size(); ++i)
421  if (qoi_indices.has_index(i))
423  (this->get_adjoint_solution(i), i);
424 #endif
425 
426  return totalrval;
427 }
428 
429 
430 
431 std::pair<unsigned int, Real>
433  const ParameterVector & weights,
434  const QoISet & qoi_indices)
435 {
436  // Log how long the linear solve takes.
437  LOG_SCOPE("weighted_sensitivity_adjoint_solve()", "ImplicitSystem");
438 
439  // We currently get partial derivatives via central differencing
440  const Real delta_p = TOLERANCE;
441 
442  ParameterVector & parameters =
443  const_cast<ParameterVector &>(parameters_in);
444 
445  // The forward system should now already be solved.
446  // The adjoint system should now already be solved.
447  // Now we're assembling a weighted sum of adjoint-adjoint systems:
448  //
449  // dR/du (u, sum_l(w_l*z^l)) = sum_l(w_l*(Q''_ul - R''_ul (u, z)))
450 
451  // FIXME: The derivation here does not yet take adjoint boundary
452  // conditions into account.
453  for (std::size_t i=0; i != this->qoi.size(); ++i)
454  if (qoi_indices.has_index(i))
456 
457  // We'll assemble the rhs first, because the R'' term will require
458  // perturbing the jacobian
459 
460  // We'll use temporary rhs vectors, because we haven't (yet) found
461  // any good reasons why users might want to save these:
462 
463  std::vector<NumericVector<Number> *> temprhs(this->qoi.size());
464  for (std::size_t i=0; i != this->qoi.size(); ++i)
465  if (qoi_indices.has_index(i))
466  temprhs[i] = this->rhs->zero_clone().release();
467 
468  // We approximate the _l partial derivatives via a central
469  // differencing perturbation in the w_l direction:
470  //
471  // sum_l(w_l*v_l) ~= (v(p + dp*w_l*e_l) - v(p - dp*w_l*e_l))/(2*dp)
472 
473  // PETSc doesn't implement SGEMX, so neither does NumericVector,
474  // so we want to avoid calculating f -= R'*z. We'll thus evaluate
475  // the above equation by first adding -v(p+dp...), then multiplying
476  // the intermediate result vectors by -1, then adding -v(p-dp...),
477  // then finally dividing by 2*dp.
478 
479  ParameterVector oldparameters, parameterperturbation;
480  parameters.deep_copy(oldparameters);
481  weights.deep_copy(parameterperturbation);
482  parameterperturbation *= delta_p;
483  parameters += parameterperturbation;
484 
485  this->assembly(false, true);
486  this->matrix->close();
487 
488  // Take the discrete adjoint, so that we can calculate R_u(u,z) with
489  // a matrix-vector product of R_u and z.
491 
492  this->assemble_qoi_derivative(qoi_indices,
493  /* include_liftfunc = */ false,
494  /* apply_constraints = */ true);
495  for (std::size_t i=0; i != this->qoi.size(); ++i)
496  if (qoi_indices.has_index(i))
497  {
498  this->get_adjoint_rhs(i).close();
499  *(temprhs[i]) -= this->get_adjoint_rhs(i);
500  this->matrix->vector_mult_add(*(temprhs[i]), this->get_adjoint_solution(i));
501  *(temprhs[i]) *= -1.0;
502  }
503 
504  oldparameters.value_copy(parameters);
505  parameterperturbation *= -1.0;
506  parameters += parameterperturbation;
507 
508  this->assembly(false, true);
509  this->matrix->close();
511 
512  this->assemble_qoi_derivative(qoi_indices,
513  /* include_liftfunc = */ false,
514  /* apply_constraints = */ true);
515  for (std::size_t i=0; i != this->qoi.size(); ++i)
516  if (qoi_indices.has_index(i))
517  {
518  this->get_adjoint_rhs(i).close();
519  *(temprhs[i]) -= this->get_adjoint_rhs(i);
520  this->matrix->vector_mult_add(*(temprhs[i]), this->get_adjoint_solution(i));
521  *(temprhs[i]) /= (2.0*delta_p);
522  }
523 
524  // Finally, assemble the jacobian at the non-perturbed parameter
525  // values. Ignore assemble_before_solve; if we had a good
526  // non-perturbed matrix before we've already overwritten it.
527  oldparameters.value_copy(parameters);
528 
529  // if (this->assemble_before_solve)
530  {
531  // Build the Jacobian
532  this->assembly(false, true);
533  this->matrix->close();
534 
535  // Take the discrete adjoint
537  }
538 
539  // The weighted adjoint-adjoint problem is linear
540  LinearSolver<Number> * linear_solver = this->get_linear_solver();
541 
542  // Our iteration counts and residuals will be sums of the individual
543  // results
544  std::pair<unsigned int, Real> solver_params =
546  std::pair<unsigned int, Real> totalrval = std::make_pair(0,0.0);
547 
548  for (std::size_t i=0; i != this->qoi.size(); ++i)
549  if (qoi_indices.has_index(i))
550  {
551  const std::pair<unsigned int, Real> rval =
552  linear_solver->solve (*matrix, this->add_weighted_sensitivity_adjoint_solution(i),
553  *(temprhs[i]),
554  solver_params.second,
555  solver_params.first);
556 
557  totalrval.first += rval.first;
558  totalrval.second += rval.second;
559  }
560 
561  this->release_linear_solver(linear_solver);
562 
563  for (std::size_t i=0; i != this->qoi.size(); ++i)
564  if (qoi_indices.has_index(i))
565  delete temprhs[i];
566 
567  // The linear solver may not have fit our constraints exactly
568 #ifdef LIBMESH_ENABLE_CONSTRAINTS
569  for (std::size_t i=0; i != this->qoi.size(); ++i)
570  if (qoi_indices.has_index(i))
573  /* homogeneous = */ true);
574 #endif
575 
576  return totalrval;
577 }
578 
579 
580 
581 std::pair<unsigned int, Real>
583  const ParameterVector & weights)
584 {
585  // Log how long the linear solve takes.
586  LOG_SCOPE("weighted_sensitivity_solve()", "ImplicitSystem");
587 
588  // We currently get partial derivatives via central differencing
589  const Real delta_p = TOLERANCE;
590 
591  ParameterVector & parameters =
592  const_cast<ParameterVector &>(parameters_in);
593 
594  // The forward system should now already be solved.
595 
596  // Now we're assembling a weighted sum of sensitivity systems:
597  //
598  // dR/du (u, v)(sum(w_l*u'_l)) = -sum_l(w_l*R'_l (u, v)) forall v
599 
600  // We'll assemble the rhs first, because the R' term will require
601  // perturbing the system, and some applications may not be able to
602  // assemble a perturbed residual without simultaneously constructing
603  // a perturbed jacobian.
604 
605  // We approximate the _l partial derivatives via a central
606  // differencing perturbation in the w_l direction:
607  //
608  // sum_l(w_l*v_l) ~= (v(p + dp*w_l*e_l) - v(p - dp*w_l*e_l))/(2*dp)
609 
610  ParameterVector oldparameters, parameterperturbation;
611  parameters.deep_copy(oldparameters);
612  weights.deep_copy(parameterperturbation);
613  parameterperturbation *= delta_p;
614  parameters += parameterperturbation;
615 
616  this->assembly(true, false, true);
617  this->rhs->close();
618 
619  UniquePtr<NumericVector<Number>> temprhs = this->rhs->clone();
620 
621  oldparameters.value_copy(parameters);
622  parameterperturbation *= -1.0;
623  parameters += parameterperturbation;
624 
625  this->assembly(true, false, true);
626  this->rhs->close();
627 
628  *temprhs -= *(this->rhs);
629  *temprhs /= (2.0*delta_p);
630 
631  // Finally, assemble the jacobian at the non-perturbed parameter
632  // values
633  oldparameters.value_copy(parameters);
634 
635  // Build the Jacobian
636  this->assembly(false, true);
637  this->matrix->close();
638 
639  // The weighted sensitivity problem is linear
640  LinearSolver<Number> * linear_solver = this->get_linear_solver();
641 
642  std::pair<unsigned int, Real> solver_params =
644 
645  const std::pair<unsigned int, Real> rval =
646  linear_solver->solve (*matrix, this->add_weighted_sensitivity_solution(),
647  *temprhs,
648  solver_params.second,
649  solver_params.first);
650 
651  this->release_linear_solver(linear_solver);
652 
653  // The linear solver may not have fit our constraints exactly
654 #ifdef LIBMESH_ENABLE_CONSTRAINTS
656  (*this, &this->get_weighted_sensitivity_solution(),
657  /* homogeneous = */ true);
658 #endif
659 
660  return rval;
661 }
662 
663 
664 
666 {
667  ParameterVector & parameters =
668  const_cast<ParameterVector &>(parameters_in);
669 
670  const unsigned int Np = cast_int<unsigned int>
671  (parameters.size());
672 
673  for (unsigned int p=0; p != Np; ++p)
674  {
675  NumericVector<Number> & sensitivity_rhs = this->add_sensitivity_rhs(p);
676 
677  // Approximate -(partial R / partial p) by
678  // (R(p-dp) - R(p+dp)) / (2*dp)
679 
680  Number old_parameter = *parameters[p];
681 
682  const Real delta_p =
683  TOLERANCE * std::max(std::abs(old_parameter), 1e-3);
684 
685  *parameters[p] -= delta_p;
686 
687  // this->assembly(true, false, true);
688  this->assembly(true, false, false);
689  this->rhs->close();
690  sensitivity_rhs = *this->rhs;
691 
692  *parameters[p] = old_parameter + delta_p;
693 
694  // this->assembly(true, false, true);
695  this->assembly(true, false, false);
696  this->rhs->close();
697 
698  sensitivity_rhs -= *this->rhs;
699  sensitivity_rhs /= (2*delta_p);
700  sensitivity_rhs.close();
701 
702  *parameters[p] = old_parameter;
703  }
704 }
705 
706 
707 
709  const ParameterVector & parameters_in,
710  SensitivityData & sensitivities)
711 {
712  ParameterVector & parameters =
713  const_cast<ParameterVector &>(parameters_in);
714 
715  const unsigned int Np = cast_int<unsigned int>
716  (parameters.size());
717  const unsigned int Nq = cast_int<unsigned int>
718  (qoi.size());
719 
720  // An introduction to the problem:
721  //
722  // Residual R(u(p),p) = 0
723  // partial R / partial u = J = system matrix
724  //
725  // This implies that:
726  // d/dp(R) = 0
727  // (partial R / partial p) +
728  // (partial R / partial u) * (partial u / partial p) = 0
729 
730  // We first do an adjoint solve:
731  // J^T * z = (partial q / partial u)
732  // if we havent already or dont have an initial condition for the adjoint
733  if (!this->is_adjoint_already_solved())
734  {
735  this->adjoint_solve(qoi_indices);
736  }
737 
738  this->assemble_residual_derivatives(parameters_in);
739 
740  // Get ready to fill in sensitivities:
741  sensitivities.allocate_data(qoi_indices, *this, parameters);
742 
743  // We use the identities:
744  // dq/dp = (partial q / partial p) + (partial q / partial u) *
745  // (partial u / partial p)
746  // dq/dp = (partial q / partial p) + (J^T * z) *
747  // (partial u / partial p)
748  // dq/dp = (partial q / partial p) + z * J *
749  // (partial u / partial p)
750 
751  // Leading to our final formula:
752  // dq/dp = (partial q / partial p) - z * (partial R / partial p)
753 
754  // In the case of adjoints with heterogenous Dirichlet boundary
755  // function phi, where
756  // q := S(u) - R(u,phi)
757  // the final formula works out to:
758  // dq/dp = (partial S / partial p) - z * (partial R / partial p)
759  // Because we currently have no direct access to
760  // (partial S / partial p), we use the identity
761  // (partial S / partial p) = (partial q / partial p) +
762  // phi * (partial R / partial p)
763  // to derive an equivalent equation:
764  // dq/dp = (partial q / partial p) - (z-phi) * (partial R / partial p)
765 
766  // Since z-phi degrees of freedom are zero for constrained indices,
767  // we can use the same constrained -(partial R / partial p) that we
768  // use for forward sensitivity solves, taking into account the
769  // differing sign convention.
770  //
771  // Since that vector is constrained, its constrained indices are
772  // zero, so its product with phi is zero, so we can neglect the
773  // evaluation of phi terms.
774 
775  for (unsigned int j=0; j != Np; ++j)
776  {
777  // We currently get partial derivatives via central differencing
778 
779  // (partial q / partial p) ~= (q(p+dp)-q(p-dp))/(2*dp)
780  // (partial R / partial p) ~= (rhs(p+dp) - rhs(p-dp))/(2*dp)
781 
782  Number old_parameter = *parameters[j];
783 
784  const Real delta_p =
785  TOLERANCE * std::max(std::abs(old_parameter), 1e-3);
786 
787  *parameters[j] = old_parameter - delta_p;
788  this->assemble_qoi(qoi_indices);
789  std::vector<Number> qoi_minus = this->qoi;
790 
791  NumericVector<Number> & neg_partialR_partialp = this->get_sensitivity_rhs(j);
792 
793  *parameters[j] = old_parameter + delta_p;
794  this->assemble_qoi(qoi_indices);
795  std::vector<Number> & qoi_plus = this->qoi;
796 
797  std::vector<Number> partialq_partialp(Nq, 0);
798  for (unsigned int i=0; i != Nq; ++i)
799  if (qoi_indices.has_index(i))
800  partialq_partialp[i] = (qoi_plus[i] - qoi_minus[i]) / (2.*delta_p);
801 
802  // Don't leave the parameter changed
803  *parameters[j] = old_parameter;
804 
805  for (unsigned int i=0; i != Nq; ++i)
806  if (qoi_indices.has_index(i))
807  sensitivities[i][j] = partialq_partialp[i] +
808  neg_partialR_partialp.dot(this->get_adjoint_solution(i));
809  }
810 
811  // All parameters have been reset.
812  // Reset the original qoi.
813 
814  this->assemble_qoi(qoi_indices);
815 }
816 
817 
818 
820  const ParameterVector & parameters_in,
821  SensitivityData & sensitivities)
822 {
823  ParameterVector & parameters =
824  const_cast<ParameterVector &>(parameters_in);
825 
826  const unsigned int Np = cast_int<unsigned int>
827  (parameters.size());
828  const unsigned int Nq = cast_int<unsigned int>
829  (qoi.size());
830 
831  // An introduction to the problem:
832  //
833  // Residual R(u(p),p) = 0
834  // partial R / partial u = J = system matrix
835  //
836  // This implies that:
837  // d/dp(R) = 0
838  // (partial R / partial p) +
839  // (partial R / partial u) * (partial u / partial p) = 0
840 
841  // We first solve for (partial u / partial p) for each parameter:
842  // J * (partial u / partial p) = - (partial R / partial p)
843 
844  this->sensitivity_solve(parameters);
845 
846  // Get ready to fill in sensitivities:
847  sensitivities.allocate_data(qoi_indices, *this, parameters);
848 
849  // We use the identity:
850  // dq/dp = (partial q / partial p) + (partial q / partial u) *
851  // (partial u / partial p)
852 
853  // We get (partial q / partial u) from the user
854  this->assemble_qoi_derivative(qoi_indices,
855  /* include_liftfunc = */ true,
856  /* apply_constraints = */ false);
857 
858  // We don't need these to be closed() in this function, but libMesh
859  // standard practice is to have them closed() by the time the
860  // function exits
861  for (std::size_t i=0; i != this->qoi.size(); ++i)
862  if (qoi_indices.has_index(i))
863  this->get_adjoint_rhs(i).close();
864 
865  for (unsigned int j=0; j != Np; ++j)
866  {
867  // We currently get partial derivatives via central differencing
868 
869  // (partial q / partial p) ~= (q(p+dp)-q(p-dp))/(2*dp)
870 
871  Number old_parameter = *parameters[j];
872 
873  const Real delta_p =
874  TOLERANCE * std::max(std::abs(old_parameter), 1e-3);
875 
876  *parameters[j] = old_parameter - delta_p;
877  this->assemble_qoi(qoi_indices);
878  std::vector<Number> qoi_minus = this->qoi;
879 
880  *parameters[j] = old_parameter + delta_p;
881  this->assemble_qoi(qoi_indices);
882  std::vector<Number> & qoi_plus = this->qoi;
883 
884  std::vector<Number> partialq_partialp(Nq, 0);
885  for (unsigned int i=0; i != Nq; ++i)
886  if (qoi_indices.has_index(i))
887  partialq_partialp[i] = (qoi_plus[i] - qoi_minus[i]) / (2.*delta_p);
888 
889  // Don't leave the parameter changed
890  *parameters[j] = old_parameter;
891 
892  for (unsigned int i=0; i != Nq; ++i)
893  if (qoi_indices.has_index(i))
894  sensitivities[i][j] = partialq_partialp[i] +
895  this->get_adjoint_rhs(i).dot(this->get_sensitivity_solution(j));
896  }
897 
898  // All parameters have been reset.
899  // We didn't cache the original rhs or matrix for memory reasons,
900  // but we can restore them to a state consistent solution -
901  // principle of least surprise.
902  this->assembly(true, true);
903  this->rhs->close();
904  this->matrix->close();
905  this->assemble_qoi(qoi_indices);
906 }
907 
908 
909 
911  const ParameterVector & parameters_in,
912  const ParameterVector & vector,
913  SensitivityData & sensitivities)
914 {
915  // We currently get partial derivatives via finite differencing
916  const Real delta_p = TOLERANCE;
917 
918  ParameterVector & parameters =
919  const_cast<ParameterVector &>(parameters_in);
920 
921  // We'll use a single temporary vector for matrix-vector-vector products
922  UniquePtr<NumericVector<Number>> tempvec = this->solution->zero_clone();
923 
924  const unsigned int Np = cast_int<unsigned int>
925  (parameters.size());
926  const unsigned int Nq = cast_int<unsigned int>
927  (qoi.size());
928 
929  // For each quantity of interest q, the parameter sensitivity
930  // Hessian is defined as q''_{kl} = {d^2 q}/{d p_k d p_l}.
931  // Given a vector of parameter perturbation weights w_l, this
932  // function evaluates the hessian-vector product sum_l(q''_{kl}*w_l)
933  //
934  // We calculate it from values and partial derivatives of the
935  // quantity of interest function Q, solution u, adjoint solution z,
936  // parameter sensitivity adjoint solutions z^l, and residual R, as:
937  //
938  // sum_l(q''_{kl}*w_l) =
939  // sum_l(w_l * Q''_{kl}) + Q''_{uk}(u)*(sum_l(w_l u'_l)) -
940  // R'_k(u, sum_l(w_l*z^l)) - R'_{uk}(u,z)*(sum_l(w_l u'_l) -
941  // sum_l(w_l*R''_{kl}(u,z))
942  //
943  // See the adjoints model document for more details.
944 
945  // We first do an adjoint solve to get z for each quantity of
946  // interest
947  // if we havent already or dont have an initial condition for the adjoint
948  if (!this->is_adjoint_already_solved())
949  {
950  this->adjoint_solve(qoi_indices);
951  }
952 
953  // Get ready to fill in sensitivities:
954  sensitivities.allocate_data(qoi_indices, *this, parameters);
955 
956  // We can't solve for all the solution sensitivities u'_l or for all
957  // of the parameter sensitivity adjoint solutions z^l without
958  // requiring O(Nq*Np) linear solves. So we'll solve directly for their
959  // weighted sum - this is just O(Nq) solves.
960 
961  // First solve for sum_l(w_l u'_l).
962  this->weighted_sensitivity_solve(parameters, vector);
963 
964  // Then solve for sum_l(w_l z^l).
965  this->weighted_sensitivity_adjoint_solve(parameters, vector, qoi_indices);
966 
967  for (unsigned int k=0; k != Np; ++k)
968  {
969  // We approximate sum_l(w_l * Q''_{kl}) with a central
970  // differencing perturbation:
971  // sum_l(w_l * Q''_{kl}) ~=
972  // (Q(p + dp*w_l*e_l + dp*e_k) - Q(p - dp*w_l*e_l + dp*e_k) -
973  // Q(p + dp*w_l*e_l - dp*e_k) + Q(p - dp*w_l*e_l - dp*e_k))/(4*dp^2)
974 
975  // The sum(w_l*R''_kl) term requires the same sort of perturbation,
976  // and so we subtract it in at the same time:
977  // sum_l(w_l * R''_{kl}) ~=
978  // (R(p + dp*w_l*e_l + dp*e_k) - R(p - dp*w_l*e_l + dp*e_k) -
979  // R(p + dp*w_l*e_l - dp*e_k) + R(p - dp*w_l*e_l - dp*e_k))/(4*dp^2)
980 
981  ParameterVector oldparameters, parameterperturbation;
982  parameters.deep_copy(oldparameters);
983  vector.deep_copy(parameterperturbation);
984  parameterperturbation *= delta_p;
985  parameters += parameterperturbation;
986 
987  Number old_parameter = *parameters[k];
988 
989  *parameters[k] = old_parameter + delta_p;
990  this->assemble_qoi(qoi_indices);
991  this->assembly(true, false, true);
992  this->rhs->close();
993  std::vector<Number> partial2q_term = this->qoi;
994  std::vector<Number> partial2R_term(this->qoi.size());
995  for (unsigned int i=0; i != Nq; ++i)
996  if (qoi_indices.has_index(i))
997  partial2R_term[i] = this->rhs->dot(this->get_adjoint_solution(i));
998 
999  *parameters[k] = old_parameter - delta_p;
1000  this->assemble_qoi(qoi_indices);
1001  this->assembly(true, false, true);
1002  this->rhs->close();
1003  for (unsigned int i=0; i != Nq; ++i)
1004  if (qoi_indices.has_index(i))
1005  {
1006  partial2q_term[i] -= this->qoi[i];
1007  partial2R_term[i] -= this->rhs->dot(this->get_adjoint_solution(i));
1008  }
1009 
1010  oldparameters.value_copy(parameters);
1011  parameterperturbation *= -1.0;
1012  parameters += parameterperturbation;
1013 
1014  // Re-center old_parameter, which may be affected by vector
1015  old_parameter = *parameters[k];
1016 
1017  *parameters[k] = old_parameter + delta_p;
1018  this->assemble_qoi(qoi_indices);
1019  this->assembly(true, false, true);
1020  this->rhs->close();
1021  for (unsigned int i=0; i != Nq; ++i)
1022  if (qoi_indices.has_index(i))
1023  {
1024  partial2q_term[i] -= this->qoi[i];
1025  partial2R_term[i] -= this->rhs->dot(this->get_adjoint_solution(i));
1026  }
1027 
1028  *parameters[k] = old_parameter - delta_p;
1029  this->assemble_qoi(qoi_indices);
1030  this->assembly(true, false, true);
1031  this->rhs->close();
1032  for (unsigned int i=0; i != Nq; ++i)
1033  if (qoi_indices.has_index(i))
1034  {
1035  partial2q_term[i] += this->qoi[i];
1036  partial2R_term[i] += this->rhs->dot(this->get_adjoint_solution(i));
1037  }
1038 
1039  for (unsigned int i=0; i != Nq; ++i)
1040  if (qoi_indices.has_index(i))
1041  {
1042  partial2q_term[i] /= (4. * delta_p * delta_p);
1043  partial2R_term[i] /= (4. * delta_p * delta_p);
1044  }
1045 
1046  for (unsigned int i=0; i != Nq; ++i)
1047  if (qoi_indices.has_index(i))
1048  sensitivities[i][k] = partial2q_term[i] - partial2R_term[i];
1049 
1050  // We get (partial q / partial u), R, and
1051  // (partial R / partial u) from the user, but centrally
1052  // difference to get q_uk, R_k, and R_uk terms:
1053  // (partial R / partial k)
1054  // 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)
1055  // (partial^2 q / partial u partial k)
1056  // q_uk = (q_u(p+dp*e_k) - q_u(p-dp*e_k))/(2*dp)
1057  // (partial^2 R / partial u partial k)
1058  // 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)
1059 
1060  // To avoid creating Nq temporary vectors for q_uk or R_uk, we add
1061  // subterms to the sensitivities output one by one.
1062  //
1063  // FIXME: this is probably a bad order of operations for
1064  // controlling floating point error.
1065 
1066  *parameters[k] = old_parameter + delta_p;
1067  this->assembly(true, true);
1068  this->rhs->close();
1069  this->matrix->close();
1070  this->assemble_qoi_derivative(qoi_indices,
1071  /* include_liftfunc = */ true,
1072  /* apply_constraints = */ false);
1073 
1074  this->matrix->vector_mult(*tempvec, this->get_weighted_sensitivity_solution());
1075 
1076  for (unsigned int i=0; i != Nq; ++i)
1077  if (qoi_indices.has_index(i))
1078  {
1079  this->get_adjoint_rhs(i).close();
1080  sensitivities[i][k] += (this->get_adjoint_rhs(i).dot(this->get_weighted_sensitivity_solution()) -
1082  this->get_adjoint_solution(i).dot(*tempvec)) / (2.*delta_p);
1083  }
1084 
1085  *parameters[k] = old_parameter - delta_p;
1086  this->assembly(true, true);
1087  this->rhs->close();
1088  this->matrix->close();
1089  this->assemble_qoi_derivative(qoi_indices,
1090  /* include_liftfunc = */ true,
1091  /* apply_constraints = */ false);
1092 
1093  this->matrix->vector_mult(*tempvec, this->get_weighted_sensitivity_solution());
1094 
1095  for (unsigned int i=0; i != Nq; ++i)
1096  if (qoi_indices.has_index(i))
1097  {
1098  this->get_adjoint_rhs(i).close();
1099  sensitivities[i][k] += (-this->get_adjoint_rhs(i).dot(this->get_weighted_sensitivity_solution()) +
1101  this->get_adjoint_solution(i).dot(*tempvec)) / (2.*delta_p);
1102  }
1103  }
1104 
1105  // All parameters have been reset.
1106  // Don't leave the qoi or system changed - principle of least
1107  // surprise.
1108  this->assembly(true, true);
1109  this->rhs->close();
1110  this->matrix->close();
1111  this->assemble_qoi(qoi_indices);
1112 }
1113 
1114 
1115 
1117  const ParameterVector & parameters_in,
1118  SensitivityData & sensitivities)
1119 {
1120  // We currently get partial derivatives via finite differencing
1121  const Real delta_p = TOLERANCE;
1122 
1123  ParameterVector & parameters =
1124  const_cast<ParameterVector &>(parameters_in);
1125 
1126  // We'll use one temporary vector for matrix-vector-vector products
1127  UniquePtr<NumericVector<Number>> tempvec = this->solution->zero_clone();
1128 
1129  // And another temporary vector to hold a copy of the true solution
1130  // so we can safely perturb this->solution.
1131  UniquePtr<NumericVector<Number>> oldsolution = this->solution->clone();
1132 
1133  const unsigned int Np = cast_int<unsigned int>
1134  (parameters.size());
1135  const unsigned int Nq = cast_int<unsigned int>
1136  (qoi.size());
1137 
1138  // For each quantity of interest q, the parameter sensitivity
1139  // Hessian is defined as q''_{kl} = {d^2 q}/{d p_k d p_l}.
1140  //
1141  // We calculate it from values and partial derivatives of the
1142  // quantity of interest function Q, solution u, adjoint solution z,
1143  // and residual R, as:
1144  //
1145  // q''_{kl} =
1146  // Q''_{kl} + Q''_{uk}(u)*u'_l + Q''_{ul}(u) * u'_k +
1147  // Q''_{uu}(u)*u'_k*u'_l -
1148  // R''_{kl}(u,z) -
1149  // R''_{uk}(u,z)*u'_l - R''_{ul}(u,z)*u'_k -
1150  // R''_{uu}(u,z)*u'_k*u'_l
1151  //
1152  // See the adjoints model document for more details.
1153 
1154  // We first do an adjoint solve to get z for each quantity of
1155  // interest
1156  // if we havent already or dont have an initial condition for the adjoint
1157  if (!this->is_adjoint_already_solved())
1158  {
1159  this->adjoint_solve(qoi_indices);
1160  }
1161 
1162  // And a sensitivity solve to get u_k for each parameter
1163  this->sensitivity_solve(parameters);
1164 
1165  // Get ready to fill in second derivatives:
1166  sensitivities.allocate_hessian_data(qoi_indices, *this, parameters);
1167 
1168  for (unsigned int k=0; k != Np; ++k)
1169  {
1170  Number old_parameterk = *parameters[k];
1171 
1172  // The Hessian is symmetric, so we just calculate the lower
1173  // triangle and the diagonal, and we get the upper triangle from
1174  // the transpose of the lower
1175 
1176  for (unsigned int l=0; l != k+1; ++l)
1177  {
1178  // The second partial derivatives with respect to parameters
1179  // are all calculated via a central finite difference
1180  // stencil:
1181  // F''_{kl} ~= (F(p+dp*e_k+dp*e_l) - F(p+dp*e_k-dp*e_l) -
1182  // F(p-dp*e_k+dp*e_l) + F(p-dp*e_k-dp*e_l))/(4*dp^2)
1183  // We will add Q''_{kl}(u) and subtract R''_{kl}(u,z) at the
1184  // same time.
1185  //
1186  // We have to be careful with the perturbations to handle
1187  // the k=l case
1188 
1189  Number old_parameterl = *parameters[l];
1190 
1191  *parameters[k] += delta_p;
1192  *parameters[l] += delta_p;
1193  this->assemble_qoi(qoi_indices);
1194  this->assembly(true, false, true);
1195  this->rhs->close();
1196  std::vector<Number> partial2q_term = this->qoi;
1197  std::vector<Number> partial2R_term(this->qoi.size());
1198  for (unsigned int i=0; i != Nq; ++i)
1199  if (qoi_indices.has_index(i))
1200  partial2R_term[i] = this->rhs->dot(this->get_adjoint_solution(i));
1201 
1202  *parameters[l] -= 2.*delta_p;
1203  this->assemble_qoi(qoi_indices);
1204  this->assembly(true, false, true);
1205  this->rhs->close();
1206  for (unsigned int i=0; i != Nq; ++i)
1207  if (qoi_indices.has_index(i))
1208  {
1209  partial2q_term[i] -= this->qoi[i];
1210  partial2R_term[i] -= this->rhs->dot(this->get_adjoint_solution(i));
1211  }
1212 
1213  *parameters[k] -= 2.*delta_p;
1214  this->assemble_qoi(qoi_indices);
1215  this->assembly(true, false, true);
1216  this->rhs->close();
1217  for (unsigned int i=0; i != Nq; ++i)
1218  if (qoi_indices.has_index(i))
1219  {
1220  partial2q_term[i] += this->qoi[i];
1221  partial2R_term[i] += this->rhs->dot(this->get_adjoint_solution(i));
1222  }
1223 
1224  *parameters[l] += 2.*delta_p;
1225  this->assemble_qoi(qoi_indices);
1226  this->assembly(true, false, true);
1227  this->rhs->close();
1228  for (unsigned int i=0; i != Nq; ++i)
1229  if (qoi_indices.has_index(i))
1230  {
1231  partial2q_term[i] -= this->qoi[i];
1232  partial2R_term[i] -= this->rhs->dot(this->get_adjoint_solution(i));
1233  partial2q_term[i] /= (4. * delta_p * delta_p);
1234  partial2R_term[i] /= (4. * delta_p * delta_p);
1235  }
1236 
1237  for (unsigned int i=0; i != Nq; ++i)
1238  if (qoi_indices.has_index(i))
1239  {
1240  Number current_terms = partial2q_term[i] - partial2R_term[i];
1241  sensitivities.second_derivative(i,k,l) += current_terms;
1242  if (k != l)
1243  sensitivities.second_derivative(i,l,k) += current_terms;
1244  }
1245 
1246  // Don't leave the parameters perturbed
1247  *parameters[l] = old_parameterl;
1248  *parameters[k] = old_parameterk;
1249  }
1250 
1251  // We get (partial q / partial u) and
1252  // (partial R / partial u) from the user, but centrally
1253  // difference to get q_uk and R_uk terms:
1254  // (partial^2 q / partial u partial k)
1255  // q_uk*u'_l = (q_u(p+dp*e_k)*u'_l - q_u(p-dp*e_k)*u'_l)/(2*dp)
1256  // 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)
1257  //
1258  // To avoid creating Nq temporary vectors, we add these
1259  // subterms to the sensitivities output one by one.
1260  //
1261  // FIXME: this is probably a bad order of operations for
1262  // controlling floating point error.
1263 
1264  *parameters[k] = old_parameterk + delta_p;
1265  this->assembly(false, true);
1266  this->matrix->close();
1267  this->assemble_qoi_derivative(qoi_indices,
1268  /* include_liftfunc = */ true,
1269  /* apply_constraints = */ false);
1270 
1271  for (unsigned int l=0; l != Np; ++l)
1272  {
1273  this->matrix->vector_mult(*tempvec, this->get_sensitivity_solution(l));
1274  for (unsigned int i=0; i != Nq; ++i)
1275  if (qoi_indices.has_index(i))
1276  {
1277  this->get_adjoint_rhs(i).close();
1278  Number current_terms =
1279  (this->get_adjoint_rhs(i).dot(this->get_sensitivity_solution(l)) -
1280  tempvec->dot(this->get_adjoint_solution(i))) / (2.*delta_p);
1281  sensitivities.second_derivative(i,k,l) += current_terms;
1282 
1283  // We use the _uk terms twice; symmetry lets us reuse
1284  // these calculations for the _ul terms.
1285 
1286  sensitivities.second_derivative(i,l,k) += current_terms;
1287  }
1288  }
1289 
1290  *parameters[k] = old_parameterk - delta_p;
1291  this->assembly(false, true);
1292  this->matrix->close();
1293  this->assemble_qoi_derivative(qoi_indices,
1294  /* include_liftfunc = */ true,
1295  /* apply_constraints = */ false);
1296 
1297  for (unsigned int l=0; l != Np; ++l)
1298  {
1299  this->matrix->vector_mult(*tempvec, this->get_sensitivity_solution(l));
1300  for (unsigned int i=0; i != Nq; ++i)
1301  if (qoi_indices.has_index(i))
1302  {
1303  this->get_adjoint_rhs(i).close();
1304  Number current_terms =
1305  (-this->get_adjoint_rhs(i).dot(this->get_sensitivity_solution(l)) +
1306  tempvec->dot(this->get_adjoint_solution(i))) / (2.*delta_p);
1307  sensitivities.second_derivative(i,k,l) += current_terms;
1308 
1309  // We use the _uk terms twice; symmetry lets us reuse
1310  // these calculations for the _ul terms.
1311 
1312  sensitivities.second_derivative(i,l,k) += current_terms;
1313  }
1314  }
1315 
1316  // Don't leave the parameter perturbed
1317  *parameters[k] = old_parameterk;
1318 
1319  // Our last remaining terms are -R_uu(u,z)*u_k*u_l and
1320  // Q_uu(u)*u_k*u_l
1321  //
1322  // We take directional central finite differences of R_u and Q_u
1323  // to approximate these terms, e.g.:
1324  //
1325  // Q_uu(u)*u_k ~= (Q_u(u+dp*u_k) - Q_u(u-dp*u_k))/(2*dp)
1326 
1327  *this->solution = this->get_sensitivity_solution(k);
1328  *this->solution *= delta_p;
1329  *this->solution += *oldsolution;
1330  this->assembly(false, true);
1331  this->matrix->close();
1332  this->assemble_qoi_derivative(qoi_indices,
1333  /* include_liftfunc = */ true,
1334  /* apply_constraints = */ false);
1335 
1336  // The Hessian is symmetric, so we just calculate the lower
1337  // triangle and the diagonal, and we get the upper triangle from
1338  // the transpose of the lower
1339  //
1340  // Note that, because we took the directional finite difference
1341  // with respect to k and not l, we've added an O(delta_p^2)
1342  // error to any permutational symmetry in the Hessian...
1343  for (unsigned int l=0; l != k+1; ++l)
1344  {
1345  this->matrix->vector_mult(*tempvec, this->get_sensitivity_solution(l));
1346  for (unsigned int i=0; i != Nq; ++i)
1347  if (qoi_indices.has_index(i))
1348  {
1349  this->get_adjoint_rhs(i).close();
1350  Number current_terms =
1351  (this->get_adjoint_rhs(i).dot(this->get_sensitivity_solution(l)) -
1352  tempvec->dot(this->get_adjoint_solution(i))) / (2.*delta_p);
1353  sensitivities.second_derivative(i,k,l) += current_terms;
1354  if (k != l)
1355  sensitivities.second_derivative(i,l,k) += current_terms;
1356  }
1357  }
1358 
1359  *this->solution = this->get_sensitivity_solution(k);
1360  *this->solution *= -delta_p;
1361  *this->solution += *oldsolution;
1362  this->assembly(false, true);
1363  this->matrix->close();
1364  this->assemble_qoi_derivative(qoi_indices,
1365  /* include_liftfunc = */ true,
1366  /* apply_constraints = */ false);
1367 
1368  for (unsigned int l=0; l != k+1; ++l)
1369  {
1370  this->matrix->vector_mult(*tempvec, this->get_sensitivity_solution(l));
1371  for (unsigned int i=0; i != Nq; ++i)
1372  if (qoi_indices.has_index(i))
1373  {
1374  this->get_adjoint_rhs(i).close();
1375  Number current_terms =
1376  (-this->get_adjoint_rhs(i).dot(this->get_sensitivity_solution(l)) +
1377  tempvec->dot(this->get_adjoint_solution(i))) / (2.*delta_p);
1378  sensitivities.second_derivative(i,k,l) += current_terms;
1379  if (k != l)
1380  sensitivities.second_derivative(i,l,k) += current_terms;
1381  }
1382  }
1383 
1384  // Don't leave the solution perturbed
1385  *this->solution = *oldsolution;
1386  }
1387 
1388  // All parameters have been reset.
1389  // Don't leave the qoi or system changed - principle of least
1390  // surprise.
1391  this->assembly(true, true);
1392  this->rhs->close();
1393  this->matrix->close();
1394  this->assemble_qoi(qoi_indices);
1395 }
1396 
1397 
1398 
1400 {
1401  // This function allocates memory and hands it back to the user as a
1402  // naked pointer. This makes it too easy to leak memory, and
1403  // therefore this function is deprecated. After a period of
1404  // deprecation, this function will eventually be marked with a
1405  // libmesh_error_msg().
1406  libmesh_deprecated();
1407  // libmesh_error_msg("This function should be overridden by derived classes. "
1408  // "It does not contain a valid LinearSolver to hand back to "
1409  // "the user, so it creates one, opening up the possibility "
1410  // "of a memory leak.");
1411 
1412  LinearSolver<Number> * new_solver =
1413  LinearSolver<Number>::build(this->comm()).release();
1414 
1415  if (libMesh::on_command_line("--solver_system_names"))
1416  new_solver->init((this->name()+"_").c_str());
1417  else
1418  new_solver->init();
1419 
1420  return new_solver;
1421 }
1422 
1423 
1424 
1425 std::pair<unsigned int, Real> ImplicitSystem::get_linear_solve_parameters() const
1426 {
1427  return std::make_pair(this->get_equation_systems().parameters.get<unsigned int>("linear solver maximum iterations"),
1428  this->get_equation_systems().parameters.get<Real>("linear solver tolerance"));
1429 }
1430 
1431 
1432 
1434 {
1435  // This is the counterpart of the get_linear_solver() function, which is now deprecated.
1436  libmesh_deprecated();
1437 
1438  delete s;
1439 }
1440 
1441 } // namespace libMesh
virtual LinearSolver< Number > * get_linear_solver() const
double abs(double a)
This is the EquationSystems class.
virtual std::pair< unsigned int, Real > adjoint_solve(const QoISet &qoi_indices=QoISet()) libmesh_override
Assembles & solves the linear system (dR/du)^T*z = dq/du, for those quantities of interest q specifie...
virtual void get_transpose(SparseMatrix< T > &dest) const =0
Copies the transpose of the matrix into dest, which may be *this.
virtual void clear() libmesh_override
Clear all the data structures associated with the system.
NumericVector< Number > & add_adjoint_solution(unsigned int i=0)
Definition: system.C:977
Data structure for specifying which Parameters should be independent variables in a parameter sensiti...
virtual void disable_cache() libmesh_override
Avoids use of any cached data that might affect any solve result.
static UniquePtr< LinearSolver< T > > build(const libMesh::Parallel::Communicator &comm_in, const SolverPackage solver_package=libMesh::default_solver_package())
Builds a LinearSolver using the linear solver package specified by solver_package.
Definition: linear_solver.C:42
Data structure for specifying which Quantities of Interest should be calculated in an adjoint or a pa...
Definition: qoi_set.h:45
void value_copy(ParameterVector &target) const
Value copy method: the target, which should already have as many parameters as I do, will now have those parameters set to my values.
void allocate_hessian_data(const QoISet &qoi_indices, const System &sys, const ParameterVector &parameter_vector)
Given QoISet and ParameterVector, allocates space for all required second derivative data...
NumericVector< Number > & get_sensitivity_solution(unsigned int i=0)
Definition: system.C:936
virtual void adjoint_qoi_parameter_sensitivity(const QoISet &qoi_indices, const ParameterVector &parameters, SensitivityData &sensitivities) libmesh_override
Solves for the derivative of each of the system&#39;s quantities of interest q in qoi[qoi_indices] with r...
NumericVector< Number > & get_sensitivity_rhs(unsigned int i=0)
Definition: system.C:1081
bool is_attached(SparseMatrix< Number > &matrix)
Matrices should not be attached more than once.
Definition: dof_map.C:282
virtual void assemble()
Prepares matrix and _dof_map for matrix assembly.
Definition: system.C:479
void attach_matrix(SparseMatrix< Number > &matrix)
Additional matrices may be handled with this DofMap.
Definition: dof_map.C:246
virtual void init_data()
Initializes the data for the system.
Definition: system.C:260
virtual void qoi_parameter_hessian_vector_product(const QoISet &qoi_indices, const ParameterVector &parameters, const ParameterVector &vector, SensitivityData &product) libmesh_override
For each of the system&#39;s quantities of interest q in qoi[qoi_indices], and for a vector of parameters...
SparseMatrix< Number > & add_matrix(const std::string &mat_name)
Adds the additional matrix mat_name to this system.
void vector_mult(NumericVector< T > &dest, const NumericVector< T > &arg) const
Multiplies the matrix by the NumericVector arg and stores the result in NumericVector dest...
const class libmesh_nullptr_t libmesh_nullptr
virtual std::pair< unsigned int, Real > sensitivity_solve(const ParameterVector &parameters) libmesh_override
Assembles & solves the linear system(s) (dR/du)*u_p = -dR/dp, for those parameters contained within p...
NumericVector< Number > * rhs
The system matrix.
virtual ~ImplicitSystem()
Destructor.
void vector_mult_add(NumericVector< T > &dest, const NumericVector< T > &arg) const
Multiplies the matrix by the NumericVector arg and adds the result to the NumericVector dest...
const T & get(const std::string &) const
Definition: parameters.h:430
NumericVector< Number > & add_sensitivity_rhs(unsigned int i=0)
Definition: system.C:1071
static const Real TOLERANCE
static UniquePtr< SparseMatrix< T > > build(const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package())
Builds a SparseMatrix<T> using the linear solver package specified by solver_package.
std::map< std::string, SparseMatrix< Number > * > _matrices
Some systems need an arbitrary number of matrices.
virtual bool initialized() const
const Number & second_derivative(unsigned int qoi_index, unsigned int parameter_index1, unsigned int parameter_index2) const
The libMesh namespace provides an interface to certain functionality in the library.
PetscErrorCode Vec Mat Mat pc
bool zero_out_matrix_and_rhs
By default, the system will zero out the matrix and the right hand side.
NumericVector< Number > & add_weighted_sensitivity_solution()
Definition: system.C:956
long double max(long double a, double b)
const std::string & name() const
Definition: system.h:1998
NumericVector< Number > & get_weighted_sensitivity_solution()
Definition: system.C:963
virtual std::pair< unsigned int, Real > get_linear_solve_parameters() const
virtual void zero()=0
Set all entries to zero.
libmesh_assert(j)
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
virtual void assemble_qoi_derivative(const QoISet &qoi_indices=QoISet(), bool include_liftfunc=true, bool apply_constraints=true) libmesh_override
Prepares adjoint_rhs for quantity of interest derivative assembly, then calls user qoi derivative fun...
bool has_adjoint_dirichlet_boundaries(unsigned int q) const
NumericVector< Number > & add_weighted_sensitivity_adjoint_solution(unsigned int i=0)
Definition: system.C:1009
virtual bool initialized() const
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:167
void allocate_data(const QoISet &qoi_indices, const System &sys, const ParameterVector &parameter_vector)
Given QoISet and ParameterVector, allocates space for all required first derivative data...
const MeshBase & get_mesh() const
Definition: system.h:2014
virtual void assembly(bool, bool, bool=false, bool=false)
Assembles a residual in rhs and/or a jacobian in matrix, as requested.
const DofMap & get_dof_map() const
Definition: system.h:2030
virtual void zero()=0
Set all entries to 0.
Data structure for holding completed parameter sensitivity calculations.
std::vector< Number > qoi
Values of the quantities of interest.
Definition: system.h:1553
virtual void init(const char *name=libmesh_nullptr)=0
Initialize data structures if not done so already.
virtual void qoi_parameter_hessian(const QoISet &qoi_indices, const ParameterVector &parameters, SensitivityData &hessian) libmesh_override
For each of the system&#39;s quantities of interest q in qoi[qoi_indices], and for a vector of parameters...
bool have_matrix(const std::string &mat_name) const
This is the base class for classes which contain information related to any physical process that mig...
Definition: system.h:76
virtual std::pair< unsigned int, Real > adjoint_solve(SparseMatrix< T > &, NumericVector< T > &, NumericVector< T > &, const double, const unsigned int)
Function to solve the adjoint system.
virtual std::pair< unsigned int, Real > solve(SparseMatrix< T > &, NumericVector< T > &, NumericVector< T > &, const double, const unsigned int)=0
This function calls the solver _solver_type preconditioned with the _preconditioner_type precondition...
UniquePtr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1523
void enforce_constraints_exactly(const System &system, NumericVector< Number > *v=libmesh_nullptr, bool homogeneous=false) const
Constrains the numeric vector v, which represents a solution defined on the mesh. ...
Definition: dof_map.h:1816
virtual void reinit() libmesh_override
Reinitializes the member data fields associated with the system, so that, e.g., assemble() may be use...
bool is_adjoint_already_solved() const
Accessor for the adjoint_already_solved boolean.
Definition: system.h:372
virtual void close()=0
Calls the NumericVector&#39;s internal assembly routines, ensuring that the values are consistent across ...
virtual void reinit() libmesh_override
Reinitializes the member data fields associated with the system, so that, e.g., assemble() may be use...
const SparseMatrix< Number > & get_matrix(const std::string &mat_name) const
std::size_t size() const
void add_system_matrix()
Add the system matrix to the _matrices data structure.
const SparseMatrix< Number > * request_matrix(const std::string &mat_name) const
ImplicitSystem(EquationSystems &es, const std::string &name, const unsigned int number)
Constructor.
virtual void close()=0
Calls the SparseMatrix&#39;s internal assembly routines, ensuring that the values are consistent across p...
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual void init_data() libmesh_override
Initializes the member data fields associated with the system, so that, e.g., assemble() may be used...
std::map< std::string, SparseMatrix< Number > * >::const_iterator const_matrices_iterator
virtual void reuse_preconditioner(bool)
Set the same_preconditioner flag, which indicates if we reuse the same preconditioner for subsequent ...
void deep_copy(ParameterVector &target) const
Deep copy constructor: the target will now own new copies of all the parameter values I&#39;m pointing to...
bool _can_add_matrices
true when additional matrices may still be added, false otherwise.
bool on_command_line(const std::string &arg)
Definition: libmesh.C:921
virtual UniquePtr< NumericVector< T > > zero_clone() const =0
void remove_matrix(const std::string &mat_name)
Removes the additional matrix mat_name from this system.
virtual std::pair< unsigned int, Real > weighted_sensitivity_adjoint_solve(const ParameterVector &parameters, const ParameterVector &weights, const QoISet &qoi_indices=QoISet()) libmesh_override
Assembles & solves the linear system(s) (dR/du)^T*z_w = sum(w_p*(d^2q/dudp - d^2R/dudp*z)), for those parameters p contained within parameters, weighted by the values w_p found within weights.
const Parallel::Communicator & comm() const
SparseMatrix< Number > * matrix
The system matrix.
std::map< std::string, SparseMatrix< Number > * >::iterator matrices_iterator
Matrix iterator typedefs.
virtual void release_linear_solver(LinearSolver< Number > *) const
Releases a pointer to a linear solver acquired by this->get_linear_solver()
void clear_sparsity()
Clears the sparsity pattern.
Definition: dof_map.C:1772
NumericVector< Number > & get_adjoint_solution(unsigned int i=0)
Definition: system.C:989
const EquationSystems & get_equation_systems() const
Definition: system.h:712
Parameters parameters
Data structure holding arbitrary parameters.
virtual std::pair< unsigned int, Real > weighted_sensitivity_solve(const ParameterVector &parameters, const ParameterVector &weights) libmesh_override
Assembles & solves the linear system(s) (dR/du)*u_w = sum(w_p*-dR/dp), for those parameters p contain...
virtual T dot(const NumericVector< T > &v) const =0
virtual void forward_qoi_parameter_sensitivity(const QoISet &qoi_indices, const ParameterVector &parameters, SensitivityData &sensitivities) libmesh_override
Solves for the derivative of each of the system&#39;s quantities of interest q in qoi[qoi_indices] with r...
virtual void assemble_qoi(const QoISet &qoi_indices=QoISet()) libmesh_override
Prepares qoi for quantity of interest assembly, then calls user qoi function.
virtual void assemble() libmesh_override
Prepares matrix and rhs for system assembly, then calls user assembly function.
virtual void assemble_residual_derivatives(const ParameterVector &parameters) libmesh_override
Residual parameter derivative function.
NumericVector< Number > & add_sensitivity_solution(unsigned int i=0)
Definition: system.C:926
void compute_sparsity(const MeshBase &)
Computes the sparsity pattern for the matrices corresponding to proc_id and sends that data to Linear...
Definition: dof_map.C:1735
bool assemble_before_solve
Flag which tells the system to whether or not to call the user assembly function during each call to ...
Definition: system.h:1477
virtual void init_matrices()
Initializes the matrices associated with this system.
bool has_index(unsigned int) const
Return whether or not this index is in the set to be calculated.
Definition: qoi_set.h:221
void enforce_adjoint_constraints_exactly(NumericVector< Number > &v, unsigned int q) const
Heterogenously constrains the numeric vector v, which represents an adjoint solution defined on the m...
Definition: dof_map.h:1820
virtual void clear() libmesh_override
Clear all the data structures associated with the system.
NumericVector< Number > & get_weighted_sensitivity_adjoint_solution(unsigned int i=0)
Definition: system.C:1021
virtual UniquePtr< NumericVector< T > > clone() const =0
NumericVector< Number > & get_adjoint_rhs(unsigned int i=0)
Definition: system.C:1051