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