Line data Source code
1 : /*
2 : * Copyright (c) 2010-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Andreas Vogel, modifications nested Newton: Markus Knodel
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"
41 : #include "lib_disc/function_spaces/grid_function_util.h"
42 : #include "common/util/string_util.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>
60 : NewtonSolver<TAlgebra>::
61 : NewtonSolver(SmartPtr<ILinearOperatorInverse<vector_type> > LinearSolver,
62 : SmartPtr<IConvergenceCheck<vector_type> > spConvCheck,
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 : m_newtonUpdater(SPNULL)
74 : //#if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
75 : // m_newtonUpdater(new NewtonUpdaterGeneric<vector_type>{})
76 : //#endif
77 : {};
78 :
79 : template <typename TAlgebra>
80 0 : NewtonSolver<TAlgebra>::
81 : NewtonSolver() :
82 : m_spLinearSolver(NULL),
83 0 : 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 0 : m_reassembe_J_freq(0),
89 0 : m_dgbCall(0),
90 0 : m_lastNumSteps(0),
91 0 : m_newtonUpdater(SPNULL)
92 : //#if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
93 : // ,
94 : // m_newtonUpdater(new NewtonUpdaterGeneric<vector_type>{})
95 : //#endif
96 0 : {};
97 :
98 : template <typename TAlgebra>
99 0 : NewtonSolver<TAlgebra>::
100 : NewtonSolver(SmartPtr<IOperator<vector_type> > N) :
101 : m_spLinearSolver(NULL),
102 0 : m_spConvCheck(new StdConvCheck<vector_type>(10, 1e-8, 1e-10, true)),
103 : m_spLineSearch(NULL),
104 : m_N(NULL),
105 : m_J(NULL),
106 : m_spAss(NULL),
107 0 : m_reassembe_J_freq(0),
108 0 : m_dgbCall(0),
109 0 : m_lastNumSteps(0),
110 0 : m_newtonUpdater(SPNULL)
111 : //#if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
112 : // ,
113 : // m_newtonUpdater(new NewtonUpdaterGeneric<vector_type>{})
114 : //#endif
115 : {
116 0 : init(N);
117 0 : };
118 :
119 : template <typename TAlgebra>
120 0 : NewtonSolver<TAlgebra>::
121 : NewtonSolver(SmartPtr<IAssemble<TAlgebra> > spAss) :
122 : m_spLinearSolver(NULL),
123 0 : m_spConvCheck(new StdConvCheck<vector_type>(10, 1e-8, 1e-10, true)),
124 : m_spLineSearch(NULL),
125 : m_N(NULL),
126 : m_J(NULL),
127 : m_spAss(NULL),
128 0 : m_reassembe_J_freq(0),
129 0 : m_dgbCall(0),
130 0 : m_lastNumSteps(0),
131 0 : m_newtonUpdater(SPNULL)
132 : //#if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
133 : // ,
134 : // m_newtonUpdater(new NewtonUpdaterGeneric<vector_type>{})
135 : //#endif
136 : {
137 0 : m_spAss = spAss;
138 0 : m_N = SmartPtr<AssembledOperator<TAlgebra> >(new AssembledOperator<TAlgebra>(m_spAss));
139 0 : };
140 :
141 : template <typename TAlgebra>
142 0 : void NewtonSolver<TAlgebra>::
143 : set_convergence_check(SmartPtr<IConvergenceCheck<vector_type> > spConvCheck)
144 : {
145 0 : m_spConvCheck = spConvCheck;
146 0 : m_spConvCheck->set_offset(3);
147 0 : m_spConvCheck->set_symbol('#');
148 0 : m_spConvCheck->set_name("Newton Solver");
149 0 : }
150 :
151 : template <typename TAlgebra>
152 0 : bool NewtonSolver<TAlgebra>::init(SmartPtr<IOperator<vector_type> > N)
153 : {
154 : NEWTON_PROFILE_BEGIN(NewtonSolver_init);
155 0 : m_N = N.template cast_dynamic<AssembledOperator<TAlgebra> >();
156 0 : if(m_N.invalid())
157 0 : UG_THROW("NewtonSolver: currently only works for AssembledDiscreteOperator.");
158 :
159 0 : m_spAss = m_N->discretization();
160 0 : return true;
161 : }
162 :
163 : template <typename TAlgebra>
164 0 : bool NewtonSolver<TAlgebra>::prepare(vector_type& u)
165 : {
166 : // todo: maybe remove this from interface
167 0 : return true;
168 : }
169 :
170 : template <typename TAlgebra>
171 0 : bool NewtonSolver<TAlgebra>::apply(vector_type& u)
172 : {
173 : NEWTON_PROFILE_BEGIN(NewtonSolver_apply);
174 : // increase call count
175 0 : m_dgbCall++;
176 :
177 : // Check for linear solver
178 0 : if(m_spLinearSolver.invalid())
179 0 : UG_THROW("NewtonSolver::apply: Linear Solver not set.");
180 :
181 : // Jacobian
182 0 : if(m_J.invalid() || m_J->discretization() != m_spAss) {
183 0 : m_J = make_sp(new AssembledLinearOperator<TAlgebra>(m_spAss));
184 : }
185 : m_J->set_level(m_N->level());
186 :
187 : // create tmp vectors
188 0 : SmartPtr<vector_type> spD = u.clone_without_values();
189 0 : SmartPtr<vector_type> spC = u.clone_without_values();
190 :
191 : // Set dirichlet values
192 : try{
193 0 : m_N->prepare(u);
194 : }
195 0 : UG_CATCH_THROW("NewtonSolver::prepare: Prepare of Operator failed.");
196 :
197 : // Compute first Defect
198 : try{
199 : NEWTON_PROFILE_BEGIN(NewtonComputeDefect1);
200 0 : m_N->apply(*spD, u);
201 : NEWTON_PROFILE_END();
202 0 : }UG_CATCH_THROW("NewtonSolver::apply: Computation of Start-Defect failed.");
203 :
204 : // loop counts (for the the convergence rate statistics etc.)
205 : int loopCnt = 0;
206 0 : m_lastNumSteps = 0;
207 :
208 : // write start defect for debug
209 : char debug_name_ext[20];
210 0 : if (this->debug_writer_valid())
211 : {
212 : snprintf(debug_name_ext, 20, "_iter%03d", loopCnt);
213 0 : write_debug(*spD, std::string("NEWTON_Defect") + debug_name_ext);
214 0 : write_debug(u, "NEWTON_StartSolution");
215 : }
216 :
217 : // increase offset of output for linear solver
218 0 : const int stdLinOffset = m_spLinearSolver->standard_offset();
219 0 : m_spLinearSolver->convergence_check()->set_offset(stdLinOffset + 3);
220 :
221 : // set info string indicating the used linear solver
222 0 : std::stringstream ss; ss << "(Linear Solver: " << m_spLinearSolver->name() << ")";
223 0 : m_spConvCheck->set_info(ss.str());
224 :
225 : // start convergence check
226 0 : m_spConvCheck->start(*spD);
227 :
228 0 : for(size_t i = 0; i < m_stepUpdate.size(); ++i)
229 0 : m_stepUpdate[i]->update();
230 :
231 : // loop iteration
232 0 : while(!m_spConvCheck->iteration_ended())
233 : {
234 0 : m_lastNumSteps = loopCnt;
235 :
236 : // set c = 0
237 : NEWTON_PROFILE_BEGIN(NewtonSetCorretionZero);
238 : spC->set(0.0);
239 : NEWTON_PROFILE_END();
240 :
241 0 : for(size_t i = 0; i < m_innerStepUpdate.size(); ++i)
242 0 : m_innerStepUpdate[i]->update();
243 :
244 : // Compute Jacobian
245 : try{
246 0 : if(m_reassembe_J_freq == 0 || loopCnt % m_reassembe_J_freq == 0) // if we need to reassemble
247 : {
248 : NEWTON_PROFILE_BEGIN(NewtonComputeJacobian);
249 0 : m_J->init(u);
250 : NEWTON_PROFILE_END();
251 : }
252 0 : }UG_CATCH_THROW("NewtonSolver::apply: Initialization of Jacobian failed.");
253 :
254 : // Write the current Jacobian for debug and prepare the section for the lin. solver
255 0 : if (this->debug_writer_valid())
256 : {
257 0 : write_debug(m_J->get_matrix(), std::string("NEWTON_Jacobian") + debug_name_ext);
258 0 : this->enter_debug_writer_section(std::string("NEWTON_LinSolver") + debug_name_ext);
259 : }
260 :
261 : // Init Jacobi Inverse
262 : try{
263 : NEWTON_PROFILE_BEGIN(NewtonPrepareLinSolver);
264 0 : if(!m_spLinearSolver->init(m_J, u))
265 : {
266 : UG_LOG("ERROR in 'NewtonSolver::apply': Cannot init Inverse Linear "
267 : "Operator for Jacobi-Operator.\n");
268 : return false;
269 : }
270 : NEWTON_PROFILE_END();
271 0 : }UG_CATCH_THROW("NewtonSolver::apply: Initialization of Linear Solver failed.");
272 :
273 : // Solve Linearized System
274 : try{
275 : NEWTON_PROFILE_BEGIN(NewtonApplyLinSolver);
276 0 : if(!m_spLinearSolver->apply(*spC, *spD))
277 : {
278 : UG_LOG("ERROR in 'NewtonSolver::apply': Cannot apply Inverse Linear "
279 : "Operator for Jacobi-Operator.\n");
280 : return false;
281 : }
282 : NEWTON_PROFILE_END();
283 0 : }UG_CATCH_THROW("NewtonSolver::apply: Application of Linear Solver failed.");
284 :
285 0 : this->leave_debug_writer_section();
286 :
287 : // store convergence history
288 0 : const int numSteps = m_spLinearSolver->step();
289 0 : if(loopCnt >= (int)m_vTotalLinSolverSteps.size()) m_vTotalLinSolverSteps.resize(loopCnt+1);
290 0 : if(loopCnt >= (int)m_vLinSolverCalls.size()) m_vLinSolverCalls.resize(loopCnt+1, 0);
291 0 : if(loopCnt >= (int)m_vLinSolverRates.size()) m_vLinSolverRates.resize(loopCnt+1, 0);
292 0 : m_vTotalLinSolverSteps[loopCnt] += numSteps;
293 0 : m_vLinSolverCalls[loopCnt] += 1;
294 0 : m_vLinSolverRates[loopCnt] += m_spLinearSolver->convergence_check()->avg_rate();
295 :
296 : try{
297 : // Line Search
298 0 : if(m_spLineSearch.valid())
299 : {
300 0 : m_spLineSearch->set_offset(" # ");
301 : NEWTON_PROFILE_BEGIN(NewtonLineSearch);
302 :
303 : //#if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
304 :
305 0 : if( m_newtonUpdater != SPNULL )
306 0 : m_spLineSearch->setNewtonUpdater(m_newtonUpdater);
307 :
308 : //#endif
309 0 : if(!m_spLineSearch->search(m_N, u, *spC, *spD, m_spConvCheck->defect()))
310 : {
311 : UG_LOG("ERROR in 'NewtonSolver::apply': "
312 : "Newton Solver did not converge.\n");
313 : return false;
314 : }
315 : NEWTON_PROFILE_END();
316 : }
317 : // No line search: Compute new defect
318 : else
319 : {
320 0 : if( m_newtonUpdater != SPNULL )
321 : {
322 : //#if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
323 :
324 0 : if( ! m_newtonUpdater->updateSolution(u,*spC) )
325 : {
326 : UG_LOG("ERROR in 'NewtonSolver::apply': "
327 : "Newton Update did not work.\n");
328 : // TODO FIXME was macht conv check update? wie kriege ich hier einfach ein riesiges Residuum rein, ohne was zu rechnen?
329 : return false;
330 : }
331 :
332 0 : if( ! m_newtonUpdater->tellAndFixUpdateEvents(u) )
333 : {
334 : UG_LOG("unable to fix local Newton updates" << std::endl );
335 : return false;
336 : }
337 :
338 : }
339 : else
340 : {
341 : //#else
342 : // update solution
343 0 : u -= *spC;
344 : }
345 : //#endif
346 :
347 : // compute new Defect
348 : NEWTON_PROFILE_BEGIN(NewtonComputeDefect);
349 0 : m_N->prepare(u);
350 0 : m_N->apply(*spD, u);
351 : NEWTON_PROFILE_END();
352 : }
353 0 : }UG_CATCH_THROW("NewtonSolver::apply: Line Search update failed.");
354 :
355 :
356 : // update counter
357 0 : loopCnt++;
358 :
359 : // check convergence
360 0 : m_spConvCheck->update(*spD);
361 0 : if(loopCnt-1 >= (int)m_vNonLinSolverRates.size()) m_vNonLinSolverRates.resize(loopCnt, 0);
362 0 : m_vNonLinSolverRates[loopCnt-1] += m_spConvCheck->rate();
363 :
364 : // write defect for debug
365 0 : if (this->debug_writer_valid())
366 : {
367 : snprintf(debug_name_ext, 20, "_iter%03d", loopCnt);
368 0 : write_debug(*spD, std::string("NEWTON_Defect") + debug_name_ext);
369 0 : write_debug(*spC, std::string("NEWTON_Correction") + debug_name_ext);
370 0 : write_debug(u, std::string("NEWTON_Solution") + debug_name_ext);
371 : }
372 : }
373 :
374 : // reset offset of output for linear solver to previous value
375 0 : m_spLinearSolver->convergence_check()->set_offset(stdLinOffset);
376 :
377 0 : return m_spConvCheck->post();
378 0 : }
379 :
380 : template <typename TAlgebra>
381 0 : void NewtonSolver<TAlgebra>::print_average_convergence() const
382 : {
383 : UG_LOG("\nNewton solver convergence history:\n");
384 : UG_LOG("Newton Step | Num Calls | Total Lin Iters | Avg Lin Iters | Avg Nonlin Rates | Avg Lin Rates \n");
385 : int allCalls = 0, allLinSteps = 0;
386 : number allLinRatesProduct = 1.0, allNonLinRatesProduct = 1.0;
387 0 : for(int call = 0; call < (int)m_vLinSolverCalls.size(); ++call)
388 : {
389 0 : UG_LOG( " " << std::setw(10) << call+1 << " | ");
390 0 : UG_LOG(std::setw(9) << m_vLinSolverCalls[call] << " | ");
391 0 : allCalls += m_vLinSolverCalls[call];
392 0 : UG_LOG(std::setw(15) << m_vTotalLinSolverSteps[call] << " | ");
393 0 : allLinSteps += m_vTotalLinSolverSteps[call];
394 0 : UG_LOG(std::setw(13) << std::setprecision(2) << std::fixed << m_vTotalLinSolverSteps[call] / (double)m_vLinSolverCalls[call] << " | ");
395 0 : allNonLinRatesProduct *= pow((number)m_vNonLinSolverRates[call]/(double)m_vLinSolverCalls[call],(double)m_vLinSolverCalls[call]);
396 0 : UG_LOG(std::setw(16) << std::setprecision(6) << std::scientific << m_vNonLinSolverRates[call] / (double)m_vLinSolverCalls[call] << " | ");
397 0 : allLinRatesProduct *= (number)std::pow((number)m_vLinSolverRates[call]/(double)m_vLinSolverCalls[call],(number)m_vTotalLinSolverSteps[call]);
398 0 : UG_LOG(std::setw(13) << std::setprecision(6) << std::scientific << m_vLinSolverRates[call] / (double)m_vLinSolverCalls[call]);
399 : UG_LOG("\n");
400 : }
401 : UG_LOG( " all | ");
402 0 : UG_LOG(std::setw(9) << allCalls << " | ");
403 0 : UG_LOG(std::setw(15) << allLinSteps << " | ");
404 0 : UG_LOG(std::setw(13) << std::setprecision(2) << std::fixed << allLinSteps / (number)allCalls << " | ");
405 0 : UG_LOG(std::setw(16) << std::setprecision(6) << std::scientific << std::pow((number)allNonLinRatesProduct,(number)1.0/(number)allCalls) << " | ");
406 0 : UG_LOG(std::setw(13) << std::setprecision(6) << std::scientific << std::pow((number)allLinRatesProduct,(number)1.0/(number)allLinSteps));
407 : UG_LOG("\n");
408 0 : }
409 :
410 : template <typename TAlgebra>
411 0 : size_t NewtonSolver<TAlgebra>::num_newton_steps() const
412 : {
413 0 : return m_vLinSolverCalls.size();
414 : }
415 :
416 : template <typename TAlgebra>
417 0 : int NewtonSolver<TAlgebra>::num_linsolver_calls(size_t call) const
418 : {
419 0 : return m_vLinSolverCalls[call];
420 : }
421 :
422 : template <typename TAlgebra>
423 0 : int NewtonSolver<TAlgebra>::num_linsolver_steps(size_t call) const
424 : {
425 0 : return m_vTotalLinSolverSteps[call];
426 : }
427 :
428 : template <typename TAlgebra>
429 0 : double NewtonSolver<TAlgebra>::average_linear_steps(size_t call) const
430 : {
431 0 : return m_vTotalLinSolverSteps[call] / ((double)m_vLinSolverCalls[call]);
432 : }
433 :
434 : template <typename TAlgebra>
435 0 : int NewtonSolver<TAlgebra>::total_linsolver_calls() const
436 : {
437 : int allCalls = 0;
438 0 : for(size_t call = 0; call < m_vLinSolverCalls.size(); ++call)
439 0 : allCalls += m_vLinSolverCalls[call];
440 0 : return allCalls;
441 : }
442 :
443 : template <typename TAlgebra>
444 0 : int NewtonSolver<TAlgebra>::total_linsolver_steps() const
445 : {
446 : int allSteps = 0;
447 0 : for(size_t call = 0; call < m_vLinSolverCalls.size(); ++call)
448 0 : allSteps += m_vTotalLinSolverSteps[call];
449 0 : return allSteps;
450 : }
451 :
452 : template <typename TAlgebra>
453 0 : double NewtonSolver<TAlgebra>::total_average_linear_steps() const
454 : {
455 0 : return total_linsolver_steps()/((double)total_linsolver_calls());
456 : }
457 :
458 :
459 : template <typename TAlgebra>
460 0 : void NewtonSolver<TAlgebra>::clear_average_convergence()
461 : {
462 : m_vLinSolverRates.clear();
463 : m_vNonLinSolverRates.clear();
464 : m_vLinSolverCalls.clear();
465 : m_vTotalLinSolverSteps.clear();
466 0 : }
467 :
468 : template <typename TAlgebra>
469 0 : void NewtonSolver<TAlgebra>::write_debug(const vector_type& vec, std::string name)
470 : {
471 : // add call count to name
472 0 : char ext[20]; snprintf(ext, 20, "_call%03d", m_dgbCall);
473 :
474 : // write
475 : typedef DebugWritingObject<TAlgebra> base_writer_type;
476 0 : base_writer_type::write_debug(vec, name + ext);
477 0 : }
478 :
479 : template <typename TAlgebra>
480 0 : void NewtonSolver<TAlgebra>::write_debug(const matrix_type& mat, std::string name)
481 : {
482 : // add call count to name
483 0 : char ext[20]; snprintf(ext, 20, "_call%03d", m_dgbCall);
484 :
485 : // write
486 : typedef DebugWritingObject<TAlgebra> base_writer_type;
487 0 : base_writer_type::write_debug(mat, name + ext);
488 0 : }
489 :
490 : template <typename TAlgebra>
491 0 : std::string NewtonSolver<TAlgebra>::config_string() const
492 : {
493 0 : std::stringstream ss;
494 0 : ss << "NewtonSolver\n";
495 0 : ss << " LinearSolver: ";
496 0 : if(m_spLinearSolver.valid()) ss << ConfigShift(m_spLinearSolver->config_string()) << "\n";
497 0 : else ss << " NOT SET!\n";
498 0 : ss << " ConvergenceCheck: ";
499 0 : if(m_spConvCheck.valid()) ss << ConfigShift(m_spConvCheck->config_string()) << "\n";
500 0 : else ss << " NOT SET!\n";
501 0 : ss << " LineSearch: ";
502 0 : if(m_spLineSearch.valid()) ss << ConfigShift(m_spLineSearch->config_string()) << "\n";
503 0 : else ss << " not set.\n";
504 0 : if(m_reassembe_J_freq != 0) ss << " Reassembling Jacobian only once per " << m_reassembe_J_freq << " step(s)\n";
505 0 : return ss.str();
506 0 : }
507 :
508 :
509 : }
510 :
511 : #endif /* __H__UG__LIB_DISC__OPERATOR__NON_LINEAR_OPERATOR__NEWTON_SOLVER__NEWTON_IMPL__ */
512 :
|
