linear_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
24 #include "libmesh/linear_solver.h"
26 #include "libmesh/numeric_vector.h" // for parameter sensitivity calcs
27 //#include "libmesh/parameter_vector.h"
28 #include "libmesh/sparse_matrix.h" // for get_transpose
29 #include "libmesh/system_subset.h"
30 
31 namespace libMesh
32 {
33 
34 
35 // ------------------------------------------------------------
36 // LinearImplicitSystem implementation
38  const std::string& name_in,
39  const unsigned int number_in) :
40 
41  Parent (es, name_in, number_in),
42  linear_solver (LinearSolver<Number>::build(es.comm())),
43  _n_linear_iterations (0),
44  _final_linear_residual (1.e20),
45  _shell_matrix(NULL),
46  _subset(NULL),
47  _subset_solve_mode(SUBSET_ZERO)
48 {
49 }
50 
51 
52 
54 {
55  // Clear data
56  this->clear();
57 }
58 
59 
60 
62 {
63  // clear the linear solver
64  linear_solver->clear();
65 
66  this->restrict_solve_to(NULL);
67 
68  // clear the parent data
69  Parent::clear();
70 }
71 
72 
73 
75 {
76  // initialize parent data
78 
79  // re-initialize the linear solver interface
80  linear_solver->clear();
81 }
82 
83 
84 
86 {
87  // re-initialize the linear solver interface
88  linear_solver->clear();
89 
90  // initialize parent data
92 }
93 
94 
95 
97  const SubsetSolveMode subset_solve_mode)
98 {
99  _subset = subset;
100  _subset_solve_mode = subset_solve_mode;
101  if(subset!=NULL)
102  {
103  libmesh_assert_equal_to (&subset->get_system(), this);
104  }
105 }
106 
107 
108 
110 {
111  if (this->assemble_before_solve)
112  // Assemble the linear system
113  this->assemble ();
114 
115  // Log how long the linear solve takes.
116  // This gets done by the LinearSolver classes now [RHS]
117  // START_LOG("solve()", "System");
118 
119  // Get a reference to the EquationSystems
120  const EquationSystems& es =
121  this->get_equation_systems();
122 
123  // Get the user-specifiied linear solver tolerance
124  const Real tol =
125  es.parameters.get<Real>("linear solver tolerance");
126 
127  // Get the user-specified maximum # of linear solver iterations
128  const unsigned int maxits =
129  es.parameters.get<unsigned int>("linear solver maximum iterations");
130 
131  if(_subset!=NULL)
132  {
133  linear_solver->restrict_solve_to(&_subset->dof_ids(),_subset_solve_mode);
134  }
135 
136  // Solve the linear system. Several cases:
137  std::pair<unsigned int, Real> rval = std::make_pair(0,0.0);
138  if(_shell_matrix)
139  // 1.) Shell matrix with or without user-supplied preconditioner.
140  rval = linear_solver->solve(*_shell_matrix, this->request_matrix("Preconditioner"), *solution, *rhs, tol, maxits);
141  else
142  // 2.) No shell matrix, with or without user-supplied preconditioner
143  rval = linear_solver->solve (*matrix, this->request_matrix("Preconditioner"), *solution, *rhs, tol, maxits);
144 
145  if(_subset!=NULL)
146  {
147  linear_solver->restrict_solve_to(NULL);
148  }
149 
150  // Store the number of linear iterations required to
151  // solve and the final residual.
152  _n_linear_iterations = rval.first;
153  _final_linear_residual = rval.second;
154 
155  // Stop logging the linear solve
156  // This gets done by the LinearSolver classes now [RHS]
157  // STOP_LOG("solve()", "System");
158 
159  // Update the system after the solve
160  this->update();
161 }
162 
163 
164 
166 {
167  _shell_matrix = shell_matrix;
168 }
169 
170 
171 /*
172 void LinearImplicitSystem::sensitivity_solve (const ParameterVector& parameters)
173 {
174  if (this->assemble_before_solve)
175  {
176  // Assemble the linear system
177  this->assemble ();
178 
179  // But now assemble right hand sides with the residual's
180  // parameter derivatives
181  this->assemble_residual_derivatives(parameters);
182  }
183 
184  // Get a reference to the EquationSystems
185  const EquationSystems& es =
186  this->get_equation_systems();
187 
188  // Get the user-specifiied linear solver tolerance
189  const Real tol =
190  es.parameters.get<Real>("sensitivity solver tolerance");
191 
192  // Get the user-specified maximum # of linear solver iterations
193  const unsigned int maxits =
194  es.parameters.get<unsigned int>("sensitivity solver maximum iterations");
195 
196  // Our iteration counts and residuals will be sums of the individual
197  // results
198  _n_linear_iterations = 0;
199  _final_linear_residual = 0.0;
200  std::pair<unsigned int, Real> rval = std::make_pair(0,0.0);
201 
202  // Solve the linear system.
203  SparseMatrix<Number> *pc = this->request_matrix("Preconditioner");
204  for (unsigned int p=0; p != parameters.size(); ++p)
205  {
206  rval = linear_solver->solve (*matrix, pc,
207  this->get_sensitivity_solution(p),
208  this->get_sensitivity_rhs(p), tol, maxits);
209  _n_linear_iterations += rval.first;
210  _final_linear_residual += rval.second;
211  }
212 
213  // Our matrix is the *negative* of the Jacobian for b-A*u, so our
214  // solutions are all inverted
215  for (unsigned int p=0; p != parameters.size(); ++p)
216  {
217  this->get_sensitivity_solution(p) *= -1.0;
218  }
219 }
220 
221 
222 
223 void LinearImplicitSystem::adjoint_solve (const QoISet &qoi_indices)
224 {
225  const unsigned int Nq = this->qoi.size();
226 
227  // We currently don't support adjoint solves of shell matrices
228  // FIXME - we should let shell matrices support
229  // vector_transpose_mult so that we can use them here.
230  if(_shell_matrix!=NULL)
231  libmesh_not_implemented();
232 
233  if (this->assemble_before_solve)
234  {
235  // Assemble the linear system
236  this->assemble ();
237 
238  // And take the adjoint
239  matrix->get_transpose(*matrix);
240 
241  // Including of any separate preconditioner
242  SparseMatrix<Number> *pc = this->request_matrix("Preconditioner");
243  if(pc)
244  pc->get_transpose(*pc);
245 
246  // But now replace the right hand sides with the quantity of
247  // interest functional's derivative
248  this->assemble_qoi_derivative(qoi_indices);
249  }
250 
251  // Get a reference to the EquationSystems
252  const EquationSystems& es =
253  this->get_equation_systems();
254 
255  // Get the user-specifiied linear solver tolerance
256  const Real tol =
257  es.parameters.get<Real>("adjoint solver tolerance");
258 
259  // Get the user-specified maximum # of linear solver iterations
260  const unsigned int maxits =
261  es.parameters.get<unsigned int>("adjoint solver maximum iterations");
262 
263  // Our iteration counts and residuals will be sums of the individual
264  // results
265  _n_linear_iterations = 0;
266  _final_linear_residual = 0.0;
267  std::pair<unsigned int, Real> rval = std::make_pair(0,0.0);
268 
269  // Solve the linear system.
270  SparseMatrix<Number> *pc = this->request_matrix("Preconditioner");
271  for (unsigned int i=0; i != Nq; ++i)
272  if (qoi_indices.has_index(i))
273  {
274  rval = linear_solver->solve (*matrix, pc,
275  this->add_adjoint_solution(i),
276  this->get_adjoint_rhs(i), tol, maxits);
277  _n_linear_iterations += rval.first;
278  _final_linear_residual += rval.second;
279  }
280 }
281 
282 
283 
284 void LinearImplicitSystem::forward_qoi_parameter_sensitivity
285  (const QoISet& qoi_indices,
286  const ParameterVector& parameters,
287  SensitivityData& sensitivities)
288 {
289  const unsigned int Np = parameters.size();
290  const unsigned int Nq = this->qoi.size();
291 
292  // An introduction to the problem:
293  //
294  // A(p)*u(p) = b(p), where x is determined implicitly.
295  // Residual R(u(p),p) := b(p) - A(p)*u(p)
296  // partial R / partial u = -A
297  //
298  // This implies that:
299  // d/dp(R) = 0
300  // (partial b / partial p) -
301  // (partial A / partial p) * u -
302  // A * (partial u / partial p) = 0
303  // A * (partial u / partial p) = (partial R / partial p)
304  // = (partial b / partial p) - (partial A / partial p) * u
305 
306  // We first solve for (partial u / partial p) for each parameter:
307  // -A * (partial u / partial p) = - (partial R / partial p)
308 
309  this->sensitivity_solve(parameters);
310 
311  // Get ready to fill in senstivities:
312  sensitivities.allocate_data(qoi_indices, *this, parameters);
313 
314  // We use the identity:
315  // dq/dp = (partial q / partial p) + (partial q / partial u) *
316  // (partial u / partial p)
317 
318  // We get (partial q / partial u) from the user
319  this->assemble_qoi_derivative(qoi_indices);
320 
321  for (unsigned int j=0; j != Np; ++j)
322  {
323  // We currently get partial derivatives via central differencing
324  Number delta_p = 1e-6;
325 
326  // (partial q / partial p) ~= (q(p+dp)-q(p-dp))/(2*dp)
327 
328  Number old_parameter = *parameters[j];
329 
330  *parameters[j] = old_parameter - delta_p;
331  this->assemble_qoi(qoi_indices);
332  std::vector<Number> qoi_minus = this->qoi;
333 
334  *parameters[j] = old_parameter + delta_p;
335  this->assemble_qoi(qoi_indices);
336  std::vector<Number> &qoi_plus = this->qoi;
337  std::vector<Number> partialq_partialp(Nq, 0);
338  for (unsigned int i=0; i != Nq; ++i)
339  if (qoi_indices.has_index(i))
340  partialq_partialp[i] = (qoi_plus[i] - qoi_minus[i]) / (2.*delta_p);
341 
342  for (unsigned int i=0; i != Nq; ++i)
343  if (qoi_indices.has_index(i))
344  sensitivities[i][j] = partialq_partialp[i] +
345  this->get_adjoint_rhs(i).dot(this->get_sensitivity_solution(i));
346  }
347 
348  // All parameters have been reset.
349  // Don't leave the qoi or system changed - principle of least
350  // surprise.
351  this->assemble();
352  this->rhs->close();
353  this->matrix->close();
354  this->assemble_qoi(qoi_indices);
355 }
356 */
357 
358 
359 
361 {
362  return linear_solver.get();
363 }
364 
365 
366 
368 {
369 }
370 
371 
372 
374  bool)
375 {
376  // Residual R(u(p),p) := A(p)*u(p) - b(p)
377  // partial R / partial u = A
378 
379  this->assemble();
380  this->rhs->close();
381  this->matrix->close();
382 
383  *(this->rhs) *= -1.0;
384  this->rhs->add_vector(*this->solution, *this->matrix);
385 }
386 
387 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo