implicit_system.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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"
27 #include "libmesh/linear_solver.h"
28 #include "libmesh/mesh_base.h"
29 #include "libmesh/numeric_vector.h"
30 #include "libmesh/parameters.h"
32 #include "libmesh/qoi_set.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 (NULL),
47  _can_add_matrices (true)
48 {
49 }
50 
51 
52 
54 {
55  // Clear data
56  this->clear();
57 }
58 
59 
60 
62 {
63  // clear the parent data
64  Parent::clear();
65 
66  // clear any user-added matrices
67  {
68  for (matrices_iterator pos = _matrices.begin();
69  pos != _matrices.end(); ++pos)
70  {
71  pos->second->clear ();
72  delete pos->second;
73  pos->second = NULL;
74  }
75 
76  _matrices.clear();
77  _can_add_matrices = true;
78  }
79 
80  // NULL the matrix.
81  matrix = NULL;
82 }
83 
84 
85 
87 {
88  // initialize parent data
90 
91  // Clear any existing matrices
92  for (matrices_iterator pos = _matrices.begin();
93  pos != _matrices.end(); ++pos)
94  pos->second->clear();
95 
96  // Add the system matrix.
97  this->add_system_matrix ();
98 
99  // Initialize the matrices for the system
100  this->init_matrices ();
101 }
102 
103 
104 
106 {
108 
109  // Check for quick return in case the system matrix
110  // (and by extension all the matrices) has already
111  // been initialized
112  if (matrix->initialized())
113  return;
114 
115  // Get a reference to the DofMap
116  DofMap& dof_map = this->get_dof_map();
117 
118  // no chance to add other matrices
119  _can_add_matrices = false;
120 
121  // Tell the matrices about the dof map, and vice versa
122  for (matrices_iterator pos = _matrices.begin();
123  pos != _matrices.end(); ++pos)
124  {
125  SparseMatrix<Number> &m = *(pos->second);
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 (matrices_iterator pos = _matrices.begin();
141  pos != _matrices.end(); ++pos)
142  pos->second->init ();
143 
144  // Set the additional matrices to 0.
145  for (matrices_iterator pos = _matrices.begin();
146  pos != _matrices.end(); ++pos)
147  pos->second->zero ();
148 }
149 
150 
151 
153 {
154  // initialize parent data
155  Parent::reinit();
156 
157  // Get a reference to the DofMap
158  DofMap& dof_map = this->get_dof_map();
159 
160  // Clear the matrices
161  for (matrices_iterator pos = _matrices.begin();
162  pos != _matrices.end(); ++pos)
163  {
164  pos->second->clear();
165  pos->second->attach_dof_map (dof_map);
166  }
167 
168  // Clear the sparsity pattern
169  this->get_dof_map().clear_sparsity();
170 
171  // Compute the sparsity pattern for the current
172  // mesh and DOF distribution. This also updates
173  // additional matrices, \p DofMap now knows them
174  dof_map.compute_sparsity (this->get_mesh());
175 
176  // Initialize matrices
177  for (matrices_iterator pos = _matrices.begin();
178  pos != _matrices.end(); ++pos)
179  pos->second->init ();
180 
181  // Set the additional matrices to 0.
182  for (matrices_iterator pos = _matrices.begin();
183  pos != _matrices.end(); ++pos)
184  pos->second->zero ();
185 }
186 
187 
188 
190 {
195 
196  // The user assembly gets to expect to accumulate on an initially
197  // empty system
198  matrix->zero ();
199  rhs->zero ();
200 
201  // Call the base class assemble function
202  Parent::assemble ();
203 }
204 
205 
206 
207 SparseMatrix<Number> & ImplicitSystem::add_matrix (const std::string& mat_name)
208 {
209  // only add matrices before initializing...
210  if (!_can_add_matrices)
211  {
212  libMesh::err << "ERROR: Too late. Cannot add matrices to the system after initialization"
213  << std::endl
214  << " any more. You should have done this earlier."
215  << std::endl;
216  libmesh_error();
217  }
218 
219  // Return the matrix if it is already there.
220  if (this->have_matrix(mat_name))
221  return *(_matrices[mat_name]);
222 
223  // Otherwise build the matrix and return it.
224  SparseMatrix<Number>* buf = SparseMatrix<Number>::build(this->comm()).release();
225  _matrices.insert (std::make_pair (mat_name, buf));
226 
227  return *buf;
228 }
229 
230 
231 
232 const SparseMatrix<Number> * ImplicitSystem::request_matrix (const std::string& mat_name) const
233 {
234  // Make sure the matrix exists
235  const_matrices_iterator pos = _matrices.find (mat_name);
236 
237  if (pos == _matrices.end())
238  return NULL;
239 
240  return pos->second;
241 }
242 
243 
244 
246 {
247  // Make sure the matrix exists
248  matrices_iterator pos = _matrices.find (mat_name);
249 
250  if (pos == _matrices.end())
251  return NULL;
252 
253  return pos->second;
254 }
255 
256 
257 
258 const SparseMatrix<Number> & ImplicitSystem::get_matrix (const std::string& mat_name) const
259 {
260  // Make sure the matrix exists
261  const_matrices_iterator pos = _matrices.find (mat_name);
262 
263  if (pos == _matrices.end())
264  {
265  libMesh::err << "ERROR: matrix "
266  << mat_name
267  << " does not exist in this system!"
268  << std::endl;
269  libmesh_error();
270  }
271 
272  return *(pos->second);
273 }
274 
275 
276 
277 SparseMatrix<Number> & ImplicitSystem::get_matrix (const std::string& mat_name)
278 {
279  // Make sure the matrix exists
280  matrices_iterator pos = _matrices.find (mat_name);
281 
282  if (pos == _matrices.end())
283  {
284  libMesh::err << "ERROR: matrix "
285  << mat_name
286  << " does not exist in this system!"
287  << std::endl;
288  libmesh_error();
289  }
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 = NULL;
301 
302 
303  // Only need to add the matrix if it isn't there
304  // already!
305  if (matrix == NULL)
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  START_LOG("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 (unsigned int 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 (unsigned int 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  // Stop logging the nonlinear solve
374  STOP_LOG("sensitivity_solve()", "ImplicitSystem");
375 
376  return totalrval;
377 }
378 
379 
380 
381 std::pair<unsigned int, Real>
383 {
384  // Log how long the linear solve takes.
385  START_LOG("adjoint_solve()", "ImplicitSystem");
386 
387  if (this->assemble_before_solve)
388  // Assemble the linear system
389  this->assembly (/* get_residual = */ false,
390  /* get_jacobian = */ true);
391 
392  // The adjoint problem is linear
393  LinearSolver<Number> *linear_solver = this->get_linear_solver();
394 
395  // Reset and build the RHS from the QOI derivative
396  this->assemble_qoi_derivative(qoi_indices);
397 
398  // Our iteration counts and residuals will be sums of the individual
399  // results
400  std::pair<unsigned int, Real> solver_params =
402  std::pair<unsigned int, Real> totalrval = std::make_pair(0,0.0);
403 
404  for (unsigned int i=0; i != this->qoi.size(); ++i)
405  if (qoi_indices.has_index(i))
406  {
407  const std::pair<unsigned int, Real> rval =
408  linear_solver->adjoint_solve (*matrix, this->add_adjoint_solution(i),
409  this->get_adjoint_rhs(i),
410  solver_params.second,
411  solver_params.first);
412 
413  totalrval.first += rval.first;
414  totalrval.second += rval.second;
415  }
416 
417  this->release_linear_solver(linear_solver);
418 
419  // The linear solver may not have fit our constraints exactly
420 #ifdef LIBMESH_ENABLE_CONSTRAINTS
421  for (unsigned int i=0; i != this->qoi.size(); ++i)
422  if (qoi_indices.has_index(i))
424  (this->get_adjoint_solution(i), i);
425 #endif
426 
427  // Stop logging the nonlinear solve
428  STOP_LOG("adjoint_solve()", "ImplicitSystem");
429 
430  return totalrval;
431 }
432 
433 
434 
435 std::pair<unsigned int, Real>
437  const ParameterVector& weights,
438  const QoISet& qoi_indices)
439 {
440  // Log how long the linear solve takes.
441  START_LOG("weighted_sensitivity_adjoint_solve()", "ImplicitSystem");
442 
443  // We currently get partial derivatives via central differencing
444  const Real delta_p = TOLERANCE;
445 
446  // The forward system should now already be solved.
447  // The adjoint system should now already be solved.
448  // Now we're assembling a weighted sum of adjoint-adjoint systems:
449  //
450  // dR/du (u, sum_l(w_l*z^l)) = sum_l(w_l*(Q''_ul - R''_ul (u, z)))
451 
452  // FIXME: The derivation here does not yet take adjoint boundary
453  // conditions into account.
454  for (unsigned int i=0; i != this->qoi.size(); ++i)
455  if (qoi_indices.has_index(i))
457 
458  // We'll assemble the rhs first, because the R'' term will require
459  // perturbing the jacobian
460 
461  // We'll use temporary rhs vectors, because we haven't (yet) found
462  // any good reasons why users might want to save these:
463 
464  std::vector<NumericVector<Number> *> temprhs(this->qoi.size());
465  for (unsigned int i=0; i != this->qoi.size(); ++i)
466  if (qoi_indices.has_index(i))
467  temprhs[i] = this->rhs->zero_clone().release();
468 
469  // We approximate the _l partial derivatives via a central
470  // differencing perturbation in the w_l direction:
471  //
472  // sum_l(w_l*v_l) ~= (v(p + dp*w_l*e_l) - v(p - dp*w_l*e_l))/(2*dp)
473 
474  // PETSc doesn't implement SGEMX, so neither does NumericVector,
475  // so we want to avoid calculating f -= R'*z. We'll thus evaluate
476  // the above equation by first adding -v(p+dp...), then multiplying
477  // the intermediate result vectors by -1, then adding -v(p-dp...),
478  // then finally dividing by 2*dp.
479 
480  ParameterVector oldparameters, parameterperturbation;
481  parameters.deep_copy(oldparameters);
482  weights.deep_copy(parameterperturbation);
483  parameterperturbation *= delta_p;
484  parameters += parameterperturbation;
485 
486  this->assembly(false, true);
487  this->matrix->close();
488 
489  // Take the discrete adjoint, so that we can calculate R_u(u,z) with
490  // a matrix-vector product of R_u and z.
492 
493  this->assemble_qoi_derivative(qoi_indices);
494  for (unsigned int i=0; i != this->qoi.size(); ++i)
495  if (qoi_indices.has_index(i))
496  {
497  this->get_adjoint_rhs(i).close();
498  *(temprhs[i]) -= this->get_adjoint_rhs(i);
499  this->matrix->vector_mult_add(*(temprhs[i]), this->get_adjoint_solution(i));
500  *(temprhs[i]) *= -1.0;
501  }
502 
503  oldparameters.value_copy(parameters);
504  parameterperturbation *= -1.0;
505  parameters += parameterperturbation;
506 
507  this->assembly(false, true);
508  this->matrix->close();
510 
511  this->assemble_qoi_derivative(qoi_indices);
512  for (unsigned int i=0; i != this->qoi.size(); ++i)
513  if (qoi_indices.has_index(i))
514  {
515  this->get_adjoint_rhs(i).close();
516  *(temprhs[i]) -= this->get_adjoint_rhs(i);
517  this->matrix->vector_mult_add(*(temprhs[i]), this->get_adjoint_solution(i));
518  *(temprhs[i]) /= (2.0*delta_p);
519  }
520 
521  // Finally, assemble the jacobian at the non-perturbed parameter
522  // values. Ignore assemble_before_solve; if we had a good
523  // non-perturbed matrix before we've already overwritten it.
524  oldparameters.value_copy(parameters);
525 
526  // if (this->assemble_before_solve)
527  {
528  // Build the Jacobian
529  this->assembly(false, true);
530  this->matrix->close();
531 
532  // Take the discrete adjoint
534  }
535 
536  // The weighted adjoint-adjoint problem is linear
537  LinearSolver<Number> *linear_solver = this->get_linear_solver();
538 
539  // Our iteration counts and residuals will be sums of the individual
540  // results
541  std::pair<unsigned int, Real> solver_params =
543  std::pair<unsigned int, Real> totalrval = std::make_pair(0,0.0);
544 
545  for (unsigned int i=0; i != this->qoi.size(); ++i)
546  if (qoi_indices.has_index(i))
547  {
548  const std::pair<unsigned int, Real> rval =
549  linear_solver->solve (*matrix, this->add_weighted_sensitivity_adjoint_solution(i),
550  *(temprhs[i]),
551  solver_params.second,
552  solver_params.first);
553 
554  totalrval.first += rval.first;
555  totalrval.second += rval.second;
556  }
557 
558  this->release_linear_solver(linear_solver);
559 
560  for (unsigned int i=0; i != this->qoi.size(); ++i)
561  if (qoi_indices.has_index(i))
562  delete temprhs[i];
563 
564  // The linear solver may not have fit our constraints exactly
565 #ifdef LIBMESH_ENABLE_CONSTRAINTS
566  for (unsigned int i=0; i != this->qoi.size(); ++i)
567  if (qoi_indices.has_index(i))
570  /* homogeneous = */ true);
571 #endif
572 
573  // Stop logging the nonlinear solve
574  STOP_LOG("weighted_sensitivity_adjoint_solve()", "ImplicitSystem");
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  START_LOG("weighted_sensitivity_solve()", "ImplicitSystem");
587 
588  // We currently get partial derivatives via central differencing
589  const Real delta_p = TOLERANCE;
590 
591  // The forward system should now already be solved.
592 
593  // Now we're assembling a weighted sum of sensitivity systems:
594  //
595  // dR/du (u, v)(sum(w_l*u'_l)) = -sum_l(w_l*R'_l (u, v)) forall v
596 
597  // We'll assemble the rhs first, because the R' term will require
598  // perturbing the system, and some applications may not be able to
599  // assemble a perturbed residual without simultaneously constructing
600  // a perturbed jacobian.
601 
602  // We approximate the _l partial derivatives via a central
603  // differencing perturbation in the w_l direction:
604  //
605  // sum_l(w_l*v_l) ~= (v(p + dp*w_l*e_l) - v(p - dp*w_l*e_l))/(2*dp)
606 
607  ParameterVector oldparameters, parameterperturbation;
608  parameters.deep_copy(oldparameters);
609  weights.deep_copy(parameterperturbation);
610  parameterperturbation *= delta_p;
611  parameters += parameterperturbation;
612 
613  this->assembly(true, false);
614  this->rhs->close();
615 
616  AutoPtr<NumericVector<Number> > temprhs = this->rhs->clone();
617 
618  oldparameters.value_copy(parameters);
619  parameterperturbation *= -1.0;
620  parameters += parameterperturbation;
621 
622  this->assembly(true, false);
623  this->rhs->close();
624 
625  *temprhs -= *(this->rhs);
626  *temprhs /= (2.0*delta_p);
627 
628  // Finally, assemble the jacobian at the non-perturbed parameter
629  // values
630  oldparameters.value_copy(parameters);
631 
632  // Build the Jacobian
633  this->assembly(false, true);
634  this->matrix->close();
635 
636  // The weighted sensitivity problem is linear
637  LinearSolver<Number> *linear_solver = this->get_linear_solver();
638 
639  std::pair<unsigned int, Real> solver_params =
641 
642  const std::pair<unsigned int, Real> rval =
643  linear_solver->solve (*matrix, this->add_weighted_sensitivity_solution(),
644  *temprhs,
645  solver_params.second,
646  solver_params.first);
647 
648  this->release_linear_solver(linear_solver);
649 
650  // The linear solver may not have fit our constraints exactly
651 #ifdef LIBMESH_ENABLE_CONSTRAINTS
653  (*this, &this->get_weighted_sensitivity_solution(),
654  /* homogeneous = */ true);
655 #endif
656 
657  // Stop logging the nonlinear solve
658  STOP_LOG("weighted_sensitivity_solve()", "ImplicitSystem");
659 
660  return rval;
661 }
662 
663 
664 
666 {
667  const unsigned int Np = libmesh_cast_int<unsigned int>
668  (parameters.size());
669  Real deltap = TOLERANCE;
670 
671  for (unsigned int p=0; p != Np; ++p)
672  {
673  NumericVector<Number> &sensitivity_rhs = this->add_sensitivity_rhs(p);
674 
675  // Approximate -(partial R / partial p) by
676  // (R(p-dp) - R(p+dp)) / (2*dp)
677 
678  Number old_parameter = *parameters[p];
679  *parameters[p] -= deltap;
680 
681  this->assembly(true, false);
682  this->rhs->close();
683  sensitivity_rhs = *this->rhs;
684 
685  *parameters[p] = old_parameter + deltap;
686 
687  this->assembly(true, false);
688  this->rhs->close();
689 
690  sensitivity_rhs -= *this->rhs;
691  sensitivity_rhs /= (2*deltap);
692  sensitivity_rhs.close();
693 
694  *parameters[p] = old_parameter;
695  }
696 }
697 
698 
699 
701  (const QoISet& qoi_indices,
702  const ParameterVector& parameters,
703  SensitivityData& sensitivities)
704 {
705  const unsigned int Np = libmesh_cast_int<unsigned int>
706  (parameters.size());
707  const unsigned int Nq = libmesh_cast_int<unsigned int>
708  (qoi.size());
709 
710  // We currently get partial derivatives via central differencing
711  const Real delta_p = TOLERANCE;
712 
713  // An introduction to the problem:
714  //
715  // Residual R(u(p),p) = 0
716  // partial R / partial u = J = system matrix
717  //
718  // This implies that:
719  // d/dp(R) = 0
720  // (partial R / partial p) +
721  // (partial R / partial u) * (partial u / partial p) = 0
722 
723  // We first do an adjoint solve:
724  // J^T * z = (partial q / partial u)
725  // if we havent already or dont have an initial condition for the adjoint
726  if (!this->is_adjoint_already_solved())
727  {
728  this->adjoint_solve(qoi_indices);
729  }
730 
731  // Get ready to fill in senstivities:
732  sensitivities.allocate_data(qoi_indices, *this, parameters);
733 
734  // We use the identities:
735  // dq/dp = (partial q / partial p) + (partial q / partial u) *
736  // (partial u / partial p)
737  // dq/dp = (partial q / partial p) + (J^T * z) *
738  // (partial u / partial p)
739  // dq/dp = (partial q / partial p) + z * J *
740  // (partial u / partial p)
741 
742  // Leading to our final formula:
743  // dq/dp = (partial q / partial p) - z * (partial R / partial p)
744 
745  // In the case of adjoints with heterogenous Dirichlet boundary
746  // function phi, where
747  // q := R(u,phi) + S(u)
748  // the final formula works out to:
749  // dq/dp = (partial S / partial p) - z * (partial R / partial p)
750  // Because we currently have no direct access to
751  // (partial S / partial p), we use the identity
752  // (partial S / partial p) = (partial q / partial p) -
753  // phi * (partial R / partial p)
754  // to derive an equivalent equation:
755  // dq/dp = (partial q / partial p) - (z+phi) * (partial R / partial p)
756 
757  for (unsigned int j=0; j != Np; ++j)
758  {
759  // (partial q / partial p) ~= (q(p+dp)-q(p-dp))/(2*dp)
760  // (partial R / partial p) ~= (rhs(p+dp) - rhs(p-dp))/(2*dp)
761 
762  Number old_parameter = *parameters[j];
763  // Number old_qoi = this->qoi;
764 
765  *parameters[j] = old_parameter - delta_p;
766  this->assemble_qoi(qoi_indices);
767  std::vector<Number> qoi_minus = this->qoi;
768 
769  this->assembly(true, false);
770  this->rhs->close();
771 
772  // FIXME - this can and should be optimized to avoid the clone()
773  AutoPtr<NumericVector<Number> > partialR_partialp = this->rhs->clone();
774  *partialR_partialp *= -1;
775 
776  *parameters[j] = old_parameter + delta_p;
777  this->assemble_qoi(qoi_indices);
778  std::vector<Number>& qoi_plus = this->qoi;
779 
780  std::vector<Number> partialq_partialp(Nq, 0);
781  for (unsigned int i=0; i != Nq; ++i)
782  if (qoi_indices.has_index(i))
783  partialq_partialp[i] = (qoi_plus[i] - qoi_minus[i]) / (2.*delta_p);
784 
785  this->assembly(true, false);
786  this->rhs->close();
787  *partialR_partialp += *this->rhs;
788  *partialR_partialp /= (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  {
796 
797  if (this->get_dof_map().has_adjoint_dirichlet_boundaries(i))
798  {
799  AutoPtr<NumericVector<Number> > lift_func =
800  this->get_adjoint_solution(i).zero_clone();
801  this->get_dof_map().enforce_constraints_exactly
802  (*this, lift_func.get(),
803  /* homogeneous = */ false);
804  sensitivities[i][j] = partialq_partialp[i] -
805  partialR_partialp->dot(*lift_func) -
806  partialR_partialp->dot(this->get_adjoint_solution(i));
807  }
808  else
809  sensitivities[i][j] = partialq_partialp[i] -
810  partialR_partialp->dot(this->get_adjoint_solution(i));
811  }
812  }
813 
814  // All parameters have been reset.
815  // We didn't cache the original rhs or matrix for memory reasons,
816  // but we can restore them to a state consistent solution -
817  // principle of least surprise.
818  this->assembly(true, true);
819  this->rhs->close();
820  this->matrix->close();
821  this->assemble_qoi(qoi_indices);
822 }
823 
824 
825 
827  (const QoISet& qoi_indices,
828  const ParameterVector& parameters,
829  SensitivityData& sensitivities)
830 {
831  const unsigned int Np = libmesh_cast_int<unsigned int>
832  (parameters.size());
833  const unsigned int Nq = libmesh_cast_int<unsigned int>
834  (qoi.size());
835 
836  // We currently get partial derivatives via central differencing
837  const Real delta_p = TOLERANCE;
838 
839  // An introduction to the problem:
840  //
841  // Residual R(u(p),p) = 0
842  // partial R / partial u = J = system matrix
843  //
844  // This implies that:
845  // d/dp(R) = 0
846  // (partial R / partial p) +
847  // (partial R / partial u) * (partial u / partial p) = 0
848 
849  // We first solve for (partial u / partial p) for each parameter:
850  // J * (partial u / partial p) = - (partial R / partial p)
851 
852  this->sensitivity_solve(parameters);
853 
854  // Get ready to fill in senstivities:
855  sensitivities.allocate_data(qoi_indices, *this, parameters);
856 
857  // We use the identity:
858  // dq/dp = (partial q / partial p) + (partial q / partial u) *
859  // (partial u / partial p)
860 
861  // We get (partial q / partial u) from the user
862  this->assemble_qoi_derivative(qoi_indices);
863 
864  // FIXME: what do we do with adjoint boundary conditions here?
865 
866  // We don't need these to be closed() in this function, but libMesh
867  // standard practice is to have them closed() by the time the
868  // function exits
869  for (unsigned int i=0; i != this->qoi.size(); ++i)
870  if (qoi_indices.has_index(i))
871  this->get_adjoint_rhs(i).close();
872 
873  for (unsigned int j=0; j != Np; ++j)
874  {
875  // (partial q / partial p) ~= (q(p+dp)-q(p-dp))/(2*dp)
876 
877  Number old_parameter = *parameters[j];
878 
879  *parameters[j] = old_parameter - delta_p;
880  this->assemble_qoi();
881  std::vector<Number> qoi_minus = this->qoi;
882 
883  *parameters[j] = old_parameter + delta_p;
884  this->assemble_qoi();
885  std::vector<Number>& qoi_plus = this->qoi;
886 
887  std::vector<Number> partialq_partialp(Nq, 0);
888  for (unsigned int i=0; i != Nq; ++i)
889  if (qoi_indices.has_index(i))
890  partialq_partialp[i] = (qoi_plus[i] - qoi_minus[i]) / (2.*delta_p);
891 
892  // Don't leave the parameter changed
893  *parameters[j] = old_parameter;
894 
895  for (unsigned int i=0; i != Nq; ++i)
896  if (qoi_indices.has_index(i))
897  sensitivities[i][j] = partialq_partialp[i] +
898  this->get_adjoint_rhs(i).dot(this->get_sensitivity_solution(j));
899  }
900 
901  // All parameters have been reset.
902  // We didn't cache the original rhs or matrix for memory reasons,
903  // but we can restore them to a state consistent solution -
904  // principle of least surprise.
905  this->assembly(true, true);
906  this->rhs->close();
907  this->matrix->close();
908  this->assemble_qoi(qoi_indices);
909 }
910 
911 
912 
914  (const QoISet& qoi_indices,
915  const ParameterVector& parameters,
916  const ParameterVector& vector,
917  SensitivityData& sensitivities)
918 {
919  // We currently get partial derivatives via finite differencing
920  const Real delta_p = TOLERANCE;
921 
922  // We'll use a single temporary vector for matrix-vector-vector products
923  AutoPtr<NumericVector<Number> > tempvec = this->solution->zero_clone();
924 
925  const unsigned int Np = libmesh_cast_int<unsigned int>
926  (parameters.size());
927  const unsigned int Nq = libmesh_cast_int<unsigned int>
928  (qoi.size());
929 
930  // For each quantity of interest q, the parameter sensitivity
931  // Hessian is defined as q''_{kl} = {d^2 q}/{d p_k d p_l}.
932  // Given a vector of parameter perturbation weights w_l, this
933  // function evaluates the hessian-vector product sum_l(q''_{kl}*w_l)
934  //
935  // We calculate it from values and partial derivatives of the
936  // quantity of interest function Q, solution u, adjoint solution z,
937  // parameter sensitivity adjoint solutions z^l, and residual R, as:
938  //
939  // sum_l(q''_{kl}*w_l) =
940  // sum_l(w_l * Q''_{kl}) + Q''_{uk}(u)*(sum_l(w_l u'_l)) -
941  // R'_k(u, sum_l(w_l*z^l)) - R'_{uk}(u,z)*(sum_l(w_l u'_l) -
942  // sum_l(w_l*R''_{kl}(u,z))
943  //
944  // See the adjoints model document for more details.
945 
946  // We first do an adjoint solve to get z for each quantity of
947  // interest
948  // if we havent already or dont have an initial condition for the adjoint
949  if (!this->is_adjoint_already_solved())
950  {
951  this->adjoint_solve(qoi_indices);
952  }
953 
954  // Get ready to fill in senstivities:
955  sensitivities.allocate_data(qoi_indices, *this, parameters);
956 
957  // We can't solve for all the solution sensitivities u'_l or for all
958  // of the parameter sensitivity adjoint solutions z^l without
959  // requiring O(Nq*Np) linear solves. So we'll solve directly for their
960  // weighted sum - this is just O(Nq) solves.
961 
962  // First solve for sum_l(w_l u'_l).
963  this->weighted_sensitivity_solve(parameters, vector);
964 
965  // Then solve for sum_l(w_l z^l).
966  this->weighted_sensitivity_adjoint_solve(parameters, vector, qoi_indices);
967 
968  for (unsigned int k=0; k != Np; ++k)
969  {
970  // We approximate sum_l(w_l * Q''_{kl}) with a central
971  // differencing pertubation:
972  // sum_l(w_l * Q''_{kl}) ~=
973  // (Q(p + dp*w_l*e_l + dp*e_k) - Q(p - dp*w_l*e_l + dp*e_k) -
974  // Q(p + dp*w_l*e_l - dp*e_k) + Q(p - dp*w_l*e_l - dp*e_k))/(4*dp^2)
975 
976  // The sum(w_l*R''_kl) term requires the same sort of pertubation,
977  // and so we subtract it in at the same time:
978  // sum_l(w_l * R''_{kl}) ~=
979  // (R(p + dp*w_l*e_l + dp*e_k) - R(p - dp*w_l*e_l + dp*e_k) -
980  // R(p + dp*w_l*e_l - dp*e_k) + R(p - dp*w_l*e_l - dp*e_k))/(4*dp^2)
981 
982  ParameterVector oldparameters, parameterperturbation;
983  parameters.deep_copy(oldparameters);
984  vector.deep_copy(parameterperturbation);
985  parameterperturbation *= delta_p;
986  parameters += parameterperturbation;
987 
988  Number old_parameter = *parameters[k];
989 
990  *parameters[k] = old_parameter + delta_p;
991  this->assemble_qoi(qoi_indices);
992  this->assembly(true, false);
993  this->rhs->close();
994  std::vector<Number> partial2q_term = this->qoi;
995  std::vector<Number> partial2R_term(this->qoi.size());
996  for (unsigned int i=0; i != Nq; ++i)
997  if (qoi_indices.has_index(i))
998  partial2R_term[i] = this->rhs->dot(this->get_adjoint_solution(i));
999 
1000  *parameters[k] = old_parameter - delta_p;
1001  this->assemble_qoi(qoi_indices);
1002  this->assembly(true, false);
1003  this->rhs->close();
1004  for (unsigned int i=0; i != Nq; ++i)
1005  if (qoi_indices.has_index(i))
1006  {
1007  partial2q_term[i] -= this->qoi[i];
1008  partial2R_term[i] -= this->rhs->dot(this->get_adjoint_solution(i));
1009  }
1010 
1011  oldparameters.value_copy(parameters);
1012  parameterperturbation *= -1.0;
1013  parameters += parameterperturbation;
1014 
1015  // Re-center old_parameter, which may be affected by vector
1016  old_parameter = *parameters[k];
1017 
1018  *parameters[k] = old_parameter + delta_p;
1019  this->assemble_qoi(qoi_indices);
1020  this->assembly(true, false);
1021  this->rhs->close();
1022  for (unsigned int i=0; i != Nq; ++i)
1023  if (qoi_indices.has_index(i))
1024  {
1025  partial2q_term[i] -= this->qoi[i];
1026  partial2R_term[i] -= this->rhs->dot(this->get_adjoint_solution(i));
1027  }
1028 
1029  *parameters[k] = old_parameter - delta_p;
1030  this->assemble_qoi(qoi_indices);
1031  this->assembly(true, false);
1032  this->rhs->close();
1033  for (unsigned int i=0; i != Nq; ++i)
1034  if (qoi_indices.has_index(i))
1035  {
1036  partial2q_term[i] += this->qoi[i];
1037  partial2R_term[i] += this->rhs->dot(this->get_adjoint_solution(i));
1038  }
1039 
1040  for (unsigned int i=0; i != Nq; ++i)
1041  if (qoi_indices.has_index(i))
1042  {
1043  partial2q_term[i] /= (4. * delta_p * delta_p);
1044  partial2R_term[i] /= (4. * delta_p * delta_p);
1045  }
1046 
1047  for (unsigned int i=0; i != Nq; ++i)
1048  if (qoi_indices.has_index(i))
1049  sensitivities[i][k] = partial2q_term[i] - partial2R_term[i];
1050 
1051  // We get (partial q / partial u), R, and
1052  // (partial R / partial u) from the user, but centrally
1053  // difference to get q_uk, R_k, and R_uk terms:
1054  // (partial R / partial k)
1055  // 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)
1056  // (partial^2 q / partial u partial k)
1057  // q_uk = (q_u(p+dp*e_k) - q_u(p-dp*e_k))/(2*dp)
1058  // (partial^2 R / partial u partial k)
1059  // 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)
1060 
1061  // To avoid creating Nq temporary vectors for q_uk or R_uk, we add
1062  // subterms to the sensitivities output one by one.
1063  //
1064  // FIXME: this is probably a bad order of operations for
1065  // controlling floating point error.
1066 
1067  *parameters[k] = old_parameter + delta_p;
1068  this->assembly(true, true);
1069  this->rhs->close();
1070  this->matrix->close();
1071  this->assemble_qoi_derivative(qoi_indices);
1072 
1073  this->matrix->vector_mult(*tempvec, this->get_weighted_sensitivity_solution());
1074 
1075  for (unsigned int i=0; i != Nq; ++i)
1076  if (qoi_indices.has_index(i))
1077  {
1078  this->get_adjoint_rhs(i).close();
1079  sensitivities[i][k] += (this->get_adjoint_rhs(i).dot(this->get_weighted_sensitivity_solution()) -
1080  this->rhs->dot(this->get_weighted_sensitivity_adjoint_solution(i)) -
1081  this->get_adjoint_solution(i).dot(*tempvec)) / (2.*delta_p);
1082  }
1083 
1084  *parameters[k] = old_parameter - delta_p;
1085  this->assembly(true, true);
1086  this->rhs->close();
1087  this->matrix->close();
1088  this->assemble_qoi_derivative(qoi_indices);
1089 
1090  this->matrix->vector_mult(*tempvec, this->get_weighted_sensitivity_solution());
1091 
1092  for (unsigned int i=0; i != Nq; ++i)
1093  if (qoi_indices.has_index(i))
1094  {
1095  this->get_adjoint_rhs(i).close();
1096  sensitivities[i][k] += (-this->get_adjoint_rhs(i).dot(this->get_weighted_sensitivity_solution()) +
1097  this->rhs->dot(this->get_weighted_sensitivity_adjoint_solution(i)) +
1098  this->get_adjoint_solution(i).dot(*tempvec)) / (2.*delta_p);
1099  }
1100  }
1101 
1102  // All parameters have been reset.
1103  // Don't leave the qoi or system changed - principle of least
1104  // surprise.
1105  this->assembly(true, true);
1106  this->rhs->close();
1107  this->matrix->close();
1108  this->assemble_qoi(qoi_indices);
1109 }
1110 
1111 
1112 
1114  (const QoISet& qoi_indices,
1115  const ParameterVector& parameters,
1116  SensitivityData& sensitivities)
1117 {
1118  // We currently get partial derivatives via finite differencing
1119  const Real delta_p = TOLERANCE;
1120 
1121  // We'll use one temporary vector for matrix-vector-vector products
1122  AutoPtr<NumericVector<Number> > tempvec = this->solution->zero_clone();
1123 
1124  // And another temporary vector to hold a copy of the true solution
1125  // so we can safely perturb this->solution.
1126  AutoPtr<NumericVector<Number> > oldsolution = this->solution->clone();
1127 
1128  const unsigned int Np = libmesh_cast_int<unsigned int>
1129  (parameters.size());
1130  const unsigned int Nq = libmesh_cast_int<unsigned int>
1131  (qoi.size());
1132 
1133  // For each quantity of interest q, the parameter sensitivity
1134  // Hessian is defined as q''_{kl} = {d^2 q}/{d p_k d p_l}.
1135  //
1136  // We calculate it from values and partial derivatives of the
1137  // quantity of interest function Q, solution u, adjoint solution z,
1138  // and residual R, as:
1139  //
1140  // q''_{kl} =
1141  // Q''_{kl} + Q''_{uk}(u)*u'_l + Q''_{ul}(u) * u'_k +
1142  // Q''_{uu}(u)*u'_k*u'_l -
1143  // R''_{kl}(u,z) -
1144  // R''_{uk}(u,z)*u'_l - R''_{ul}(u,z)*u'_k -
1145  // R''_{uu}(u,z)*u'_k*u'_l
1146  //
1147  // See the adjoints model document for more details.
1148 
1149  // We first do an adjoint solve to get z for each quantity of
1150  // interest
1151  // if we havent already or dont have an initial condition for the adjoint
1152  if (!this->is_adjoint_already_solved())
1153  {
1154  this->adjoint_solve(qoi_indices);
1155  }
1156 
1157  // And a sensitivity solve to get u_k for each parameter
1158  this->sensitivity_solve(parameters);
1159 
1160  // Get ready to fill in second derivatives:
1161  sensitivities.allocate_hessian_data(qoi_indices, *this, parameters);
1162 
1163  for (unsigned int k=0; k != Np; ++k)
1164  {
1165  Number old_parameterk = *parameters[k];
1166 
1167  // The Hessian is symmetric, so we just calculate the lower
1168  // triangle and the diagonal, and we get the upper triangle from
1169  // the transpose of the lower
1170 
1171  for (unsigned int l=0; l != k+1; ++l)
1172  {
1173  // The second partial derivatives with respect to parameters
1174  // are all calculated via a central finite difference
1175  // stencil:
1176  // F''_{kl} ~= (F(p+dp*e_k+dp*e_l) - F(p+dp*e_k-dp*e_l) -
1177  // F(p-dp*e_k+dp*e_l) + F(p-dp*e_k-dp*e_l))/(4*dp^2)
1178  // We will add Q''_{kl}(u) and subtract R''_{kl}(u,z) at the
1179  // same time.
1180  //
1181  // We have to be careful with the perturbations to handle
1182  // the k=l case
1183 
1184  Number old_parameterl = *parameters[l];
1185 
1186  *parameters[k] += delta_p;
1187  *parameters[l] += delta_p;
1188  this->assemble_qoi(qoi_indices);
1189  this->assembly(true, false);
1190  this->rhs->close();
1191  std::vector<Number> partial2q_term = this->qoi;
1192  std::vector<Number> partial2R_term(this->qoi.size());
1193  for (unsigned int i=0; i != Nq; ++i)
1194  if (qoi_indices.has_index(i))
1195  partial2R_term[i] = this->rhs->dot(this->get_adjoint_solution(i));
1196 
1197  *parameters[l] -= 2.*delta_p;
1198  this->assemble_qoi(qoi_indices);
1199  this->assembly(true, false);
1200  this->rhs->close();
1201  for (unsigned int i=0; i != Nq; ++i)
1202  if (qoi_indices.has_index(i))
1203  {
1204  partial2q_term[i] -= this->qoi[i];
1205  partial2R_term[i] -= this->rhs->dot(this->get_adjoint_solution(i));
1206  }
1207 
1208  *parameters[k] -= 2.*delta_p;
1209  this->assemble_qoi(qoi_indices);
1210  this->assembly(true, false);
1211  this->rhs->close();
1212  for (unsigned int i=0; i != Nq; ++i)
1213  if (qoi_indices.has_index(i))
1214  {
1215  partial2q_term[i] += this->qoi[i];
1216  partial2R_term[i] += this->rhs->dot(this->get_adjoint_solution(i));
1217  }
1218 
1219  *parameters[l] += 2.*delta_p;
1220  this->assemble_qoi(qoi_indices);
1221  this->assembly(true, false);
1222  this->rhs->close();
1223  for (unsigned int i=0; i != Nq; ++i)
1224  if (qoi_indices.has_index(i))
1225  {
1226  partial2q_term[i] -= this->qoi[i];
1227  partial2R_term[i] -= this->rhs->dot(this->get_adjoint_solution(i));
1228  partial2q_term[i] /= (4. * delta_p * delta_p);
1229  partial2R_term[i] /= (4. * delta_p * delta_p);
1230  }
1231 
1232  for (unsigned int i=0; i != Nq; ++i)
1233  if (qoi_indices.has_index(i))
1234  {
1235  Number current_terms = partial2q_term[i] - partial2R_term[i];
1236  sensitivities.second_derivative(i,k,l) += current_terms;
1237  if (k != l)
1238  sensitivities.second_derivative(i,l,k) += current_terms;
1239  }
1240 
1241  // Don't leave the parameters perturbed
1242  *parameters[l] = old_parameterl;
1243  *parameters[k] = old_parameterk;
1244  }
1245 
1246  // We get (partial q / partial u) and
1247  // (partial R / partial u) from the user, but centrally
1248  // difference to get q_uk and R_uk terms:
1249  // (partial^2 q / partial u partial k)
1250  // q_uk*u'_l = (q_u(p+dp*e_k)*u'_l - q_u(p-dp*e_k)*u'_l)/(2*dp)
1251  // 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)
1252  //
1253  // To avoid creating Nq temporary vectors, we add these
1254  // subterms to the sensitivities output one by one.
1255  //
1256  // FIXME: this is probably a bad order of operations for
1257  // controlling floating point error.
1258 
1259  *parameters[k] = old_parameterk + delta_p;
1260  this->assembly(false, true);
1261  this->matrix->close();
1262  this->assemble_qoi_derivative(qoi_indices);
1263 
1264  for (unsigned int l=0; l != Np; ++l)
1265  {
1266  this->matrix->vector_mult(*tempvec, this->get_sensitivity_solution(l));
1267  for (unsigned int i=0; i != Nq; ++i)
1268  if (qoi_indices.has_index(i))
1269  {
1270  this->get_adjoint_rhs(i).close();
1271  Number current_terms =
1272  (this->get_adjoint_rhs(i).dot(this->get_sensitivity_solution(l)) -
1273  tempvec->dot(this->get_adjoint_solution(i))) / (2.*delta_p);
1274  sensitivities.second_derivative(i,k,l) += current_terms;
1275 
1276  // We use the _uk terms twice; symmetry lets us reuse
1277  // these calculations for the _ul terms.
1278 
1279  sensitivities.second_derivative(i,l,k) += current_terms;
1280  }
1281  }
1282 
1283  *parameters[k] = old_parameterk - delta_p;
1284  this->assembly(false, true);
1285  this->matrix->close();
1286  this->assemble_qoi_derivative(qoi_indices);
1287 
1288  for (unsigned int l=0; l != Np; ++l)
1289  {
1290  this->matrix->vector_mult(*tempvec, this->get_sensitivity_solution(l));
1291  for (unsigned int i=0; i != Nq; ++i)
1292  if (qoi_indices.has_index(i))
1293  {
1294  this->get_adjoint_rhs(i).close();
1295  Number current_terms =
1296  (-this->get_adjoint_rhs(i).dot(this->get_sensitivity_solution(l)) +
1297  tempvec->dot(this->get_adjoint_solution(i))) / (2.*delta_p);
1298  sensitivities.second_derivative(i,k,l) += current_terms;
1299 
1300  // We use the _uk terms twice; symmetry lets us reuse
1301  // these calculations for the _ul terms.
1302 
1303  sensitivities.second_derivative(i,l,k) += current_terms;
1304  }
1305  }
1306 
1307  // Don't leave the parameter perturbed
1308  *parameters[k] = old_parameterk;
1309 
1310  // Our last remaining terms are -R_uu(u,z)*u_k*u_l and
1311  // Q_uu(u)*u_k*u_l
1312  //
1313  // We take directional central finite differences of R_u and Q_u
1314  // to approximate these terms, e.g.:
1315  //
1316  // Q_uu(u)*u_k ~= (Q_u(u+dp*u_k) - Q_u(u-dp*u_k))/(2*dp)
1317 
1318  *this->solution = this->get_sensitivity_solution(k);
1319  *this->solution *= delta_p;
1320  *this->solution += *oldsolution;
1321  this->assembly(false, true);
1322  this->matrix->close();
1323  this->assemble_qoi_derivative(qoi_indices);
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  this->assembly(false, true);
1352  this->matrix->close();
1353  this->assemble_qoi_derivative(qoi_indices);
1354 
1355  for (unsigned int l=0; l != k+1; ++l)
1356  {
1357  this->matrix->vector_mult(*tempvec, this->get_sensitivity_solution(l));
1358  for (unsigned int i=0; i != Nq; ++i)
1359  if (qoi_indices.has_index(i))
1360  {
1361  this->get_adjoint_rhs(i).close();
1362  Number current_terms =
1363  (-this->get_adjoint_rhs(i).dot(this->get_sensitivity_solution(l)) +
1364  tempvec->dot(this->get_adjoint_solution(i))) / (2.*delta_p);
1365  sensitivities.second_derivative(i,k,l) += current_terms;
1366  if (k != l)
1367  sensitivities.second_derivative(i,l,k) += current_terms;
1368  }
1369  }
1370 
1371  // Don't leave the solution perturbed
1372  *this->solution = *oldsolution;
1373  }
1374 
1375  // All parameters have been reset.
1376  // Don't leave the qoi or system changed - principle of least
1377  // surprise.
1378  this->assembly(true, true);
1379  this->rhs->close();
1380  this->matrix->close();
1381  this->assemble_qoi(qoi_indices);
1382 }
1383 
1384 
1385 
1387 {
1388  return LinearSolver<Number>::build(this->comm()).release();
1389 }
1390 
1391 
1392 
1393 std::pair<unsigned int, Real> ImplicitSystem::get_linear_solve_parameters() const
1394 {
1395  return std::make_pair(this->get_equation_systems().parameters.get<unsigned int>("linear solver maximum iterations"),
1396  this->get_equation_systems().parameters.get<Real>("linear solver tolerance"));
1397 }
1398 
1399 
1400 
1402 {
1403  delete s;
1404 }
1405 
1406 } // namespace libMesh

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:05 UTC

Hosted By:
SourceForge.net Logo