ug4
newton_impl.h
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2010-2015: G-CSC, Goethe University Frankfurt
3  * Author: Andreas Vogel
4  *
5  * This file is part of UG4.
6  *
7  * UG4 is free software: you can redistribute it and/or modify it under the
8  * terms of the GNU Lesser General Public License version 3 (as published by the
9  * Free Software Foundation) with the following additional attribution
10  * requirements (according to LGPL/GPL v3 §7):
11  *
12  * (1) The following notice must be displayed in the Appropriate Legal Notices
13  * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14  *
15  * (2) The following notice must be displayed at a prominent place in the
16  * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17  *
18  * (3) The following bibliography is recommended for citation and must be
19  * preserved in all covered files:
20  * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21  * parallel geometric multigrid solver on hierarchically distributed grids.
22  * Computing and visualization in science 16, 4 (2013), 151-164"
23  * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24  * flexible software system for simulating pde based models on high performance
25  * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26  *
27  * This program is distributed in the hope that it will be useful,
28  * but WITHOUT ANY WARRANTY; without even the implied warranty of
29  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30  * GNU Lesser General Public License for more details.
31  */
32 
33 #ifndef __H__UG__LIB_DISC__OPERATOR__NON_LINEAR_OPERATOR__NEWTON_SOLVER__NEWTON_IMPL__
34 #define __H__UG__LIB_DISC__OPERATOR__NON_LINEAR_OPERATOR__NEWTON_SOLVER__NEWTON_IMPL__
35 
36 #include <iostream>
37 #include <sstream>
38 #include <limits>
39 
40 #include "newton.h"
43 
44 #include "lib_disc/operator/non_linear_operator/newton_solver/nestedNewtonRFSwitch.h"
45 
46 #define PROFILE_NEWTON
47 #ifdef PROFILE_NEWTON
48  #define NEWTON_PROFILE_FUNC() PROFILE_FUNC_GROUP("Newton")
49  #define NEWTON_PROFILE_BEGIN(name) PROFILE_BEGIN_GROUP(name, "Newton")
50  #define NEWTON_PROFILE_END() PROFILE_END()
51 #else
52  #define NEWTON_PROFILE_FUNC()
53  #define NEWTON_PROFILE_BEGIN(name)
54  #define NEWTON_PROFILE_END()
55 #endif
56 
57 namespace ug{
58 
59 template <typename TAlgebra>
63  SmartPtr<ILineSearch<vector_type> > spLineSearch) :
64  m_spLinearSolver(LinearSolver),
65  m_spConvCheck(spConvCheck),
66  m_spLineSearch(spLineSearch),
67  m_N(NULL),
68  m_J(NULL),
69  m_spAss(NULL),
70  m_reassembe_J_freq(0),
71  m_dgbCall(0),
72  m_lastNumSteps(0)
73 #if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
74  ,
75  m_newtonUpdater(new NewtonUpdaterGeneric<vector_type>{})
76 #endif
77 {};
78 
79 template <typename TAlgebra>
81 NewtonSolver() :
82  m_spLinearSolver(NULL),
83  m_spConvCheck(new StdConvCheck<vector_type>(10, 1e-8, 1e-10, true)),
84  m_spLineSearch(NULL),
85  m_N(NULL),
86  m_J(NULL),
87  m_spAss(NULL),
88  m_reassembe_J_freq(0),
89  m_dgbCall(0),
90  m_lastNumSteps(0)
91 #if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
92  ,
93  m_newtonUpdater(new NewtonUpdaterGeneric<vector_type>{})
94 #endif
95 {};
96 
97 template <typename TAlgebra>
100  m_spLinearSolver(NULL),
101  m_spConvCheck(new StdConvCheck<vector_type>(10, 1e-8, 1e-10, true)),
102  m_spLineSearch(NULL),
103  m_N(NULL),
104  m_J(NULL),
105  m_spAss(NULL),
106  m_reassembe_J_freq(0),
107  m_dgbCall(0),
108  m_lastNumSteps(0)
109 #if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
110  ,
111  m_newtonUpdater(new NewtonUpdaterGeneric<vector_type>{})
112 #endif
113 {
114  init(N);
115 };
116 
117 template <typename TAlgebra>
120  m_spLinearSolver(NULL),
121  m_spConvCheck(new StdConvCheck<vector_type>(10, 1e-8, 1e-10, true)),
122  m_spLineSearch(NULL),
123  m_N(NULL),
124  m_J(NULL),
125  m_spAss(NULL),
126  m_reassembe_J_freq(0),
127  m_dgbCall(0),
128  m_lastNumSteps(0)
129 #if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
130  ,
131  m_newtonUpdater(new NewtonUpdaterGeneric<vector_type>{})
132 #endif
133 {
134  m_spAss = spAss;
136 };
137 
138 template <typename TAlgebra>
141 {
142  m_spConvCheck = spConvCheck;
143  m_spConvCheck->set_offset(3);
144  m_spConvCheck->set_symbol('#');
145  m_spConvCheck->set_name("Newton Solver");
146 }
147 
148 template <typename TAlgebra>
150 {
151  NEWTON_PROFILE_BEGIN(NewtonSolver_init);
152  m_N = N.template cast_dynamic<AssembledOperator<TAlgebra> >();
153  if(m_N.invalid())
154  UG_THROW("NewtonSolver: currently only works for AssembledDiscreteOperator.");
155 
156  m_spAss = m_N->discretization();
157  return true;
158 }
159 
160 template <typename TAlgebra>
162 {
163 // todo: maybe remove this from interface
164  return true;
165 }
166 
167 template <typename TAlgebra>
169 {
170  NEWTON_PROFILE_BEGIN(NewtonSolver_apply);
171 // increase call count
172  m_dgbCall++;
173 
174 // Check for linear solver
175  if(m_spLinearSolver.invalid())
176  UG_THROW("NewtonSolver::apply: Linear Solver not set.");
177 
178 // Jacobian
179  if(m_J.invalid() || m_J->discretization() != m_spAss) {
180  m_J = make_sp(new AssembledLinearOperator<TAlgebra>(m_spAss));
181  }
182  m_J->set_level(m_N->level());
183 
184 // create tmp vectors
185  SmartPtr<vector_type> spD = u.clone_without_values();
186  SmartPtr<vector_type> spC = u.clone_without_values();
187 
188 // Set dirichlet values
189  try{
190  m_N->prepare(u);
191  }
192  UG_CATCH_THROW("NewtonSolver::prepare: Prepare of Operator failed.");
193 
194 // Compute first Defect
195  try{
196  NEWTON_PROFILE_BEGIN(NewtonComputeDefect1);
197  m_N->apply(*spD, u);
199  }UG_CATCH_THROW("NewtonSolver::apply: Computation of Start-Defect failed.");
200 
201 // loop counts (for the the convergence rate statistics etc.)
202  int loopCnt = 0;
203  m_lastNumSteps = 0;
204 
205 // write start defect for debug
206  char debug_name_ext[20];
207  if (this->debug_writer_valid())
208  {
209  snprintf(debug_name_ext, 20, "_iter%03d", loopCnt);
210  write_debug(*spD, std::string("NEWTON_Defect") + debug_name_ext);
211  write_debug(u, "NEWTON_StartSolution");
212  }
213 
214 // increase offset of output for linear solver
215  const int stdLinOffset = m_spLinearSolver->standard_offset();
216  m_spLinearSolver->convergence_check()->set_offset(stdLinOffset + 3);
217 
218 // set info string indicating the used linear solver
219  std::stringstream ss; ss << "(Linear Solver: " << m_spLinearSolver->name() << ")";
220  m_spConvCheck->set_info(ss.str());
221 
222 // start convergence check
223  m_spConvCheck->start(*spD);
224 
225  for(size_t i = 0; i < m_stepUpdate.size(); ++i)
226  m_stepUpdate[i]->update();
227 
228 // loop iteration
229  while(!m_spConvCheck->iteration_ended())
230  {
231  m_lastNumSteps = loopCnt;
232 
233  // set c = 0
234  NEWTON_PROFILE_BEGIN(NewtonSetCorretionZero);
235  spC->set(0.0);
237 
238  for(size_t i = 0; i < m_innerStepUpdate.size(); ++i)
239  m_innerStepUpdate[i]->update();
240 
241  // Compute Jacobian
242  try{
243  if(m_reassembe_J_freq == 0 || loopCnt % m_reassembe_J_freq == 0) // if we need to reassemble
244  {
245  NEWTON_PROFILE_BEGIN(NewtonComputeJacobian);
246  m_J->init(u);
248  }
249  }UG_CATCH_THROW("NewtonSolver::apply: Initialization of Jacobian failed.");
250 
251  // Write the current Jacobian for debug and prepare the section for the lin. solver
252  if (this->debug_writer_valid())
253  {
254  write_debug(m_J->get_matrix(), std::string("NEWTON_Jacobian") + debug_name_ext);
255  this->enter_debug_writer_section(std::string("NEWTON_LinSolver") + debug_name_ext);
256  }
257 
258  // Init Jacobi Inverse
259  try{
260  NEWTON_PROFILE_BEGIN(NewtonPrepareLinSolver);
261  if(!m_spLinearSolver->init(m_J, u))
262  {
263  UG_LOG("ERROR in 'NewtonSolver::apply': Cannot init Inverse Linear "
264  "Operator for Jacobi-Operator.\n");
265  return false;
266  }
268  }UG_CATCH_THROW("NewtonSolver::apply: Initialization of Linear Solver failed.");
269 
270  // Solve Linearized System
271  try{
272  NEWTON_PROFILE_BEGIN(NewtonApplyLinSolver);
273  if(!m_spLinearSolver->apply(*spC, *spD))
274  {
275  UG_LOG("ERROR in 'NewtonSolver::apply': Cannot apply Inverse Linear "
276  "Operator for Jacobi-Operator.\n");
277  return false;
278  }
280  }UG_CATCH_THROW("NewtonSolver::apply: Application of Linear Solver failed.");
281 
282  this->leave_debug_writer_section();
283 
284  // store convergence history
285  const int numSteps = m_spLinearSolver->step();
286  if(loopCnt >= (int)m_vTotalLinSolverSteps.size()) m_vTotalLinSolverSteps.resize(loopCnt+1);
287  if(loopCnt >= (int)m_vLinSolverCalls.size()) m_vLinSolverCalls.resize(loopCnt+1, 0);
288  if(loopCnt >= (int)m_vLinSolverRates.size()) m_vLinSolverRates.resize(loopCnt+1, 0);
289  m_vTotalLinSolverSteps[loopCnt] += numSteps;
290  m_vLinSolverCalls[loopCnt] += 1;
291  m_vLinSolverRates[loopCnt] += m_spLinearSolver->convergence_check()->avg_rate();
292 
293  try{
294  // Line Search
295  if(m_spLineSearch.valid())
296  {
297  m_spLineSearch->set_offset(" # ");
298  NEWTON_PROFILE_BEGIN(NewtonLineSearch);
299 
300 #if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
301 
302  m_spLineSearch->setNewtonUpdater(m_newtonUpdater);
303 
304 #endif
305  if(!m_spLineSearch->search(m_N, u, *spC, *spD, m_spConvCheck->defect()))
306  {
307  UG_LOG("ERROR in 'NewtonSolver::apply': "
308  "Newton Solver did not converge.\n");
309  return false;
310  }
312  }
313  // No line search: Compute new defect
314  else
315  {
316 #if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
317  // update solution
318  //u -= *spC;
319 // if( ! (*m_newtonUpdater)() )
320 
321  //bool acceptedNewtonUpdate = m_newtonUpdater->updateNewton(u,*spC);
322 
323  //if( ! acceptedNewtonUpdate )
324 
325  if( ! m_newtonUpdater->updateSolution(u,*spC) )
326  {
327  UG_LOG("ERROR in 'NewtonSolver::apply': "
328  "Newton Update did not work.\n");
329  // TODO FIXME was macht conv check update? wie kriege ich hier einfach ein riesiges Residuum rein, ohne was zu rechnen?
330  return false;
331  }
332 
333  if( ! m_newtonUpdater->tellAndFixUpdateEvents(u) )
334  {
335  UG_LOG("unable to fix local Newton updates" << std::endl );
336  return false;
337  }
338 
339 
340 #else
341  // update solution
342  u -= *spC;
343 #endif
344 
345  // compute new Defect
346  NEWTON_PROFILE_BEGIN(NewtonComputeDefect);
347  m_N->prepare(u);
348  m_N->apply(*spD, u);
350  }
351  }UG_CATCH_THROW("NewtonSolver::apply: Line Search update failed.");
352 
353 
354  // update counter
355  loopCnt++;
356 
357  // check convergence
358  m_spConvCheck->update(*spD);
359  if(loopCnt-1 >= (int)m_vNonLinSolverRates.size()) m_vNonLinSolverRates.resize(loopCnt, 0);
360  m_vNonLinSolverRates[loopCnt-1] += m_spConvCheck->rate();
361 
362  // write defect for debug
363  if (this->debug_writer_valid())
364  {
365  snprintf(debug_name_ext, 20, "_iter%03d", loopCnt);
366  write_debug(*spD, std::string("NEWTON_Defect") + debug_name_ext);
367  write_debug(*spC, std::string("NEWTON_Correction") + debug_name_ext);
368  write_debug(u, std::string("NEWTON_Solution") + debug_name_ext);
369  }
370  }
371 
372  // reset offset of output for linear solver to previous value
373  m_spLinearSolver->convergence_check()->set_offset(stdLinOffset);
374 
375  return m_spConvCheck->post();
376 }
377 
378 template <typename TAlgebra>
380 {
381  UG_LOG("\nNewton solver convergence history:\n");
382  UG_LOG("Newton Step | Num Calls | Total Lin Iters | Avg Lin Iters | Avg Nonlin Rates | Avg Lin Rates \n");
383  int allCalls = 0, allLinSteps = 0;
384  number allLinRatesProduct = 1.0, allNonLinRatesProduct = 1.0;
385  for(int call = 0; call < (int)m_vLinSolverCalls.size(); ++call)
386  {
387  UG_LOG( " " << std::setw(10) << call+1 << " | ");
388  UG_LOG(std::setw(9) << m_vLinSolverCalls[call] << " | ");
389  allCalls += m_vLinSolverCalls[call];
390  UG_LOG(std::setw(15) << m_vTotalLinSolverSteps[call] << " | ");
391  allLinSteps += m_vTotalLinSolverSteps[call];
392  UG_LOG(std::setw(13) << std::setprecision(2) << std::fixed << m_vTotalLinSolverSteps[call] / (double)m_vLinSolverCalls[call] << " | ");
393  allNonLinRatesProduct *= pow((number)m_vNonLinSolverRates[call]/(double)m_vLinSolverCalls[call],(double)m_vLinSolverCalls[call]);
394  UG_LOG(std::setw(16) << std::setprecision(6) << std::scientific << m_vNonLinSolverRates[call] / (double)m_vLinSolverCalls[call] << " | ");
395  allLinRatesProduct *= (number)std::pow((number)m_vLinSolverRates[call]/(double)m_vLinSolverCalls[call],(number)m_vTotalLinSolverSteps[call]);
396  UG_LOG(std::setw(13) << std::setprecision(6) << std::scientific << m_vLinSolverRates[call] / (double)m_vLinSolverCalls[call]);
397  UG_LOG("\n");
398  }
399  UG_LOG( " all | ");
400  UG_LOG(std::setw(9) << allCalls << " | ");
401  UG_LOG(std::setw(15) << allLinSteps << " | ");
402  UG_LOG(std::setw(13) << std::setprecision(2) << std::fixed << allLinSteps / (number)allCalls << " | ");
403  UG_LOG(std::setw(16) << std::setprecision(6) << std::scientific << std::pow((number)allNonLinRatesProduct,(number)1.0/(number)allCalls) << " | ");
404  UG_LOG(std::setw(13) << std::setprecision(6) << std::scientific << std::pow((number)allLinRatesProduct,(number)1.0/(number)allLinSteps));
405  UG_LOG("\n");
406 }
407 
408 template <typename TAlgebra>
410 {
411  return m_vLinSolverCalls.size();
412 }
413 
414 template <typename TAlgebra>
416 {
417  return m_vLinSolverCalls[call];
418 }
419 
420 template <typename TAlgebra>
422 {
423  return m_vTotalLinSolverSteps[call];
424 }
425 
426 template <typename TAlgebra>
428 {
429  return m_vTotalLinSolverSteps[call] / ((double)m_vLinSolverCalls[call]);
430 }
431 
432 template <typename TAlgebra>
434 {
435  int allCalls = 0;
436  for(size_t call = 0; call < m_vLinSolverCalls.size(); ++call)
437  allCalls += m_vLinSolverCalls[call];
438  return allCalls;
439 }
440 
441 template <typename TAlgebra>
443 {
444  int allSteps = 0;
445  for(size_t call = 0; call < m_vLinSolverCalls.size(); ++call)
446  allSteps += m_vTotalLinSolverSteps[call];
447  return allSteps;
448 }
449 
450 template <typename TAlgebra>
452 {
453  return total_linsolver_steps()/((double)total_linsolver_calls());
454 }
455 
456 
457 template <typename TAlgebra>
459 {
460  m_vLinSolverRates.clear();
461  m_vNonLinSolverRates.clear();
462  m_vLinSolverCalls.clear();
463  m_vTotalLinSolverSteps.clear();
464 }
465 
466 template <typename TAlgebra>
468 {
469 // add call count to name
470  char ext[20]; snprintf(ext, 20, "_call%03d", m_dgbCall);
471 
472 // write
473  typedef DebugWritingObject<TAlgebra> base_writer_type;
474  base_writer_type::write_debug(vec, name + ext);
475 }
476 
477 template <typename TAlgebra>
479 {
480 // add call count to name
481  char ext[20]; snprintf(ext, 20, "_call%03d", m_dgbCall);
482 
483 // write
484  typedef DebugWritingObject<TAlgebra> base_writer_type;
485  base_writer_type::write_debug(mat, name + ext);
486 }
487 
488 template <typename TAlgebra>
490 {
491  std::stringstream ss;
492  ss << "NewtonSolver\n";
493  ss << " LinearSolver: ";
494  if(m_spLinearSolver.valid()) ss << ConfigShift(m_spLinearSolver->config_string()) << "\n";
495  else ss << " NOT SET!\n";
496  ss << " ConvergenceCheck: ";
497  if(m_spConvCheck.valid()) ss << ConfigShift(m_spConvCheck->config_string()) << "\n";
498  else ss << " NOT SET!\n";
499  ss << " LineSearch: ";
500  if(m_spLineSearch.valid()) ss << ConfigShift(m_spLineSearch->config_string()) << "\n";
501  else ss << " not set.\n";
502  if(m_reassembe_J_freq != 0) ss << " Reassembling Jacobian only once per " << m_reassembe_J_freq << " step(s)\n";
503  return ss.str();
504 }
505 
506 
507 }
508 
509 #endif /* __H__UG__LIB_DISC__OPERATOR__NON_LINEAR_OPERATOR__NEWTON_SOLVER__NEWTON_IMPL__ */
510 
location name
Definition: checkpoint_util.lua:128
Definition: smart_pointer.h:108
matrix operator based on the assembling of a problem
Definition: assembled_linear_operator.h:60
Definition: assembled_non_linear_operator.h:43
Definition: debug_writer.h:354
TAlgebra::vector_type vector_type
type of vector
Definition: debug_writer.h:360
TAlgebra::matrix_type matrix_type
type of matrix
Definition: debug_writer.h:363
Interface providing Jacobian and Defect of a discretization.
Definition: assemble_interface.h:110
Definition: convergence_check.h:72
Definition: line_search.h:71
describes an inverse linear mapping X->Y
Definition: linear_operator_inverse.h:80
describes a mapping X->Y
Definition: operator.h:86
linear solver using abstract preconditioner interface
Definition: linear_solver.h:57
int total_linsolver_calls() const
Definition: newton_impl.h:433
void clear_average_convergence()
resets average linear solver convergence
Definition: newton_impl.h:458
virtual bool prepare(vector_type &u)
prepare Operator
Definition: newton_impl.h:161
virtual bool apply(vector_type &u)
apply Operator, i.e. N^{-1}(0) = u
Definition: newton_impl.h:168
NewtonSolver()
default constructor
Definition: newton_impl.h:81
virtual bool init(SmartPtr< IOperator< vector_type > > N)
This operator inverts the Operator N: Y -> X.
Definition: newton_impl.h:149
int num_linsolver_steps(size_t call) const
Definition: newton_impl.h:421
void write_debug(const vector_type &vec, std::string filename)
Definition: newton_impl.h:467
double total_average_linear_steps() const
Definition: newton_impl.h:451
virtual std::string config_string() const
returns information about configuration parameters
Definition: newton_impl.h:489
double average_linear_steps(size_t call) const
Definition: newton_impl.h:427
void print_average_convergence() const
prints average linear solver convergence
Definition: newton_impl.h:379
size_t num_newton_steps() const
Definition: newton_impl.h:409
void set_convergence_check(SmartPtr< IConvergenceCheck< vector_type > > spConvCheck)
sets the convergence check
Definition: newton_impl.h:140
int num_linsolver_calls(size_t call) const
Definition: newton_impl.h:415
int total_linsolver_steps() const
Definition: newton_impl.h:442
Definition: newtonUpdaterGeneric.h:25
Definition: convergence_check.h:179
virtual void init()
#define UG_CATCH_THROW(msg)
Definition: error.h:64
#define UG_THROW(msg)
Definition: error.h:57
#define UG_LOG(msg)
Definition: log.h:367
double number
Definition: types.h:124
if(!(yy_init))
Definition: lexer.cpp:997
the ug namespace
string ConfigShift(string s)
Definition: string_util.cpp:457
#define NEWTON_PROFILE_END()
Definition: newton_impl.h:50
#define NEWTON_PROFILE_BEGIN(name)
Definition: newton_impl.h:49
SmartPtr< T, FreePolicy > make_sp(T *inst)
returns a SmartPtr for the passed raw pointer
Definition: smart_pointer.h:836
function ProblemDisc new(problemDesc, dom)