Hermes2D  2.0
newton_solver.cpp
1 // This file is part of Hermes2D
2 //
3 // Copyright (c) 2009 hp-FEM group at the University of Nevada, Reno (UNR).
4 // Email: hpfem-group@unr.edu, home page: http://hpfem.org/.
5 //
6 // Hermes2D is free software; you can redistribute it and/or modify
7 // it under the terms of the GNU General Public License as published
8 // by the Free Software Foundation; either version 2 of the License,
9 // or (at your option) any later version.
10 //
11 // Hermes2D is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 // GNU General Public License for more details.
15 //
16 // You should have received a copy of the GNU General Public License
17 // along with Hermes2D; if not, write to the Free Software
18 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
19 #include "newton_solver.h"
20 #include "projections/ogprojection.h"
21 #include "hermes_common.h"
22 
23 namespace Hermes
24 {
25  namespace Hermes2D
26  {
27  template<typename Scalar>
28  NewtonSolver<Scalar>::NewtonSolver() : NonlinearSolver<Scalar>(new DiscreteProblem<Scalar>()), own_dp(true), kept_jacobian(NULL)
29  {
30  init_attributes();
31  init_linear_solver();
32  }
33 
34  template<typename Scalar>
35  NewtonSolver<Scalar>::NewtonSolver(DiscreteProblem<Scalar>* dp) : NonlinearSolver<Scalar>(dp), own_dp(false), kept_jacobian(NULL)
36  {
37  init_attributes();
38  init_linear_solver();
39  }
40 
41  template<typename Scalar>
42  NewtonSolver<Scalar>::NewtonSolver(const WeakForm<Scalar>* wf, const Space<Scalar>* space) : NonlinearSolver<Scalar>(new DiscreteProblem<Scalar>(wf, space)), own_dp(true), kept_jacobian(NULL)
43  {
44  init_attributes();
45  init_linear_solver();
46  }
47 
48  template<typename Scalar>
49  NewtonSolver<Scalar>::NewtonSolver(const WeakForm<Scalar>* wf, Hermes::vector<const Space<Scalar> *> spaces) : NonlinearSolver<Scalar>(new DiscreteProblem<Scalar>(wf, spaces)), own_dp(true), kept_jacobian(NULL)
50  {
51  init_attributes();
52  init_linear_solver();
53  }
54 
55  template<typename Scalar>
56  bool NewtonSolver<Scalar>::isOkay() const
57  {
58  if(static_cast<DiscreteProblem<Scalar>*>(this->dp)->get_weak_formulation() == NULL)
59  return false;
60  if(static_cast<DiscreteProblem<Scalar>*>(this->dp)->get_spaces().size() == 0)
61  return false;
62  return true;
63  }
64 
65  template<typename Scalar>
66  void NewtonSolver<Scalar>::init_attributes()
67  {
68  this->newton_tol = 1e-8;
69  this->newton_max_iter = 15;
70  this->residual_as_function = false;
71  this->max_allowed_residual_norm = 1E9;
72  this->min_allowed_damping_coeff = 1E-4;
73  this->currentDampingCofficient = 1.0;
74  this->manual_damping = false;
75  this->auto_damping_ratio = 2.0;
76  this->initial_auto_damping_ratio = 1.0;
77  this->sufficient_improvement_factor = 0.95;
78  this->necessary_successful_steps_to_increase = 1;
79  }
80 
81  template<typename Scalar>
82  void NewtonSolver<Scalar>::set_newton_tol(double newton_tol)
83  {
84  this->newton_tol = newton_tol;
85  }
86 
87  template<typename Scalar>
88  void NewtonSolver<Scalar>::set_weak_formulation(const WeakForm<Scalar>* wf)
89  {
90  (static_cast<DiscreteProblem<Scalar>*>(this->dp))->set_weak_formulation(wf);
91  }
92 
93  template<typename Scalar>
94  void NewtonSolver<Scalar>::set_newton_max_iter(int newton_max_iter)
95  {
96  if(newton_max_iter < 1)
97  throw Exceptions::ValueException("newton_max_iter", newton_max_iter, 1);
98  this->newton_max_iter = newton_max_iter;
99  }
100 
101  template<typename Scalar>
102  void NewtonSolver<Scalar>::set_max_allowed_residual_norm(double max_allowed_residual_norm_to_set)
103  {
104  if(max_allowed_residual_norm_to_set <= 0.0)
105  throw Exceptions::ValueException("max_allowed_residual_norm_to_set", max_allowed_residual_norm_to_set, 0.0);
106  this->max_allowed_residual_norm = max_allowed_residual_norm_to_set;
107  }
108 
109  template<typename Scalar>
110  void NewtonSolver<Scalar>::set_min_allowed_damping_coeff(double min_allowed_damping_coeff_to_set)
111  {
112  if(min_allowed_damping_coeff_to_set <= 0.0)
113  throw Exceptions::ValueException("min_allowed_damping_coeff_to_set", min_allowed_damping_coeff_to_set, 0.0);
114  this->min_allowed_damping_coeff = min_allowed_damping_coeff_to_set;
115  }
116 
117  template<typename Scalar>
118  void NewtonSolver<Scalar>::set_residual_as_function()
119  {
120  this->residual_as_function = true;
121  }
122 
123  template<typename Scalar>
124  void NewtonSolver<Scalar>::set_time(double time)
125  {
126  Hermes::vector<Space<Scalar>*> spaces;
127  for(unsigned int i = 0; i < static_cast<DiscreteProblem<Scalar>*>(this->dp)->get_spaces().size(); i++)
128  spaces.push_back(const_cast<Space<Scalar>*>(static_cast<DiscreteProblem<Scalar>*>(this->dp)->get_space(i)));
129 
130  Space<Scalar>::update_essential_bc_values(spaces, time);
131  const_cast<WeakForm<Scalar>*>(static_cast<DiscreteProblem<Scalar>*>(this->dp)->wf)->set_current_time(time);
132  }
133 
134  template<typename Scalar>
135  void NewtonSolver<Scalar>::set_time_step(double time_step)
136  {
137  const_cast<WeakForm<Scalar>*>(static_cast<DiscreteProblem<Scalar>*>(this->dp)->wf)->set_current_time_step(time_step);
138  }
139 
140  template<typename Scalar>
141  void NewtonSolver<Scalar>::set_spaces(Hermes::vector<const Space<Scalar>*> spaces)
142  {
143  static_cast<DiscreteProblem<Scalar>*>(this->dp)->set_spaces(spaces);
144  if(kept_jacobian != NULL)
145  delete kept_jacobian;
146  kept_jacobian = NULL;
147  }
148 
149  template<typename Scalar>
150  void NewtonSolver<Scalar>::set_space(const Space<Scalar>* space)
151  {
152  static_cast<DiscreteProblem<Scalar>*>(this->dp)->set_space(space);
153  if(kept_jacobian != NULL)
154  delete kept_jacobian;
155  kept_jacobian = NULL;
156  }
157 
158  template<typename Scalar>
159  Hermes::vector<const Space<Scalar>*> NewtonSolver<Scalar>::get_spaces() const
160  {
161  return static_cast<DiscreteProblem<Scalar>*>(this->dp)->get_spaces();
162  }
163 
164  template<typename Scalar>
165  void NewtonSolver<Scalar>::init_linear_solver()
166  {
167  // Set up the solver, jacobian, and residual according to the solver selection.
168  jacobian = create_matrix<Scalar>();
169  residual = create_vector<Scalar>();
170  linear_solver = create_linear_solver<Scalar>(jacobian, residual);
171  }
172 
173  template<typename Scalar>
174  NewtonSolver<Scalar>::~NewtonSolver()
175  {
176  if(kept_jacobian != NULL)
177  delete kept_jacobian;
178  delete jacobian;
179  delete residual;
180  delete linear_solver;
181  if(own_dp)
182  delete this->dp;
183  else
184  static_cast<DiscreteProblem<Scalar>*>(this->dp)->have_matrix = false;
185  }
186 
187  template<typename Scalar>
188  void NewtonSolver<Scalar>::set_manual_damping_coeff(bool onOff, double coeff)
189  {
190  if(coeff <= 0.0 || coeff > 1.0)
191  throw Exceptions::ValueException("coeff", coeff, 0.0, 1.0);
192  if(onOff)
193  {
194  this->manual_damping = true;
195  this->currentDampingCofficient = coeff;
196  }
197  else
198  this->manual_damping = false;
199  }
200 
201  template<typename Scalar>
202  void NewtonSolver<Scalar>::set_initial_auto_damping_coeff(double coeff)
203  {
204  if(coeff <= 0.0 || coeff > 1.0)
205  throw Exceptions::ValueException("coeff", coeff, 0.0, 1.0);
206  if(this->manual_damping)
207  this->warn("Manual damping is turned on and you called set_initial_auto_damping_coeff(), turn off manual damping first by set_manual_damping_coeff(false);");
208  else
209  this->currentDampingCofficient = coeff;
210  }
211 
212  template<typename Scalar>
213  void NewtonSolver<Scalar>::set_auto_damping_ratio(double ratio)
214  {
215  if(ratio <= 0.0)
216  throw Exceptions::ValueException("ratio", ratio, 0.0);
217  if(this->manual_damping)
218  this->warn("Manual damping is turned on and you called set_initial_auto_damping_coeff(), turn off manual damping first by set_manual_damping_coeff(false);");
219  this->auto_damping_ratio = ratio;
220  }
221 
222  template<typename Scalar>
223  void NewtonSolver<Scalar>::set_sufficient_improvement_factor(double ratio)
224  {
225  if(ratio <= 0.0)
226  throw Exceptions::ValueException("ratio", ratio, 0.0);
227  if(this->manual_damping)
228  this->warn("Manual damping is turned on and you called set_initial_auto_damping_coeff(), turn off manual damping first by set_manual_damping_coeff(false);");
229  this->sufficient_improvement_factor = ratio;
230  }
231 
232  template<typename Scalar>
233  void NewtonSolver<Scalar>::set_necessary_successful_steps_to_increase(unsigned int steps)
234  {
235  if(steps < 1)
236  throw Exceptions::ValueException("steps", steps, 0.0);
237  if(this->manual_damping)
238  this->warn("Manual damping is turned on and you called set_initial_auto_damping_coeff(), turn off manual damping first by set_manual_damping_coeff(false);");
239  this->necessary_successful_steps_to_increase = steps;
240  }
241 
242  template<typename Scalar>
243  void NewtonSolver<Scalar>::solve(Solution<Scalar>* initial_guess)
244  {
245  Hermes::vector<Solution<Scalar>*> vectorToPass;
246  vectorToPass.push_back(initial_guess);
247  this->solve(vectorToPass);
248  }
249 
250  template<typename Scalar>
251  void NewtonSolver<Scalar>::solve(Hermes::vector<Solution<Scalar>*> initial_guess)
252  {
253  int ndof = this->dp->get_num_dofs();
254  Scalar* coeff_vec = new Scalar[ndof];
255  OGProjection<Scalar> ogProjection;
256  ogProjection.project_global(static_cast<DiscreteProblem<Scalar>*>(this->dp)->get_spaces(), initial_guess, coeff_vec);
257  this->solve(coeff_vec);
258  }
259 
260  template<typename Scalar>
261  void NewtonSolver<Scalar>::solve(Scalar* coeff_vec)
262  {
263  this->check();
264 
265  this->tick();
266 
267  // Obtain the number of degrees of freedom.
268  int ndof = this->dp->get_num_dofs();
269 
270  // If coeff_vec == NULL then create a zero vector.
271  bool delete_coeff_vec = false;
272  if(coeff_vec == NULL)
273  {
274  coeff_vec = new Scalar[ndof];
275  memset(coeff_vec, 0, ndof*sizeof(Scalar));
276  delete_coeff_vec = true;
277  }
278 
279  // Vector for restarting steps.
280  Scalar* coeff_vec_back = new Scalar[ndof];
281  memset(coeff_vec_back, 0, ndof*sizeof(Scalar));
282 
283  // Delete the old solution vector, if there is any.
284  if(this->sln_vector != NULL)
285  {
286  delete [] this->sln_vector;
287  this->sln_vector = NULL;
288  }
289 
290  // The Newton's loop.
291  double residual_norm;
292  double last_residual_norm;
293  int it = 1;
294  int successfulSteps = 0;
295 
296  this->on_initialization();
297 
298  while (true)
299  {
300  this->on_step_begin();
301 
302  // Assemble just the residual vector.
303  this->dp->assemble(coeff_vec, residual);
304  if(this->output_rhsOn && (this->output_rhsIterations == -1 || this->output_rhsIterations >= it))
305  {
306  char* fileName = new char[this->RhsFilename.length() + 5];
307  if(this->RhsFormat == Hermes::Algebra::DF_MATLAB_SPARSE)
308  sprintf(fileName, "%s%i.m", this->RhsFilename.c_str(), it);
309  else
310  sprintf(fileName, "%s%i", this->RhsFilename.c_str(), it);
311  FILE* rhs_file = fopen(fileName, "w+");
312  residual->dump(rhs_file, this->RhsVarname.c_str(), this->RhsFormat, this->rhs_number_format);
313  fclose(rhs_file);
314  delete [] fileName;
315  }
316 
317  Element* e;
318  for(unsigned int i = 0; i < static_cast<DiscreteProblem<Scalar>*>(this->dp)->get_spaces().size(); i++)
319  for_all_active_elements(e, static_cast<DiscreteProblem<Scalar>*>(this->dp)->get_spaces()[i]->get_mesh())
320  static_cast<DiscreteProblem<Scalar>*>(this->dp)->get_spaces()[i]->edata[e->id].changed_in_last_adaptation = false;
321 
322  // Measure the residual norm.
323  if(residual_as_function)
324  {
325  // Prepare solutions for measuring residual norm.
326  Hermes::vector<Solution<Scalar>*> solutions;
327  Hermes::vector<bool> dir_lift_false;
328  for (unsigned int i = 0; i < static_cast<DiscreteProblem<Scalar>*>(this->dp)->get_spaces().size(); i++)
329  {
330  solutions.push_back(new Solution<Scalar>());
331  dir_lift_false.push_back(false);
332  }
333 
334  Solution<Scalar>::vector_to_solutions(residual, static_cast<DiscreteProblem<Scalar>*>(this->dp)->get_spaces(), solutions, dir_lift_false);
335 
336  // Calculate the norm.
337  residual_norm = Global<Scalar>::calc_norms(solutions);
338 
339  // Clean up.
340  for (unsigned int i = 0; i < static_cast<DiscreteProblem<Scalar>*>(this->dp)->get_spaces().size(); i++)
341  delete solutions[i];
342  }
343  else
344  {
345  // Calculate the l2-norm of residual vector, this is the traditional way.
346  if(it == 1)
347  last_residual_norm = residual_norm = Global<Scalar>::get_l2_norm(residual);
348  else
349  {
350  last_residual_norm = residual_norm;
351  residual_norm = Global<Scalar>::get_l2_norm(residual);
352  }
353  }
354 
355  // Info for the user.
356  if(it == 1)
357  this->info("\tNewton: initial residual norm: %g", residual_norm);
358  else
359  {
360  this->info("\tNewton: iteration %d, residual norm: %g", it - 1, residual_norm);
361  if(!this->manual_damping && !((residual_norm > max_allowed_residual_norm) || (residual_norm < newton_tol && it > 1)))
362  {
363  if(residual_norm < last_residual_norm * this->sufficient_improvement_factor)
364  {
365  if(++successfulSteps >= this->necessary_successful_steps_to_increase)
366  this->currentDampingCofficient = std::min(this->initial_auto_damping_ratio, 2 * this->currentDampingCofficient);
367  if(residual_norm < last_residual_norm)
368  this->info("\t Newton: step successful, calculation continues with damping coefficient: %g.", this->currentDampingCofficient);
369  }
370  else
371  {
372  successfulSteps = 0;
373  if(this->currentDampingCofficient < this->min_allowed_damping_coeff)
374  {
375  this->warn("\t Newton: results NOT improved, current damping coefficient is at the minimum possible level: %g.", min_allowed_damping_coeff);
376  this->info("\t If you want to decrease the minimum level, use the method set_min_allowed_damping_coeff()");
377  if(delete_coeff_vec)
378  {
379  delete [] coeff_vec;
380  coeff_vec = NULL;
381  }
382 
383  delete [] coeff_vec_back;
384 
385  throw Exceptions::Exception("Newton NOT converged because of damping coefficient could not be decreased anymore to possibly handle non-converging process.");
386  }
387  else
388  {
389  this->currentDampingCofficient = (1 / this->auto_damping_ratio) * this->currentDampingCofficient;
390  this->warn(" Newton: results NOT improved, step restarted with damping coefficient: %g.", this->currentDampingCofficient);
391  }
392  }
393  }
394  }
395 
396  // If maximum allowed residual norm is exceeded, fail.
397  if(residual_norm > max_allowed_residual_norm)
398  {
399  this->tick();
400  this->info("\tNewton: solution duration: %f s.\n", this->last());
401  this->on_finish();
402 
403  if(delete_coeff_vec)
404  {
405  delete [] coeff_vec;
406  coeff_vec = NULL;
407  }
408 
409  delete [] coeff_vec_back;
410 
411  throw Exceptions::ValueException("residual norm", residual_norm, max_allowed_residual_norm);
412  }
413 
414  // If residual norm is within tolerance, return 'true'.
415  // This is the only correct way of ending.
416  if(residual_norm < newton_tol && it > 1)
417  {
418  // We want to return the solution in a different structure.
419  this->sln_vector = new Scalar[ndof];
420  for (int i = 0; i < ndof; i++)
421  this->sln_vector[i] = coeff_vec[i];
422 
423  if(delete_coeff_vec)
424  {
425  delete [] coeff_vec;
426  coeff_vec = NULL;
427  }
428 
429  delete [] coeff_vec_back;
430 
431  this->on_finish();
432 
433  this->tick();
434  this->info("\tNewton: solution duration: %f s.\n", this->last());
435 
436  return;
437  }
438 
439  // Assemble just the jacobian.
440  this->dp->assemble(coeff_vec, jacobian);
441  if(this->output_matrixOn && (this->output_matrixIterations == -1 || this->output_matrixIterations >= it))
442  {
443  char* fileName = new char[this->matrixFilename.length() + 5];
444  if(this->matrixFormat == Hermes::Algebra::DF_MATLAB_SPARSE)
445  sprintf(fileName, "%s%i.m", this->matrixFilename.c_str(), it);
446  else
447  sprintf(fileName, "%s%i", this->matrixFilename.c_str(), it);
448  FILE* matrix_file = fopen(fileName, "w+");
449 
450  jacobian->dump(matrix_file, this->matrixVarname.c_str(), this->matrixFormat, this->matrix_number_format);
451  fclose(matrix_file);
452  delete [] fileName;
453  }
454 
455  this->on_step_end();
456 
457  // Multiply the residual vector with -1 since the matrix
458  // equation reads J(Y^n) \deltaY^{n + 1} = -F(Y^n).
459  residual->change_sign();
460 
461  // Solve the linear system.
462  if(!linear_solver->solve())
463  throw Exceptions::LinearMatrixSolverException();
464 
465  // Add \deltaY^{n + 1} to Y^n.
466  // The good case.
467  if(residual_norm < last_residual_norm * this->sufficient_improvement_factor || this->manual_damping || it == 1)
468  {
469  memcpy(coeff_vec_back, coeff_vec, sizeof(Scalar)*ndof);
470  for (int i = 0; i < ndof; i++)
471  coeff_vec[i] += currentDampingCofficient * linear_solver->get_sln_vector()[i];
472  }
473  else
474  {
475  for (int i = 0; i < ndof; i++)
476  coeff_vec[i] = coeff_vec_back[i] + currentDampingCofficient * (coeff_vec[i] - coeff_vec_back[i]);
477  }
478 
479  // Increase the number of iterations and test if we are still under the limit.
480  if(it++ >= newton_max_iter)
481  {
482  // We want to return the solution in a different structure.
483  this->sln_vector = new Scalar[ndof];
484  for (int i = 0; i < ndof; i++)
485  this->sln_vector[i] = coeff_vec[i];
486 
487  if(delete_coeff_vec)
488  {
489  delete [] coeff_vec;
490  coeff_vec = NULL;
491  }
492 
493  delete [] coeff_vec_back;
494  coeff_vec_back = NULL;
495 
496  this->tick();
497  this->info("\tNewton: solution duration: %f s.\n", this->last());
498 
499  this->on_finish();
500  throw Exceptions::ValueException("iterations", it, newton_max_iter);
501  }
502  }
503  }
504 
505  template<typename Scalar>
506  void NewtonSolver<Scalar>::solve_keep_jacobian(Solution<Scalar>* initial_guess)
507  {
508  Hermes::vector<Solution<Scalar>*> vectorToPass;
509  vectorToPass.push_back(initial_guess);
510  this->solve_keep_jacobian(vectorToPass);
511  }
512 
513  template<typename Scalar>
514  void NewtonSolver<Scalar>::solve_keep_jacobian(Hermes::vector<Solution<Scalar>*> initial_guess)
515  {
516  int ndof = this->dp->get_num_dofs();
517  Scalar* coeff_vec = new Scalar[ndof];
518  OGProjection<Scalar> ogProjection;
519  ogProjection.project_global(static_cast<DiscreteProblem<Scalar>*>(this->dp)->get_spaces(), initial_guess, coeff_vec);
520  this->solve_keep_jacobian(coeff_vec);
521  }
522 
523  template<typename Scalar>
524  void NewtonSolver<Scalar>::solve_keep_jacobian(Scalar* coeff_vec)
525  {
526  this->check();
527 
528  this->tick();
529 
530  // Obtain the number of degrees of freedom.
531  int ndof = this->dp->get_num_dofs();
532 
533  // If coeff_vec == NULL then create a zero vector.
534  bool delete_coeff_vec = false;
535  if(coeff_vec == NULL)
536  {
537  coeff_vec = new Scalar[ndof];
538  memset(coeff_vec, 0, ndof*sizeof(Scalar));
539  delete_coeff_vec = true;
540  }
541 
542  // Vector for restarting steps.
543  Scalar* coeff_vec_back = new Scalar[ndof];
544  memset(coeff_vec_back, 0, ndof*sizeof(Scalar));
545 
546 
547  // The Newton's loop.
548  double residual_norm;
549  double last_residual_norm;
550  int it = 1;
551  int successfulSteps = 0;
552 
553  this->on_initialization();
554 
555  while (true)
556  {
557  this->on_step_begin();
558 
559  // Assemble the residual vector.
560  this->dp->assemble(coeff_vec, residual);
561  if(this->output_rhsOn && (this->output_rhsIterations == -1 || this->output_rhsIterations >= it))
562  {
563  char* fileName = new char[this->RhsFilename.length() + 5];
564  if(this->RhsFormat == Hermes::Algebra::DF_MATLAB_SPARSE)
565  sprintf(fileName, "%s%i.m", this->RhsFilename.c_str(), it);
566  else
567  sprintf(fileName, "%s%i", this->RhsFilename.c_str(), it);
568  FILE* rhs_file = fopen(fileName, "w+");
569  residual->dump(rhs_file, this->RhsVarname.c_str(), this->RhsFormat, this->rhs_number_format);
570  fclose(rhs_file);
571  }
572 
573  Element* e;
574  for(unsigned int i = 0; i < static_cast<DiscreteProblem<Scalar>*>(this->dp)->get_spaces().size(); i++)
575  for_all_active_elements(e, static_cast<DiscreteProblem<Scalar>*>(this->dp)->get_spaces()[i]->get_mesh())
576  static_cast<DiscreteProblem<Scalar>*>(this->dp)->get_spaces()[i]->edata[e->id].changed_in_last_adaptation = false;
577 
578  // Measure the residual norm.
579  if(residual_as_function)
580  {
581  // Prepare solutions for measuring residual norm.
582  Hermes::vector<Solution<Scalar>* > solutions;
583  Hermes::vector<bool> dir_lift_false;
584  for (unsigned int i = 0; i < static_cast<DiscreteProblem<Scalar>*>(this->dp)->get_spaces().size(); i++) {
585  solutions.push_back(new Solution<Scalar>());
586  dir_lift_false.push_back(false);
587  }
588  Solution<Scalar>::vector_to_solutions(residual,
589  static_cast<DiscreteProblem<Scalar>*>(this->dp)->get_spaces(), solutions, dir_lift_false);
590 
591  // Calculate the norm.
592  residual_norm = Global<Scalar>::calc_norms(solutions);
593 
594  // Clean up.
595  for (unsigned int i = 0; i < static_cast<DiscreteProblem<Scalar>*>(this->dp)->get_spaces().size(); i++)
596  delete solutions[i];
597  }
598  else
599  {
600  // Calculate the l2-norm of residual vector, this is the traditional way.
601  if(it == 1)
602  last_residual_norm = residual_norm = Global<Scalar>::get_l2_norm(residual);
603  else
604  {
605  last_residual_norm = residual_norm;
606  residual_norm = Global<Scalar>::get_l2_norm(residual);
607  }
608  }
609 
610  // Info for the user.
611  if(it == 1)
612  this->info("\tNewton: initial residual norm: %g", residual_norm);
613  else
614  {
615  this->info("\tNewton: iteration %d, residual norm: %g", it - 1, residual_norm);
616  if(!this->manual_damping && !((residual_norm > max_allowed_residual_norm) || (residual_norm < newton_tol && it > 1)))
617  {
618  if(residual_norm < last_residual_norm * this->sufficient_improvement_factor)
619  {
620  if(++successfulSteps >= this->necessary_successful_steps_to_increase)
621  this->currentDampingCofficient = std::min(this->initial_auto_damping_ratio, 2 * this->currentDampingCofficient);
622  if(residual_norm < last_residual_norm)
623  this->info("\t Newton: step successful, calculation continues with damping coefficient: %g.", this->currentDampingCofficient);
624  }
625  else
626  {
627  successfulSteps = 0;
628  if(this->currentDampingCofficient < this->min_allowed_damping_coeff)
629  {
630  this->warn("\t Newton: results NOT improved, current damping coefficient is at the minimum possible level: %g.", min_allowed_damping_coeff);
631  this->info("\t If you want to decrease the minimum level, use the method set_min_allowed_damping_coeff()");
632 
633  if(delete_coeff_vec)
634  {
635  delete [] coeff_vec;
636  coeff_vec = NULL;
637  }
638 
639  delete [] coeff_vec_back;
640 
641  throw Exceptions::Exception("Newton NOT converged because of damping coefficient could not be decreased anymore to possibly handle non-converging process.");
642  }
643  else
644  {
645  this->currentDampingCofficient = (1 / this->auto_damping_ratio) * this->currentDampingCofficient;
646  this->warn("\t Newton: results NOT improved, step restarted with damping coefficient: %g.", this->currentDampingCofficient);
647  }
648  }
649  }
650  }
651 
652  // If maximum allowed residual norm is exceeded, fail.
653  if(residual_norm > max_allowed_residual_norm)
654  {
655  this->tick();
656  this->info("\tNewton: solution duration: %f s.", this->last());
657 
658  if(delete_coeff_vec)
659  {
660  delete [] coeff_vec;
661  coeff_vec = NULL;
662  }
663 
664  delete [] coeff_vec_back;
665 
666  this->on_finish();
667 
668  throw Exceptions::ValueException("residual norm", residual_norm, max_allowed_residual_norm);
669  }
670 
671  // If residual norm is within tolerance, return 'true'.
672  // This is the only correct way of ending.
673  if(residual_norm < newton_tol && it > 1)
674  {
675  // We want to return the solution in a different structure.
676  this->sln_vector = new Scalar[ndof];
677  for (int i = 0; i < ndof; i++)
678  this->sln_vector[i] = coeff_vec[i];
679 
680  if(delete_coeff_vec)
681  {
682  delete [] coeff_vec;
683  coeff_vec = NULL;
684  }
685 
686  delete [] coeff_vec_back;
687 
688  this->tick();
689  this->info("\tNewton: solution duration: %f s.", this->last());
690 
691  this->on_finish();
692  return;
693  }
694 
695  // Assemble and keep the jacobian if this has not been done before.
696  // Also declare that LU-factorization in case of a direct solver will be done only once and reused afterwards.
697  if(kept_jacobian == NULL || !(static_cast<DiscreteProblem<Scalar>*>(this->dp))->have_matrix)
698  {
699  if(kept_jacobian != NULL)
700  delete kept_jacobian;
701 
702  kept_jacobian = create_matrix<Scalar>();
703 
704  // Give the matrix solver the correct Jacobian. NOTE: It would be cleaner if the whole decision whether to keep
705  // Jacobian or not was made in the constructor.
706  //
707  // Delete the matrix solver created in the constructor.
708  delete linear_solver;
709  // Create new matrix solver with correct matrix.
710  linear_solver = create_linear_solver<Scalar>(kept_jacobian, residual);
711 
712  this->dp->assemble(coeff_vec, kept_jacobian);
713 
714  if(this->output_matrixOn && (this->output_matrixIterations == -1 || this->output_matrixIterations >= it))
715  {
716  char* fileName = new char[this->matrixFilename.length() + 5];
717  if(this->matrixFormat == Hermes::Algebra::DF_MATLAB_SPARSE)
718  sprintf(fileName, "%s%i.m", this->matrixFilename.c_str(), it);
719  else
720  sprintf(fileName, "%s%i", this->matrixFilename.c_str(), it);
721  FILE* matrix_file = fopen(fileName, "w+");
722 
723  kept_jacobian->dump(matrix_file, this->matrixVarname.c_str(), this->matrixFormat, this->matrix_number_format);
724  fclose(matrix_file);
725  }
726 
727  linear_solver->set_factorization_scheme(HERMES_REUSE_FACTORIZATION_COMPLETELY);
728  }
729 
730  this->on_step_end();
731 
732  // Multiply the residual vector with -1 since the matrix
733  // equation reads J(Y^n) \deltaY^{n + 1} = -F(Y^n).
734  residual->change_sign();
735 
736  // Solve the linear system.
737  if(!linear_solver->solve())
738  {
739  throw Exceptions::LinearMatrixSolverException();
740  }
741 
742  // Add \deltaY^{n + 1} to Y^n.
743  // The good case.
744  if(residual_norm < last_residual_norm * this->sufficient_improvement_factor || this->manual_damping || it == 1)
745  {
746  memcpy(coeff_vec_back, coeff_vec, sizeof(Scalar)*ndof);
747  for (int i = 0; i < ndof; i++)
748  coeff_vec[i] += currentDampingCofficient * linear_solver->get_sln_vector()[i];
749  }
750  else
751  {
752  for (int i = 0; i < ndof; i++)
753  coeff_vec[i] = coeff_vec_back[i] + currentDampingCofficient * (coeff_vec[i] - coeff_vec_back[i]);
754  }
755 
756  // Increase the number of iterations and test if we are still under the limit.
757  if(it++ >= newton_max_iter)
758  {
759  // We want to return the solution in a different structure.
760  this->sln_vector = new Scalar[ndof];
761  for (int i = 0; i < ndof; i++)
762  this->sln_vector[i] = coeff_vec[i];
763 
764  if(delete_coeff_vec)
765  {
766  delete [] coeff_vec;
767  coeff_vec = NULL;
768  }
769 
770  delete [] coeff_vec_back;
771 
772  this->tick();
773  this->info("\tNewton: solution duration: %f s.\n", this->last());
774 
775  this->on_finish();
776  throw Exceptions::ValueException("iterations", it, newton_max_iter);
777  }
778  }
779  }
780 
781  template<typename Scalar>
782  void NewtonSolver<Scalar>::set_iterative_method(const char* iterative_method_name)
783  {
784  NonlinearSolver<Scalar>::set_iterative_method(iterative_method_name);
785  // Set iterative method in case of iterative solver AztecOO.
786 #ifdef HAVE_AZTECOO
787  dynamic_cast<Hermes::Solvers::AztecOOSolver<Scalar>*>(linear_solver)->set_solver(iterative_method_name);
788 #else
789  this->warn("Trying to set iterative method without AztecOO present.");
790 #endif
791  }
792 
793  template<typename Scalar>
794  void NewtonSolver<Scalar>::set_preconditioner(const char* preconditioner_name)
795  {
796  NonlinearSolver<Scalar>::set_preconditioner(preconditioner_name);
797  // Set preconditioner in case of iterative solver AztecOO.
798 #ifdef HAVE_AZTECOO
799  dynamic_cast<Hermes::Solvers::AztecOOSolver<Scalar> *>(linear_solver)->set_precond(preconditioner_name);
800 #else
801  this->warn("Trying to set iterative method without AztecOO present.");
802 #endif
803  }
804 
805  template class HERMES_API NewtonSolver<double>;
806  template class HERMES_API NewtonSolver<std::complex<double> >;
807  }
808 }