MAST
nonlinear_implicit_assembly.cpp
Go to the documentation of this file.
1 /*
2  * MAST: Multidisciplinary-design Adaptation and Sensitivity Toolkit
3  * Copyright (C) 2013-2019 Manav Bhatia
4  *
5  * This library is free software; you can redistribute it and/or
6  * modify it under the terms of the GNU Lesser General Public
7  * License as published by the Free Software Foundation; either
8  * version 2.1 of the License, or (at your option) any later version.
9  *
10  * This library is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13  * Lesser General Public License for more details.
14  *
15  * You should have received a copy of the GNU Lesser General Public
16  * License along with this library; if not, write to the Free Software
17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
18  */
19 
20 // MAST includes
25 #include "base/nonlinear_system.h"
28 #include "numerics/utility.h"
29 #include "mesh/geom_elem.h"
30 
31 // libMesh includes
32 #include "libmesh/nonlinear_solver.h"
33 #include "libmesh/numeric_vector.h"
34 #include "libmesh/sparse_matrix.h"
35 #include "libmesh/dof_map.h"
36 
37 
38 
41 _post_assembly (nullptr),
42 _res_l2_norm (0.),
43 _first_iter_res_l2_norm (-1.) {
44 
45 }
46 
47 
48 
50 
51 }
52 
53 
54 
55 void
58 
59  _post_assembly = &post;
60 }
61 
62 
63 
64 void
66 residual_and_jacobian (const libMesh::NumericVector<Real>& X,
67  libMesh::NumericVector<Real>* R,
68  libMesh::SparseMatrix<Real>* J,
69  libMesh::NonlinearImplicitSystem& S) {
70 
71  libmesh_assert(_system);
72  libmesh_assert(_discipline);
73  libmesh_assert(_elem_ops);
74 
75  MAST::NonlinearSystem& nonlin_sys = _system->system();
76 
77  // make sure that the system for which this object was created,
78  // and the system passed through the function call are the same
79  libmesh_assert_equal_to(&S, &(nonlin_sys));
80 
81  if (R) R->zero();
82  if (J) J->zero();
83 
84  // iterate over each element, initialize it and get the relevant
85  // analysis quantities
86  RealVectorX vec, sol;
87  RealMatrixX mat;
88 
89  std::vector<libMesh::dof_id_type> dof_indices;
90  const libMesh::DofMap& dof_map = _system->system().get_dof_map();
91 
92 
93  std::unique_ptr<libMesh::NumericVector<Real> > localized_solution;
94  localized_solution.reset(build_localized_vector(nonlin_sys,
95  X).release());
96 
97 
98  // if a solution function is attached, initialize it
99  if (_sol_function)
100  _sol_function->init( X);
101 
102 
103  libMesh::MeshBase::const_element_iterator el =
104  nonlin_sys.get_mesh().active_local_elements_begin();
105  const libMesh::MeshBase::const_element_iterator end_el =
106  nonlin_sys.get_mesh().active_local_elements_end();
107 
109  ops = dynamic_cast<MAST::NonlinearImplicitAssemblyElemOperations&>(*_elem_ops);
110 
111  for ( ; el != end_el; ++el) {
112 
113  const libMesh::Elem* elem = *el;
114 
115  dof_map.dof_indices (elem, dof_indices);
116 
117  MAST::GeomElem geom_elem;
118  ops.set_elem_data(elem->dim(), *elem, geom_elem);
119  geom_elem.init(*elem, *_system);
120 
121  ops.init(geom_elem);
122 
123  // get the solution
124  unsigned int ndofs = (unsigned int)dof_indices.size();
125  sol.setZero(ndofs);
126  vec.setZero(ndofs);
127  mat.setZero(ndofs, ndofs);
128 
129  for (unsigned int i=0; i<dof_indices.size(); i++)
130  sol(i) = (*localized_solution)(dof_indices[i]);
131 
132  ops.set_elem_solution(sol);
133 
134 
135 // if (_sol_function)
136 // physics_elem->attach_active_solution_function(*_sol_function);
137 
138  //_check_element_numerical_jacobian(*physics_elem, sol);
139 
140  // perform the element level calculations
141  ops.elem_calculations(J!=nullptr?true:false,
142  vec, mat);
143 
144 // physics_elem->detach_active_solution_function();
145 
146  ops.clear_elem();
147 
148  // copy to the libMesh matrix for further processing
149  DenseRealVector v;
150  DenseRealMatrix m;
151  if (R)
152  MAST::copy(v, vec);
153  if (J)
154  MAST::copy(m, mat);
155 
156  // constrain the quantities to account for hanging dofs,
157  // Dirichlet constraints, etc.
158  if (R && J)
159  dof_map.constrain_element_matrix_and_vector(m, v, dof_indices);
160  else if (R)
161  dof_map.constrain_element_vector(v, dof_indices);
162  else
163  dof_map.constrain_element_matrix(m, dof_indices);
164 
165  // add to the global matrices
166  if (R) R->add_vector(v, dof_indices);
167  if (J) J->add_matrix(m, dof_indices);
168  dof_indices.clear();
169  }
170 
171 
172  // add the point loads if any in the discipline
173  if (R && _discipline->point_loads().size()) {
174 
176  loads = _discipline->point_loads();
177 
178  vec = RealVectorX::Zero(_system->n_vars());
179 
180  MAST::PointLoadSetType::const_iterator
181  it = loads.begin(),
182  end = loads.end();
183 
184  const libMesh::dof_id_type
185  first_dof = dof_map.first_dof(nonlin_sys.comm().rank()),
186  last_dof = dof_map.last_dof(nonlin_sys.comm().rank());
187 
188  for ( ; it != end; it++) {
189 
190  // get the point load function
192  &func = (*it)->get<MAST::FieldFunction<RealVectorX>>("load");
193 
194  // get the nodes on which this object defines the load
195  const std::set<const libMesh::Node*>
196  nodes = (*it)->get_nodes();
197 
198  std::set<const libMesh::Node*>::const_iterator
199  n_it = nodes.begin(),
200  n_end = nodes.end();
201 
202  for (; n_it != n_end; n_it++) {
203 
204  // load at the node
205  vec.setZero();
206  func(**n_it, nonlin_sys.time, vec);
207  // multiply with -1 to be consistent with res(X) = 0, which
208  // requires taking the force vector on RHS to the LHS
209  vec *= -1.;
210 
211  dof_map.dof_indices(*n_it, dof_indices);
212 
213  libmesh_assert_equal_to(dof_indices.size(), vec.rows());
214 
215  // zero the components of the vector if they do not
216  // belong to this processor
217  for (unsigned int i=0; i<dof_indices.size(); i++)
218  if (dof_indices[i] < first_dof ||
219  dof_indices[i] >= last_dof)
220  vec(i) = 0.;
221 
222  DenseRealVector v;
223  MAST::copy(v, vec);
224 
225  dof_map.constrain_element_vector(v, dof_indices);
226  R->add_vector(v, dof_indices);
227  dof_indices.clear();
228  }
229  }
230  }
231 
232  // call the post assembly object, if provided by user
233  if (_post_assembly)
234  _post_assembly->post_assembly(X, R, J, S);
235 
236 
237  // if a solution function is attached, clear it
238  if (_sol_function)
239  _sol_function->clear();
240 
241  if (R) {
242 
243  R->close();
244  _res_l2_norm = R->l2_norm();
245  if (_first_iter_res_l2_norm < 0.)
247  }
248  if (J && close_matrix) J->close();
249 }
250 
251 
252 
253 
254 void
256 linearized_jacobian_solution_product (const libMesh::NumericVector<Real>& X,
257  const libMesh::NumericVector<Real>& dX,
258  libMesh::NumericVector<Real>& JdX,
259  libMesh::NonlinearImplicitSystem& S) {
260 
261  libmesh_assert(_system);
262  libmesh_assert(_discipline);
263  libmesh_assert(_elem_ops);
264 
265  // zero the solution vector
266  JdX.zero();
267 
268  MAST::NonlinearSystem& nonlin_sys = _system->system();
269 
270  // make sure that the system for which this object was created,
271  // and the system passed through the function call are the same
272  libmesh_assert_equal_to(&S, &(nonlin_sys));
273 
274  // iterate over each element, initialize it and get the relevant
275  // analysis quantities
276  RealVectorX vec, sol, dsol;
277  RealMatrixX mat;
278 
279  std::vector<libMesh::dof_id_type> dof_indices;
280  const libMesh::DofMap& dof_map = _system->system().get_dof_map();
281 
282 
283  std::unique_ptr<libMesh::NumericVector<Real> >
284  localized_solution,
285  localized_perturbed_solution;
286 
287  localized_solution.reset(build_localized_vector(nonlin_sys,
288  X).release());
289  localized_perturbed_solution.reset(build_localized_vector(nonlin_sys,
290  dX).release());
291 
292 
293  // if a solution function is attached, initialize it
294  if (_sol_function)
295  _sol_function->init( X);
296 
297 
298  libMesh::MeshBase::const_element_iterator el =
299  nonlin_sys.get_mesh().active_local_elements_begin();
300  const libMesh::MeshBase::const_element_iterator end_el =
301  nonlin_sys.get_mesh().active_local_elements_end();
302 
304  ops = dynamic_cast<MAST::NonlinearImplicitAssemblyElemOperations&>(*_elem_ops);
305 
306  for ( ; el != end_el; ++el) {
307 
308  const libMesh::Elem* elem = *el;
309 
310  dof_map.dof_indices (elem, dof_indices);
311 
312  MAST::GeomElem geom_elem;
313  ops.set_elem_data(elem->dim(), *elem, geom_elem);
314  geom_elem.init(*elem, *_system);
315 
316  ops.init(geom_elem);
317 
318  // get the solution
319  unsigned int ndofs = (unsigned int)dof_indices.size();
320  sol.setZero(ndofs);
321  dsol.setZero(ndofs);
322  vec.setZero(ndofs);
323  mat.setZero(ndofs, ndofs);
324 
325  for (unsigned int i=0; i<dof_indices.size(); i++) {
326  sol (i) = (*localized_solution) (dof_indices[i]);
327  dsol(i) = (*localized_perturbed_solution)(dof_indices[i]);
328  }
329 
330  ops.set_elem_solution(sol);
331  ops.set_elem_perturbed_solution(dsol);
332 
333 // if (_sol_function)
334 // physics_elem->attach_active_solution_function(*_sol_function);
335 
336  //_check_element_numerical_jacobian(*physics_elem, sol);
337 
338  // perform the element level calculations
340 
341  //physics_elem->detach_active_solution_function();
342  ops.clear_elem();
343 
344  // copy to the libMesh matrix for further processing
345  DenseRealVector v;
346  MAST::copy(v, vec);
347 
348  // constrain the quantities to account for hanging dofs,
349  // Dirichlet constraints, etc.
350  dof_map.constrain_element_vector(v, dof_indices);
351 
352  // add to the global matrices
353  JdX.add_vector(v, dof_indices);
354  dof_indices.clear();
355  }
356 
357 
358  // if a solution function is attached, clear it
359  if (_sol_function)
360  _sol_function->clear();
361 
362  JdX.close();
363 }
364 
365 
366 
367 void
369 second_derivative_dot_solution_assembly (const libMesh::NumericVector<Real>& X,
370  const libMesh::NumericVector<Real>& dX,
371  libMesh::SparseMatrix<Real>& d_JdX_dX,
372  libMesh::NonlinearImplicitSystem& S) {
373 
374  libmesh_assert(_system);
375  libmesh_assert(_discipline);
376  libmesh_assert(_elem_ops);
377 
378  // zero the matrix
379  d_JdX_dX.zero();
380 
381  MAST::NonlinearSystem& nonlin_sys = _system->system();
382 
383  // make sure that the system for which this object was created,
384  // and the system passed through the function call are the same
385  libmesh_assert_equal_to(&S, &(nonlin_sys));
386 
387  // iterate over each element, initialize it and get the relevant
388  // analysis quantities
389  RealVectorX sol, dsol;
390  RealMatrixX mat;
391 
392  std::vector<libMesh::dof_id_type> dof_indices;
393  const libMesh::DofMap& dof_map = _system->system().get_dof_map();
394 
395 
396  std::unique_ptr<libMesh::NumericVector<Real> >
397  localized_solution,
398  localized_perturbed_solution;
399 
400  localized_solution.reset(build_localized_vector(nonlin_sys,
401  X).release());
402  localized_perturbed_solution.reset(build_localized_vector(nonlin_sys,
403  dX).release());
404 
405 
406  // if a solution function is attached, initialize it
407  if (_sol_function)
408  _sol_function->init( X);
409 
410 
411  libMesh::MeshBase::const_element_iterator el =
412  nonlin_sys.get_mesh().active_local_elements_begin();
413  const libMesh::MeshBase::const_element_iterator end_el =
414  nonlin_sys.get_mesh().active_local_elements_end();
415 
417  ops = dynamic_cast<MAST::NonlinearImplicitAssemblyElemOperations&>(*_elem_ops);
418 
419  for ( ; el != end_el; ++el) {
420 
421  const libMesh::Elem* elem = *el;
422 
423  dof_map.dof_indices (elem, dof_indices);
424 
425  MAST::GeomElem geom_elem;
426  ops.set_elem_data(elem->dim(), *elem, geom_elem);
427  geom_elem.init(*elem, *_system);
428 
429  ops.init(geom_elem);
430 
431  // get the solution
432  unsigned int ndofs = (unsigned int)dof_indices.size();
433  sol.setZero(ndofs);
434  dsol.setZero(ndofs);
435  mat.setZero(ndofs, ndofs);
436 
437  for (unsigned int i=0; i<dof_indices.size(); i++) {
438  sol (i) = (*localized_solution) (dof_indices[i]);
439  dsol(i) = (*localized_perturbed_solution)(dof_indices[i]);
440  }
441 
442  ops.set_elem_solution(sol);
444 
445 // if (_sol_function)
446 // physics_elem->attach_active_solution_function(*_sol_function);
447 
448  // perform the element level calculations
450 
451 // physics_elem->detach_active_solution_function();
452  ops.clear_elem();
453 
454  // copy to the libMesh matrix for further processing
455  DenseRealMatrix m;
456  MAST::copy(m, mat);
457 
458  // constrain the quantities to account for hanging dofs,
459  // Dirichlet constraints, etc.
460  dof_map.constrain_element_matrix(m, dof_indices);
461 
462  // add to the global matrices
463  d_JdX_dX.add_matrix(m, dof_indices);
464  }
465 
466 
467  // if a solution function is attached, clear it
468  if (_sol_function)
469  _sol_function->clear();
470 
471  d_JdX_dX.close();
472 }
473 
474 
475 
476 
477 
478 bool
481  libMesh::NumericVector<Real>& sensitivity_rhs) {
482 
483  libmesh_assert(_system);
484  libmesh_assert(_discipline);
485  libmesh_assert(_elem_ops);
486 
487  MAST::NonlinearSystem& nonlin_sys = _system->system();
488 
489  sensitivity_rhs.zero();
490 
491  // iterate over each element, initialize it and get the relevant
492  // analysis quantities
493  RealVectorX vec, sol;
494 
495  std::vector<libMesh::dof_id_type> dof_indices;
496  const libMesh::DofMap& dof_map = nonlin_sys.get_dof_map();
497 
498 
499  std::unique_ptr<libMesh::NumericVector<Real> > localized_solution;
500  localized_solution.reset(build_localized_vector(nonlin_sys,
501  *nonlin_sys.solution).release());
502 
503  // if a solution function is attached, initialize it
504  if (_sol_function)
505  _sol_function->init( *nonlin_sys.solution);
506 
507  libMesh::MeshBase::const_element_iterator el =
508  nonlin_sys.get_mesh().active_local_elements_begin();
509  const libMesh::MeshBase::const_element_iterator end_el =
510  nonlin_sys.get_mesh().active_local_elements_end();
511 
513  ops = dynamic_cast<MAST::NonlinearImplicitAssemblyElemOperations&>(*_elem_ops);
514 
515  for ( ; el != end_el; ++el) {
516 
517  const libMesh::Elem* elem = *el;
518 
519  // no sensitivity computation assembly is neeed in these cases
520  if (_param_dependence &&
521  // if object is specified and elem does not depend on it
523  continue;
524 
525  dof_map.dof_indices (elem, dof_indices);
526 
527  MAST::GeomElem geom_elem;
528  ops.set_elem_data(elem->dim(), *elem, geom_elem);
529  geom_elem.init(*elem, *_system);
530 
531  ops.init(geom_elem);
532 
533  // get the solution
534  unsigned int ndofs = (unsigned int)dof_indices.size();
535  sol.setZero(ndofs);
536  vec.setZero(ndofs);
537 
538  for (unsigned int i=0; i<dof_indices.size(); i++)
539  sol(i) = (*localized_solution)(dof_indices[i]);
540 
541  ops.set_elem_solution(sol);
542 
543 // if (_sol_function)
544 // physics_elem->attach_active_solution_function(*_sol_function);
545 
546  ops.elem_sensitivity_calculations(f, vec);
547 
548 // physics_elem->detach_active_solution_function();
549  ops.clear_elem();
550 
551  // copy to the libMesh matrix for further processing
552  DenseRealVector v;
553  MAST::copy(v, vec);
554 
555  // constrain the quantities to account for hanging dofs,
556  // Dirichlet constraints, etc.
557  dof_map.constrain_element_vector(v, dof_indices);
558 
559  // add to the global matrices
560  sensitivity_rhs.add_vector(v, dof_indices);
561  dof_indices.clear();
562  }
563 
564  // add the point loads if any in the discipline
565  if (_discipline->point_loads().size()) {
566 
568  loads = _discipline->point_loads();
569 
570  vec = RealVectorX::Zero(_system->n_vars());
571 
572  MAST::PointLoadSetType::const_iterator
573  it = loads.begin(),
574  end = loads.end();
575 
576  const libMesh::dof_id_type
577  first_dof = dof_map.first_dof(nonlin_sys.comm().rank()),
578  last_dof = dof_map.last_dof(nonlin_sys.comm().rank());
579 
580  for ( ; it != end; it++) {
581 
582  // get the point load function
584  &func = (*it)->get<MAST::FieldFunction<RealVectorX>>("load");
585 
586  // get the nodes on which this object defines the load
587  const std::set<const libMesh::Node*>
588  nodes = (*it)->get_nodes();
589 
590  std::set<const libMesh::Node*>::const_iterator
591  n_it = nodes.begin(),
592  n_end = nodes.end();
593 
594  for (; n_it != n_end; n_it++) {
595 
596  // load at the node
597  vec.setZero();
598  func.derivative(f, **n_it, nonlin_sys.time, vec);
599  vec *= -1.;
600 
601  dof_map.dof_indices(*n_it, dof_indices);
602 
603  libmesh_assert_equal_to(dof_indices.size(), vec.rows());
604 
605  // zero the components of the vector if they do not
606  // belong to this processor
607  for (unsigned int i=0; i<dof_indices.size(); i++)
608  if (dof_indices[i] < first_dof ||
609  dof_indices[i] >= last_dof)
610  vec(i) = 0.;
611 
612  DenseRealVector v;
613  MAST::copy(v, vec);
614 
615  dof_map.constrain_element_vector(v, dof_indices);
616  sensitivity_rhs.add_vector(v, dof_indices);
617  dof_indices.clear();
618  }
619  }
620  }
621 
622  // if a solution function is attached, initialize it
623  if (_sol_function)
624  _sol_function->clear();
625 
626  sensitivity_rhs.close();
627 
628  return true;
629 }
630 
631 
MAST::AssemblyElemOperations * _elem_ops
provides assembly elem operations for use by this class
virtual void set_elem_perturbed_solution(const RealVectorX &sol)
sets the element perturbed solution
MAST::NonlinearSystem & system()
virtual bool if_elem_depends_on_parameter(const libMesh::Elem &e, const MAST::FunctionBase &p) const =0
virtual void residual_and_jacobian(const libMesh::NumericVector< Real > &X, libMesh::NumericVector< Real > *R, libMesh::SparseMatrix< Real > *J, libMesh::NonlinearImplicitSystem &S)
function that assembles the matrices and vectors quantities for nonlinear solution ...
libMesh::DenseMatrix< Real > DenseRealMatrix
Real _res_l2_norm
L2 norm of the last-assembled residual.
This class implements a system for solution of nonlinear systems.
virtual void set_elem_solution_sensitivity(const RealVectorX &sol)
sets the element solution sensitivity
bool close_matrix
flag to control the closing fo the Jacobian after assembly
virtual bool sensitivity_assemble(const MAST::FunctionBase &f, libMesh::NumericVector< Real > &sensitivity_rhs)
Assembly function.
virtual void set_elem_solution(const RealVectorX &sol)
sets the element solution
virtual void second_derivative_dot_solution_assembly(const libMesh::NumericVector< Real > &X, const libMesh::NumericVector< Real > &dX, libMesh::SparseMatrix< Real > &d_JdX_dX, libMesh::NonlinearImplicitSystem &S)
calculates .
virtual void linearized_jacobian_solution_product(const libMesh::NumericVector< Real > &X, const libMesh::NumericVector< Real > &dX, libMesh::NumericVector< Real > &JdX, libMesh::NonlinearImplicitSystem &S)
calculates the product of the Jacobian and a perturbation in solution vector .
MAST::SystemInitialization * _system
System for which this assembly is performed.
virtual void elem_sensitivity_calculations(const MAST::FunctionBase &f, RealVectorX &vec)=0
performs the element sensitivity calculations over elem, and returns the element residual sensitivity...
virtual void elem_calculations(bool if_jac, RealVectorX &vec, RealMatrixX &mat)=0
performs the element calculations over elem, and returns the element vector and matrix quantities in ...
MAST::PhysicsDisciplineBase * _discipline
PhysicsDisciplineBase object for which this class is assembling.
virtual void set_elem_data(unsigned int dim, const libMesh::Elem &ref_elem, MAST::GeomElem &elem) const =0
some analyses may want to set additional element data before initialization of the GeomElem...
virtual void post_assembly(const libMesh::NumericVector< Real > &X, libMesh::NumericVector< Real > *R, libMesh::SparseMatrix< Real > *J, libMesh::NonlinearImplicitSystem &S)=0
virtual void derivative(const MAST::FunctionBase &f, ValType &v) const
calculates the value of the function derivative and returns it in v.
libMesh::DenseVector< Real > DenseRealVector
virtual void elem_second_derivative_dot_solution_assembly(RealMatrixX &mat)=0
calculates over elem, and returns the matrix in vec .
NonlinearImplicitAssembly()
constructor associates this assembly object with the system
MAST::NonlinearImplicitAssembly::PostAssemblyOperation * _post_assembly
this object, if non-NULL is user-provided to perform actions after assembly and before returning to t...
Matrix< Real, Dynamic, Dynamic > RealMatrixX
void copy(DenseRealMatrix &m1, const RealMatrixX &m2)
Definition: utility.h:167
const MAST::PointLoadSetType & point_loads() const
void clear()
clear the solution
Matrix< Real, Dynamic, 1 > RealVectorX
virtual void init(const MAST::GeomElem &elem)=0
initializes the object for calculation of element quantities for the specified elem.
std::unique_ptr< libMesh::NumericVector< Real > > build_localized_vector(const libMesh::System &sys, const libMesh::NumericVector< Real > &global) const
localizes the parallel vector so that the local copy stores all values necessary for calculation of t...
This class acts as a wrapper around libMesh::Elem for the purpose of providing a uniform interface fo...
Definition: geom_elem.h:59
user-provided object to perform actions after assembly and before returning to the solver...
void init(const libMesh::NumericVector< Real > &sol, const libMesh::NumericVector< Real > *dsol=nullptr)
initializes the data structures to perform the interpolation function of sol.
MAST::AssemblyBase::ElemParameterDependence * _param_dependence
If provided by user, this object is used by sensitiivty analysis to check for whether or the current ...
std::set< MAST::PointLoadCondition * > PointLoadSetType
virtual void clear_elem()
clears the element initialization
virtual void init(const libMesh::Elem &elem, const MAST::SystemInitialization &sys_init)
initialize the object for the specified reference elem.
Definition: geom_elem.cpp:128
virtual void elem_linearized_jacobian_solution_product(RealVectorX &vec)=0
performs the element calculations over elem, and returns the element vector quantity in vec...
void set_post_assembly_operation(MAST::NonlinearImplicitAssembly::PostAssemblyOperation &post)
sets the PostAssemblyOperation object for use after assembly.
MAST::MeshFieldFunction * _sol_function
system solution that will be initialized before each solution
virtual ~NonlinearImplicitAssembly()
destructor resets the association of this assembly object with the system