libMesh
linear_implicit_system.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 // C++ includes
21 
22 // Local includes
23 #include "libmesh/linear_implicit_system.h"
24 #include "libmesh/linear_solver.h"
25 #include "libmesh/equation_systems.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(libmesh_nullptr),
46  _subset(libmesh_nullptr),
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 
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 
102  if (subset != libmesh_nullptr)
103  libmesh_assert_equal_to (&subset->get_system(), this);
104 }
105 
106 
107 
109 {
110  if (this->assemble_before_solve)
111  // Assemble the linear system
112  this->assemble ();
113 
114  // Get a reference to the EquationSystems
115  const EquationSystems & es =
116  this->get_equation_systems();
117 
118  // If the linear solver hasn't been initialized, we do so here.
119  if (libMesh::on_command_line("--solver_system_names"))
120  linear_solver->init((this->name()+"_").c_str());
121  else
122  linear_solver->init();
123 
124  // Get the user-specified linear solver tolerance
125  const Real tol =
126  es.parameters.get<Real>("linear solver tolerance");
127 
128  // Get the user-specified maximum # of linear solver iterations
129  const unsigned int maxits =
130  es.parameters.get<unsigned int>("linear solver maximum iterations");
131 
132  if (_subset != libmesh_nullptr)
133  linear_solver->restrict_solve_to(&_subset->dof_ids(),_subset_solve_mode);
134 
135  // Solve the linear system. Several cases:
136  std::pair<unsigned int, Real> rval = std::make_pair(0,0.0);
137  if (_shell_matrix)
138  // 1.) Shell matrix with or without user-supplied preconditioner.
139  rval = linear_solver->solve(*_shell_matrix, this->request_matrix("Preconditioner"), *solution, *rhs, tol, maxits);
140  else
141  // 2.) No shell matrix, with or without user-supplied preconditioner
142  rval = linear_solver->solve (*matrix, this->request_matrix("Preconditioner"), *solution, *rhs, tol, maxits);
143 
144  if (_subset != libmesh_nullptr)
145  linear_solver->restrict_solve_to(libmesh_nullptr);
146 
147  // Store the number of linear iterations required to
148  // solve and the final residual.
149  _n_linear_iterations = rval.first;
150  _final_linear_residual = rval.second;
151 
152  // Update the system after the solve
153  this->update();
154 }
155 
156 
157 
159 {
160  _shell_matrix = shell_matrix;
161 }
162 
163 
164 /*
165  void LinearImplicitSystem::sensitivity_solve (const ParameterVector & parameters)
166  {
167  if (this->assemble_before_solve)
168  {
169  // Assemble the linear system
170  this->assemble ();
171 
172  // But now assemble right hand sides with the residual's
173  // parameter derivatives
174  this->assemble_residual_derivatives(parameters);
175  }
176 
177  // Get a reference to the EquationSystems
178  const EquationSystems & es =
179  this->get_equation_systems();
180 
181  // Get the user-specified linear solver tolerance
182  const Real tol =
183  es.parameters.get<Real>("sensitivity solver tolerance");
184 
185  // Get the user-specified maximum # of linear solver iterations
186  const unsigned int maxits =
187  es.parameters.get<unsigned int>("sensitivity solver maximum iterations");
188 
189  // Our iteration counts and residuals will be sums of the individual
190  // results
191  _n_linear_iterations = 0;
192  _final_linear_residual = 0.0;
193  std::pair<unsigned int, Real> rval = std::make_pair(0,0.0);
194 
195  // Solve the linear system.
196  SparseMatrix<Number> * pc = this->request_matrix("Preconditioner");
197  for (std::size_t p=0; p != parameters.size(); ++p)
198  {
199  rval = linear_solver->solve (*matrix, pc,
200  this->get_sensitivity_solution(p),
201  this->get_sensitivity_rhs(p), tol, maxits);
202  _n_linear_iterations += rval.first;
203  _final_linear_residual += rval.second;
204  }
205 
206  // Our matrix is the *negative* of the Jacobian for b-A*u, so our
207  // solutions are all inverted
208  for (std::size_t p=0; p != parameters.size(); ++p)
209  {
210  this->get_sensitivity_solution(p) *= -1.0;
211  }
212  }
213 
214 
215 
216  void LinearImplicitSystem::adjoint_solve (const QoISet & qoi_indices)
217  {
218  const unsigned int Nq = this->qoi.size();
219 
220  // We currently don't support adjoint solves of shell matrices
221  // FIXME - we should let shell matrices support
222  // vector_transpose_mult so that we can use them here.
223  if (_shell_matrix!=libmesh_nullptr)
224  libmesh_not_implemented();
225 
226  if (this->assemble_before_solve)
227  {
228  // Assemble the linear system
229  this->assemble ();
230 
231  // And take the adjoint
232  matrix->get_transpose(*matrix);
233 
234  // Including of any separate preconditioner
235  SparseMatrix<Number> * pc = this->request_matrix("Preconditioner");
236  if (pc)
237  pc->get_transpose(*pc);
238 
239  // But now replace the right hand sides with the quantity of
240  // interest functional's derivative
241  this->assemble_qoi_derivative(qoi_indices);
242  }
243 
244  // Get a reference to the EquationSystems
245  const EquationSystems & es =
246  this->get_equation_systems();
247 
248  // Get the user-specified linear solver tolerance
249  const Real tol =
250  es.parameters.get<Real>("adjoint solver tolerance");
251 
252  // Get the user-specified maximum # of linear solver iterations
253  const unsigned int maxits =
254  es.parameters.get<unsigned int>("adjoint solver maximum iterations");
255 
256  // Our iteration counts and residuals will be sums of the individual
257  // results
258  _n_linear_iterations = 0;
259  _final_linear_residual = 0.0;
260  std::pair<unsigned int, Real> rval = std::make_pair(0,0.0);
261 
262  // Solve the linear system.
263  SparseMatrix<Number> * pc = this->request_matrix("Preconditioner");
264  for (unsigned int i=0; i != Nq; ++i)
265  if (qoi_indices.has_index(i))
266  {
267  rval = linear_solver->solve (*matrix, pc,
268  this->add_adjoint_solution(i),
269  this->get_adjoint_rhs(i), tol, maxits);
270  _n_linear_iterations += rval.first;
271  _final_linear_residual += rval.second;
272  }
273  }
274 
275 
276 
277  void LinearImplicitSystem::forward_qoi_parameter_sensitivity
278  (const QoISet & qoi_indices,
279  const ParameterVector & parameters,
280  SensitivityData & sensitivities)
281  {
282  const unsigned int Np = parameters.size();
283  const unsigned int Nq = this->qoi.size();
284 
285  // An introduction to the problem:
286  //
287  // A(p)*u(p) = b(p), where x is determined implicitly.
288  // Residual R(u(p),p) := b(p) - A(p)*u(p)
289  // partial R / partial u = -A
290  //
291  // This implies that:
292  // d/dp(R) = 0
293  // (partial b / partial p) -
294  // (partial A / partial p) * u -
295  // A * (partial u / partial p) = 0
296  // A * (partial u / partial p) = (partial R / partial p)
297  // = (partial b / partial p) - (partial A / partial p) * u
298 
299  // We first solve for (partial u / partial p) for each parameter:
300  // -A * (partial u / partial p) = - (partial R / partial p)
301 
302  this->sensitivity_solve(parameters);
303 
304  // Get ready to fill in sensitivities:
305  sensitivities.allocate_data(qoi_indices, *this, parameters);
306 
307  // We use the identity:
308  // dq/dp = (partial q / partial p) + (partial q / partial u) *
309  // (partial u / partial p)
310 
311  // We get (partial q / partial u) from the user
312  this->assemble_qoi_derivative(qoi_indices);
313 
314  for (unsigned int j=0; j != Np; ++j)
315  {
316  // We currently get partial derivatives via central differencing
317  Number delta_p = 1e-6;
318 
319  // (partial q / partial p) ~= (q(p+dp)-q(p-dp))/(2*dp)
320 
321  Number old_parameter = *parameters[j];
322 
323  *parameters[j] = old_parameter - delta_p;
324  this->assemble_qoi(qoi_indices);
325  std::vector<Number> qoi_minus = this->qoi;
326 
327  *parameters[j] = old_parameter + delta_p;
328  this->assemble_qoi(qoi_indices);
329  std::vector<Number> & qoi_plus = this->qoi;
330  std::vector<Number> partialq_partialp(Nq, 0);
331  for (unsigned int i=0; i != Nq; ++i)
332  if (qoi_indices.has_index(i))
333  partialq_partialp[i] = (qoi_plus[i] - qoi_minus[i]) / (2.*delta_p);
334 
335  for (unsigned int i=0; i != Nq; ++i)
336  if (qoi_indices.has_index(i))
337  sensitivities[i][j] = partialq_partialp[i] +
338  this->get_adjoint_rhs(i).dot(this->get_sensitivity_solution(i));
339  }
340 
341  // All parameters have been reset.
342  // Don't leave the qoi or system changed - principle of least
343  // surprise.
344  this->assemble();
345  this->rhs->close();
346  this->matrix->close();
347  this->assemble_qoi(qoi_indices);
348  }
349 */
350 
351 
352 
354 {
355  return linear_solver.get();
356 }
357 
358 
359 
361 {
362 }
363 
364 
365 
367  bool,
368  bool,
369  bool)
370 {
371  // Residual R(u(p),p) := A(p)*u(p) - b(p)
372  // partial R / partial u = A
373 
374  this->assemble();
375  this->rhs->close();
376  this->matrix->close();
377 
378  *(this->rhs) *= -1.0;
379  this->rhs->add_vector(*this->solution, *this->matrix);
380 }
381 
382 } // namespace libMesh
void attach_shell_matrix(ShellMatrix< Number > *shell_matrix)
This function enables the user to provide a shell matrix, i.e.
This is the EquationSystems class.
ShellMatrix< Number > * _shell_matrix
User supplies shell matrix or NULL if no shell matrix is used.
unsigned int _n_linear_iterations
The number of linear iterations required to solve the linear system Ax=b.
virtual const std::vector< unsigned int > & dof_ids() const =0
Set dofs outside the subset to zero.
virtual void add_vector(const T *v, const std::vector< numeric_index_type > &dof_indices)
Computes , where v is a pointer and each dof_indices[i] specifies where to add value v[i]...
LinearImplicitSystem(EquationSystems &es, const std::string &name, const unsigned int number)
Constructor.
const class libmesh_nullptr_t libmesh_nullptr
NumericVector< Number > * rhs
The system matrix.
const T & get(const std::string &) const
Definition: parameters.h:430
SubsetSolveMode _subset_solve_mode
If restrict-solve-to-subset mode is active, this member decides what happens with the dofs outside th...
The libMesh namespace provides an interface to certain functionality in the library.
const std::string & name() const
Definition: system.h:1998
virtual void clear() libmesh_override
Clear all the data structures associated with the system.
virtual void assembly(bool get_residual, bool get_jacobian, bool apply_heterogeneous_constraints=false, bool apply_no_constraints=false) libmesh_override
Assembles a residual in rhs and/or a jacobian in matrix, as requested.
This base class can be inherited from to provide interfaces to linear solvers from different packages...
Definition: linear_solver.h:58
virtual LinearSolver< Number > * get_linear_solver() const libmesh_override
virtual void assemble() libmesh_override
Prepares matrix and _dof_map for matrix assembly.
const System & get_system() const
Definition: system_subset.C:43
virtual void reinit() libmesh_override
Reinitializes the member data fields associated with the system, so that, e.g., assemble() may be use...
This is a base class for classes which represent subsets of the dofs of a System. ...
Definition: system_subset.h:42
This is the base class for classes which contain information related to any physical process that mig...
Definition: system.h:76
virtual void restrict_solve_to(const SystemSubset *subset, const SubsetSolveMode subset_solve_mode=SUBSET_ZERO) libmesh_override
After calling this method, any solve will be limited to the given subset.
UniquePtr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1523
virtual void reinit() libmesh_override
Reinitializes the member data fields associated with the system, so that, e.g., assemble() may be use...
virtual void close()=0
Calls the NumericVector&#39;s internal assembly routines, ensuring that the values are consistent across ...
Real _final_linear_residual
The final residual for the linear system Ax=b.
virtual void init_data() libmesh_override
Initializes new data members of the system.
const SparseMatrix< Number > * request_matrix(const std::string &mat_name) const
virtual void release_linear_solver(LinearSolver< Number > *) const libmesh_override
Releases a pointer to a linear solver acquired by this->get_linear_solver()
UniquePtr< LinearSolver< Number > > linear_solver
The LinearSolver defines the interface used to solve the linear_implicit system.
virtual void update()
Update the local values to reflect the solution on neighboring processors.
Definition: system.C:425
virtual void close()=0
Calls the SparseMatrix&#39;s internal assembly routines, ensuring that the values are consistent across p...
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual void init_data() libmesh_override
Initializes the member data fields associated with the system, so that, e.g., assemble() may be used...
bool on_command_line(const std::string &arg)
Definition: libmesh.C:921
SparseMatrix< Number > * matrix
The system matrix.
const EquationSystems & get_equation_systems() const
Definition: system.h:712
Parameters parameters
Data structure holding arbitrary parameters.
SubsetSolveMode
defines an enum for the question what happens to the dofs outside the given subset when a system is s...
const SystemSubset * _subset
The current subset on which to solve (or NULL if none).
bool assemble_before_solve
Flag which tells the system to whether or not to call the user assembly function during each call to ...
Definition: system.h:1477
virtual void solve() libmesh_override
Assembles & solves the linear system A*x=b.
virtual ~LinearImplicitSystem()
Destructor.
virtual void clear() libmesh_override
Clear all the data structures associated with the system.