GetFEM  5.4.2
getfem_model_solvers.cc
1 /*===========================================================================
2 
3  Copyright (C) 2009-2020 Yves Renard
4 
5  This file is a part of GetFEM
6 
7  GetFEM is free software; you can redistribute it and/or modify it
8  under the terms of the GNU Lesser General Public License as published
9  by the Free Software Foundation; either version 3 of the License, or
10  (at your option) any later version along with the GCC Runtime Library
11  Exception either version 3.1 or (at your option) any later version.
12  This program is distributed in the hope that it will be useful, but
13  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15  License and GCC Runtime Library Exception for more details.
16  You should have received a copy of the GNU Lesser General Public License
17  along with this program; if not, write to the Free Software Foundation,
18  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
19 
20 ===========================================================================*/
21 
22 #include "getfem/getfem_model_solvers.h"
23 #include "gmm/gmm_inoutput.h"
24 #include <iomanip>
25 
26 namespace getfem {
27 
28  #define TRACE_SOL 0
29 
30  /* ***************************************************************** */
31  /* Intermediary structure for Newton algorithms with getfem::model. */
32  /* ***************************************************************** */
33 
34  template <typename PLSOLVER>
35  class pb_base {
36  public:
37  typedef typename PLSOLVER::element_type::MATRIX MATRIX;
38  typedef typename PLSOLVER::element_type::VECTOR VECTOR;
39  typedef typename gmm::linalg_traits<VECTOR>::value_type T;
40  typedef typename gmm::number_traits<T>::magnitude_type R;
41 
42  protected:
43  PLSOLVER linear_solver;
44  const MATRIX &K;
45  VECTOR &rhs;
46  VECTOR state;
47 
48  public:
49  virtual const VECTOR &residual() const { return rhs; }
50  // A norm1 seems to be better than a norm2 (at least for contact problems).
51  virtual R residual_norm() { return gmm::vect_norm1(residual()); }
52  virtual const VECTOR &state_vector() const { return state; }
53  virtual VECTOR &state_vector() { return state; }
54  virtual R state_norm() const { return gmm::vect_norm1(state_vector()); }
55 
56  virtual void perturbation() {
57  R res = gmm::vect_norm2(state), ampl = std::max(res*R(1E-20), R(1E-50));
58  std::vector<R> V(gmm::vect_size(state));
59  gmm::fill_random(V);
60  gmm::add(gmm::scaled(V, ampl), state);
61  }
62 
63  virtual void add_to_residual(VECTOR &extra_rhs, R mult=1.) {
64  if (mult == R(1)) gmm::add(extra_rhs, rhs);
65  else if (mult != R(0)) gmm::add(gmm::scaled(extra_rhs, mult), rhs);
66  }
67 
68  virtual void linear_solve(VECTOR &dr, gmm::iteration &iter) {
69  (*linear_solver)(K, dr, rhs, iter);
70  }
71 
72  pb_base(PLSOLVER linsolv, const MATRIX &K_, VECTOR &rhs_)
73  : linear_solver(linsolv), K(K_), rhs(rhs_), state(gmm::vect_size(rhs_)) {}
74  virtual ~pb_base() {}
75  };
76 
77  /* ***************************************************************** */
78  /* Linear model problem. */
79  /* ***************************************************************** */
80  template <typename PLSOLVER>
81  class lin_model_pb : pb_base<PLSOLVER> {
82  model &md;
83  public:
84  using typename pb_base<PLSOLVER>::VECTOR;
85  using typename pb_base<PLSOLVER>::R;
86  using pb_base<PLSOLVER>::state_vector;
87  using pb_base<PLSOLVER>::linear_solve;
88  void compute_all() { md.assembly(model::BUILD_ALL); }
89  lin_model_pb(model &, PLSOLVER) = delete;
90  };
91 
92  template <>
93  lin_model_pb<rmodel_plsolver_type>::lin_model_pb
94  (model &md_, rmodel_plsolver_type linsolv)
95  : pb_base<rmodel_plsolver_type>
96  (linsolv, md_.real_tangent_matrix(), md_.set_real_rhs()),
97  md(md_)
98  { md.from_variables(state_vector()); }
99  template <>
100  lin_model_pb<cmodel_plsolver_type>::lin_model_pb
101  (model &md_, cmodel_plsolver_type linsolv)
102  : pb_base<cmodel_plsolver_type>
103  (linsolv, md_.complex_tangent_matrix(), md_.set_complex_rhs()),
104  md(md_)
105  { md.from_variables(state_vector()); }
106 
107 
108  /* ***************************************************************** */
109  /* Non-linear model problem. */
110  /* ***************************************************************** */
111  template <typename PLSOLVER>
112  class nonlin_model_pb : protected pb_base<PLSOLVER> {
113  public:
114  using typename pb_base<PLSOLVER>::VECTOR;
115  using typename pb_base<PLSOLVER>::R;
116  protected:
117  model &md;
118  abstract_newton_line_search &ls;
119  private:
120  VECTOR stateinit; // just a temp used in line_search, could also be mutable
121  public:
122  using pb_base<PLSOLVER>::residual;
123  using pb_base<PLSOLVER>::residual_norm;
124  using pb_base<PLSOLVER>::state_vector;
125  using pb_base<PLSOLVER>::state_norm;
126  using pb_base<PLSOLVER>::add_to_residual;
127  using pb_base<PLSOLVER>::perturbation;
128  using pb_base<PLSOLVER>::linear_solve;
129 
130  virtual R approx_external_load_norm() { return md.approx_external_load(); }
131 
132  virtual void compute_tangent_matrix() {
133  md.to_variables(state_vector());
134  md.assembly(model::BUILD_MATRIX);
135  }
136 
137  virtual void compute_residual() {
138  md.to_variables(state_vector());
139  md.assembly(model::BUILD_RHS);
140  }
141 
142  virtual R line_search(VECTOR &dr, const gmm::iteration &iter) {
143  size_type nit = 0;
144  gmm::resize(stateinit, gmm::vect_size(state_vector()));
145  gmm::copy(state_vector(), stateinit);
146  R alpha(1), res, /* res_init, */ R0;
147 
148  /* res_init = */ res = residual_norm();
149  // cout << "first residual = " << residual() << endl << endl;
150  R0 = gmm::real(gmm::vect_sp(dr, residual()));
151  ls.init_search(res, iter.get_iteration(), R0);
152  do {
153  alpha = ls.next_try();
154  gmm::add(stateinit, gmm::scaled(dr, alpha), state_vector());
155 
156  compute_residual();
157  res = residual_norm();
158  // cout << "residual = " << residual() << endl << endl;
159  R0 = gmm::real(gmm::vect_sp(dr, residual()));
160  ++ nit;
161  } while (!ls.is_converged(res, R0));
162 
163  if (alpha != ls.converged_value()) {
164  alpha = ls.converged_value();
165  gmm::add(stateinit, gmm::scaled(dr, alpha), state_vector());
166  res = ls.converged_residual();
167  compute_residual();
168  }
169 
170  return alpha;
171  }
172 
173  nonlin_model_pb(model &md_, abstract_newton_line_search &ls_,
174  PLSOLVER linear_solver_) = delete;
175  };
176 
177  template <>
178  nonlin_model_pb<rmodel_plsolver_type>::nonlin_model_pb
179  (model &md_, abstract_newton_line_search &ls_, rmodel_plsolver_type linsolv)
180  : pb_base<rmodel_plsolver_type>
181  (linsolv, md_.real_tangent_matrix(), md_.set_real_rhs()),
182  md(md_), ls(ls_)
183  { md.from_variables(state_vector()); }
184  template <>
185  nonlin_model_pb<cmodel_plsolver_type>::nonlin_model_pb
186  (model &md_, abstract_newton_line_search &ls_, cmodel_plsolver_type linsolv)
187  : pb_base<cmodel_plsolver_type>
188  (linsolv, md_.complex_tangent_matrix(), md_.set_complex_rhs()),
189  md(md_), ls(ls_)
190  { md.from_variables(state_vector()); }
191 
192 
193 
194  /* ***************************************************************** */
195  /* Non-linear model problem with internal variable condensation. */
196  /* ***************************************************************** */
197  template <typename PLSOLVER>
198  class nonlin_condensed_model_pb : public nonlin_model_pb<PLSOLVER> {
199  public:
200  using typename pb_base<PLSOLVER>::MATRIX;
201  using typename pb_base<PLSOLVER>::VECTOR;
202  using typename pb_base<PLSOLVER>::R;
203 
204  private:
205  using nonlin_model_pb<PLSOLVER>::md;
206  bool condensed_vars;
207  const MATRIX &internal_K;
208  VECTOR &full_rhs;
209 
210  public:
211  virtual const VECTOR &residual() const { return full_rhs; }
212  // A norm1 seems to be better than a norm2 (at least for contact problems).
213  virtual R residual_norm() { return gmm::vect_norm1(full_rhs); }
214 
215  using pb_base<PLSOLVER>::state_vector;
216  using pb_base<PLSOLVER>::state_norm;
217 
218  virtual void add_to_residual(VECTOR &extra_rhs, R mult=1.) {
219  if (mult == R(1)) gmm::add(extra_rhs, full_rhs);
220  else if (mult != R(0)) gmm::add(gmm::scaled(extra_rhs, mult), full_rhs);
221  }
222 
223  using pb_base<PLSOLVER>::perturbation;
224 
225  virtual void linear_solve(VECTOR &dr, gmm::iteration &iter) {
226  VECTOR dr0(md.nb_primary_dof(), R(0));
227  pb_base<PLSOLVER>::linear_solve(dr0, iter);
228  gmm::sub_interval I_prim(0, md.nb_primary_dof()),
229  I_intern(md.nb_primary_dof(), md.nb_internal_dof());
230  gmm::copy(dr0, gmm::sub_vector(dr, I_prim));
231  // recover solution of condensed variables: intern_K*prim_sol + intern_sol
232  gmm::mult(internal_K,
233  gmm::scaled(gmm::sub_vector(dr, I_prim), scalar_type(-1)),
234  md.internal_solution(),
235  gmm::sub_vector(dr, I_intern));
236  }
237 
238  virtual void compute_tangent_matrix() {
239  md.to_variables(state_vector(), condensed_vars);
240  md.assembly(condensed_vars ? model::BUILD_MATRIX_CONDENSED
241  : model::BUILD_MATRIX);
242  }
243 
244  virtual void compute_residual() {
245  md.to_variables(state_vector(), condensed_vars);
246  md.assembly(condensed_vars ? model::BUILD_RHS_WITH_INTERNAL // --> full_rhs
247  : model::BUILD_RHS); // ---> rhs (= full_rhs)
248  }
249 
250  nonlin_condensed_model_pb(model &md_, abstract_newton_line_search &ls_,
251  PLSOLVER linear_solver_) = delete;
252  };
253 
254  template <>
255  nonlin_condensed_model_pb<rmodel_plsolver_type>::nonlin_condensed_model_pb
256  (model &md_, abstract_newton_line_search &ls_, rmodel_plsolver_type linsolv)
257  : nonlin_model_pb<rmodel_plsolver_type>(md_, ls_, linsolv),
258  condensed_vars(md_.has_internal_variables()),
259  internal_K(md_.real_tangent_matrix(condensed_vars)),
260  full_rhs(md_.set_real_rhs(condensed_vars))
261  {
262  gmm::resize(state_vector(), md.nb_dof(condensed_vars));
263  md.from_variables(state_vector(), condensed_vars);
264  }
265 
266  template <>
267  nonlin_condensed_model_pb<cmodel_plsolver_type>::nonlin_condensed_model_pb
268  (model &md_, abstract_newton_line_search &ls_, cmodel_plsolver_type linsolv)
269  : nonlin_model_pb<cmodel_plsolver_type>(md_, ls_, linsolv),
270  condensed_vars(false), internal_K(md_.complex_tangent_matrix()),
271  full_rhs(md_.set_complex_rhs())
272  {
273  GMM_ASSERT1(false, "No support for internal variable condensation in "
274  "complex valued models.");
275  }
276 
277 
278 
279  /* ***************************************************************** */
280  /* Linear solvers. */
281  /* ***************************************************************** */
282  static rmodel_plsolver_type rdefault_linear_solver(const model &md) {
283  return default_linear_solver<model_real_sparse_matrix,
284  model_real_plain_vector>(md);
285  }
286 
287  static cmodel_plsolver_type cdefault_linear_solver(const model &md) {
288  return default_linear_solver<model_complex_sparse_matrix,
289  model_complex_plain_vector>(md);
290  }
291 
292  void default_newton_line_search::init_search(double r, size_t git, double) {
293  alpha_min_ratio = 0.9;
294  alpha_min = 1e-10;
295  alpha_max_ratio = 10.0;
296  alpha_mult = 0.25;
297  itmax = size_type(-1);
298  glob_it = git; if (git <= 1) count_pat = 0;
299  conv_alpha = alpha = alpha_old = 1.;
300  conv_r = first_res = r; it = 0;
301  count = 0;
302  max_ratio_reached = false;
303  }
304 
305  double default_newton_line_search::next_try() {
306  alpha_old = alpha; ++it;
307  // alpha *= 0.5;
308  if (alpha >= 0.4) alpha *= 0.5; else alpha *= alpha_mult;
309  return alpha_old;
310  }
311 
312  bool default_newton_line_search::is_converged(double r, double) {
313  // cout << "r = " << r << " alpha = " << alpha_old << " count_pat = " << count_pat << endl;
314  if (!max_ratio_reached && r < first_res * alpha_max_ratio) {
315  alpha_max_ratio_reached = alpha_old; r_max_ratio_reached = r;
316  it_max_ratio_reached = it; max_ratio_reached = true;
317  }
318  if (max_ratio_reached &&
319  r < r_max_ratio_reached * 0.5 &&
320  r > first_res * 1.1 && it <= it_max_ratio_reached+1) {
321  alpha_max_ratio_reached = alpha_old; r_max_ratio_reached = r;
322  it_max_ratio_reached = it;
323  }
324  if (count == 0 || r < conv_r)
325  { conv_r = r; conv_alpha = alpha_old; count = 1; }
326  if (conv_r < first_res) ++count;
327 
328  if (r < first_res * alpha_min_ratio)
329  { count_pat = 0; return true; }
330  if (count>=5 || (alpha < alpha_min && max_ratio_reached) || alpha<1e-15) {
331  if (conv_r < first_res * 0.99) count_pat = 0;
332  if (/*gmm::random() * 50. < -log(conv_alpha)-4.0 ||*/ count_pat >= 3)
333  { conv_r=r_max_ratio_reached; conv_alpha=alpha_max_ratio_reached; }
334  if (conv_r >= first_res * 0.999) count_pat++;
335  return true;
336  }
337  return false;
338  }
339 
340 
341  /* ***************************************************************** */
342  /* Computation of initial values of velocity/acceleration for */
343  /* time integration schemes. */
344  /* ***************************************************************** */
345 
346  template <typename PLSOLVER>
347  void compute_init_values(model &md, gmm::iteration &iter,
348  PLSOLVER lsolver,
349  abstract_newton_line_search &ls) {
350 
351  typename PLSOLVER::element_type::VECTOR state(md.nb_dof());
352  md.from_variables(state);
353  md.cancel_init_step();
354  md.set_time_integration(2);
355  scalar_type dt = md.get_time_step();
356  scalar_type ddt = md.get_init_time_step();
357  scalar_type t = md.get_time();
358 
359  // Solve for ddt
360  md.set_time_step(ddt);
361  gmm::iteration iter1 = iter;
362  standard_solve(md, iter1, lsolver, ls);
363  md.copy_init_time_derivative();
364 
365  // Restore the model state
366  md.set_time_step(dt);
367  md.set_time(t);
368  md.to_variables(state);
369  md.set_time_integration(1);
370  }
371 
372  /* ***************************************************************** */
373  /* Standard solve. */
374  /* ***************************************************************** */
375 
376  template <typename PLSOLVER>
377  void standard_solve(model &md, gmm::iteration &iter,
378  PLSOLVER lsolver,
379  abstract_newton_line_search &ls) {
380 
381  int time_integration = md.is_time_integration();
382  if (time_integration) {
383  if (time_integration == 1 && md.is_init_step()) {
384  compute_init_values(md, iter, lsolver, ls);
385  return;
386  }
387  md.set_time(md.get_time() + md.get_time_step());
388  md.call_init_affine_dependent_variables(time_integration);
389  }
390 
391  if (md.is_linear()) {
392  lin_model_pb<PLSOLVER> mdpb(md, lsolver);
393  mdpb.compute_all();
394  mdpb.linear_solve(mdpb.state_vector(), iter);
395  md.to_variables(mdpb.state_vector()); // copy the state vector into the model variables
396  } else {
397  std::unique_ptr<nonlin_model_pb<PLSOLVER>> mdpb;
398  if (md.has_internal_variables())
399  mdpb = std::make_unique<nonlin_condensed_model_pb<PLSOLVER>>(md, ls, lsolver);
400  else
401  mdpb = std::make_unique<nonlin_model_pb<PLSOLVER>>(md, ls, lsolver);
402  if (dynamic_cast<newton_search_with_step_control *>(&ls))
403  Newton_with_step_control(*mdpb, iter);
404  else
405  classical_Newton(*mdpb, iter);
406  md.to_variables(mdpb->state_vector()); // copy the state vector into the model variables
407  }
408  }
409 
410  void standard_solve(model &md, gmm::iteration &iter,
411  rmodel_plsolver_type lsolver,
412  abstract_newton_line_search &ls) {
413  standard_solve<rmodel_plsolver_type>(md, iter, lsolver, ls);
414  }
415 
416  void standard_solve(model &md, gmm::iteration &iter,
417  cmodel_plsolver_type lsolver,
418  abstract_newton_line_search &ls) {
419  standard_solve<cmodel_plsolver_type>(md, iter, lsolver, ls);
420  }
421 
422 
423  void standard_solve(model &md, gmm::iteration &iter,
424  rmodel_plsolver_type lsolver) {
425  default_newton_line_search ls;
426  standard_solve(md, iter, lsolver, ls);
427  }
428 
429  void standard_solve(model &md, gmm::iteration &iter,
430  cmodel_plsolver_type lsolver) {
431  newton_search_with_step_control ls;
432  // default_newton_line_search ls;
433  standard_solve(md, iter, lsolver, ls);
434  }
435 
436  void standard_solve(model &md, gmm::iteration &iter) {
437  newton_search_with_step_control ls;
438  // default_newton_line_search ls;
439  if (md.is_complex())
440  standard_solve(md, iter, cdefault_linear_solver(md), ls);
441  else
442  standard_solve(md, iter, rdefault_linear_solver(md), ls);
443  }
444 
445 
446 
447 } /* end of namespace getfem. */
448 
getfem_model_solvers.h
Standard solvers for model bricks.
gmm::resize
void resize(M &v, size_type m, size_type n)
*‍/
Definition: gmm_blas.h:231
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
getfem::model
`‘Model’' variables store the variables, the data and the description of a model.
Definition: getfem_models.h:114
gmm::iteration
The Iteration object calculates whether the solution has reached the desired accuracy,...
Definition: gmm_iter.h:53
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:46
gmm::vect_norm2
number_traits< typename linalg_traits< V >::value_type >::magnitude_type vect_norm2(const V &v)
Euclidean norm of a vector.
Definition: gmm_blas.h:557
gmm_inoutput.h
Input/output on sparse matrices.
gmm::fill_random
void fill_random(L &l, double cfill)
*‍/
Definition: gmm_blas.h:154
gmm::vect_norm1
number_traits< typename linalg_traits< V >::value_type >::magnitude_type vect_norm1(const V &v)
1-norm of a vector
Definition: gmm_blas.h:646
bgeot::alpha
size_type alpha(short_type n, short_type d)
Return the value of which is the number of monomials of a polynomial of variables and degree .
Definition: bgeot_poly.cc:47
getfem::standard_solve
void standard_solve(model &md, gmm::iteration &iter, rmodel_plsolver_type lsolver, abstract_newton_line_search &ls)
A default solver for the model brick system.
Definition: getfem_model_solvers.cc:410