ug4
Loading...
Searching...
No Matches
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
57namespace ug{
58
59template <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
79template <typename TAlgebra>
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
97template <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
117template <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
138template <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
148template <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
160template <typename TAlgebra>
162{
163// todo: maybe remove this from interface
164 return true;
165}
166
167template <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
378template <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
408template <typename TAlgebra>
410{
411 return m_vLinSolverCalls.size();
412}
413
414template <typename TAlgebra>
416{
417 return m_vLinSolverCalls[call];
418}
419
420template <typename TAlgebra>
422{
423 return m_vTotalLinSolverSteps[call];
424}
425
426template <typename TAlgebra>
428{
429 return m_vTotalLinSolverSteps[call] / ((double)m_vLinSolverCalls[call]);
430}
431
432template <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
441template <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
450template <typename TAlgebra>
452{
453 return total_linsolver_steps()/((double)total_linsolver_calls());
454}
455
456
457template <typename TAlgebra>
459{
460 m_vLinSolverRates.clear();
461 m_vNonLinSolverRates.clear();
462 m_vLinSolverCalls.clear();
463 m_vTotalLinSolverSteps.clear();
464}
465
466template <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
477template <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
488template <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
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
TAlgebra::matrix_type matrix_type
Matrix type.
Definition newton.h:74
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
TAlgebra::vector_type vector_type
Vector type.
Definition newton.h:71
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)