DOpE
fractional_step_theta_step_newtonsolver.h
Go to the documentation of this file.
1 
24 #ifndef FRACTIONAL_STEP_THETA_STEP_NEWTON_SOLVER_H_
25 #define FRACTIONAL_STEP_THETA_STEP_NEWTON_SOLVER_H_
26 
27 #include <deal.II/lac/vector.h>
28 #include <deal.II/lac/block_sparsity_pattern.h>
29 #include <deal.II/lac/block_sparse_matrix.h>
30 #include <deal.II/numerics/vector_tools.h>
31 
32 #include <vector>
33 #include <iostream>
34 #include <fstream>
35 #include <iomanip>
36 
37 #include "parameterreader.h"
38 
39 
40 
41 
42 namespace DOpE
43 {
54  template <typename INTEGRATOR, typename LINEARSOLVER, typename VECTOR>
55  class FractionalStepThetaStepNewtonSolver : public LINEARSOLVER
56  {
57  public:
65  FractionalStepThetaStepNewtonSolver(INTEGRATOR &integrator, ParameterReader &param_reader);
66  ~FractionalStepThetaStepNewtonSolver();
67 
68  /******************************************************/
69 
75  static void declare_params(ParameterReader &param_reader);
76 
77  /******************************************************/
78 
83  template<typename PROBLEM>
84  void ReInit(PROBLEM& pde);
85 
86  /******************************************************/
87 
117  template<typename PROBLEM>
118  bool NonlinearSolve(PROBLEM& pde, const VECTOR &last_time_solution, VECTOR &solution,
119  bool apply_boundary_values=true,
120  bool force_matrix_build=false, int priority = 5, std::string algo_level = "\t\t ");
121 
122  /******************************************************/
123 
152  template<typename PROBLEM>
153  bool NonlinearSolve_Initial(PROBLEM& pde, VECTOR &solution, bool apply_boundary_values=true,
154  bool force_matrix_build=false, int priority = 5, std::string algo_level = "\t\t ");
155 
156  /******************************************************/
157 
168  template<typename PROBLEM>
169  void NonlinearLastTimeEvals(PROBLEM& pde, const VECTOR &last_time_solution, VECTOR &residual);
170 
171  protected:
172 
173  inline INTEGRATOR& GetIntegrator();
174 
175  private:
176  INTEGRATOR &integrator_;
177 
178  bool build_matrix_;
179 
180  double nonlinear_global_tol_, nonlinear_tol_, nonlinear_rho_;
181  double linesearch_rho_;
182  int nonlinear_maxiter_, line_maxiter_;
183  };
184 
185  /**********************************Implementation*******************************************/
186 
187 template <typename INTEGRATOR, typename LINEARSOLVER, typename VECTOR>
188  void FractionalStepThetaStepNewtonSolver<INTEGRATOR,LINEARSOLVER, VECTOR>
189  ::declare_params(ParameterReader &param_reader)
190  {
191  param_reader.SetSubsection("newtonsolver parameters");
192  param_reader.declare_entry("nonlinear_global_tol", "1.e-12",Patterns::Double(0),"global tolerance for the newton iteration");
193  param_reader.declare_entry("nonlinear_tol", "1.e-10",Patterns::Double(0),"relative tolerance for the newton iteration");
194  param_reader.declare_entry("nonlinear_maxiter", "10",Patterns::Integer(0),"maximal number of newton iterations");
195  param_reader.declare_entry("nonlinear_rho", "0.1",Patterns::Double(0),"minimal newton reduction, if actual reduction is less, matrix is rebuild ");
196 
197  param_reader.declare_entry("line_maxiter", "4",Patterns::Integer(0),"maximal number of linesearch steps");
198  param_reader.declare_entry("linesearch_rho", "0.9",Patterns::Double(0),"reduction rate for the linesearch damping paramete");
199 
200  LINEARSOLVER::declare_params(param_reader);
201  }
202 
203  /*******************************************************************************************/
204 
205  template <typename INTEGRATOR, typename LINEARSOLVER, typename VECTOR>
206  FractionalStepThetaStepNewtonSolver<INTEGRATOR,LINEARSOLVER, VECTOR>
207  ::FractionalStepThetaStepNewtonSolver(INTEGRATOR &integrator, ParameterReader &param_reader)
208  : LINEARSOLVER(param_reader), integrator_(integrator)
209  {
210  param_reader.SetSubsection("newtonsolver parameters");
211  nonlinear_global_tol_ = param_reader.get_double ("nonlinear_global_tol");
212  nonlinear_tol_ = param_reader.get_double ("nonlinear_tol");
213  nonlinear_maxiter_ = param_reader.get_integer ("nonlinear_maxiter");
214  nonlinear_rho_ = param_reader.get_double ("nonlinear_rho");
215 
216  line_maxiter_ = param_reader.get_integer ("line_maxiter");
217  linesearch_rho_ = param_reader.get_double ("linesearch_rho");
218 
219  }
220 
221  /*******************************************************************************************/
222 
223  template <typename INTEGRATOR, typename LINEARSOLVER, typename VECTOR>
224  FractionalStepThetaStepNewtonSolver<INTEGRATOR,LINEARSOLVER, VECTOR>
225  ::~FractionalStepThetaStepNewtonSolver()
226  {
227  }
228 
229  /*******************************************************************************************/
230  template <typename INTEGRATOR, typename LINEARSOLVER, typename VECTOR>
231  template<typename PROBLEM>
232  void FractionalStepThetaStepNewtonSolver<INTEGRATOR,LINEARSOLVER, VECTOR>
233  ::ReInit(PROBLEM& pde)
234  {
235  LINEARSOLVER::ReInit(pde);
236  }
237 
238  /*******************************************************************************************/
239  template <typename INTEGRATOR, typename LINEARSOLVER, typename VECTOR>
240  template<typename PROBLEM>
241  void FractionalStepThetaStepNewtonSolver<INTEGRATOR,LINEARSOLVER, VECTOR>
242  ::NonlinearLastTimeEvals(PROBLEM& pde, const VECTOR &last_time_solution
243  ,VECTOR &residual)
244  {
245  VECTOR tmp_residual;
246  tmp_residual.reinit(residual);
247  residual =0.;
248  GetIntegrator().AddDomainData("last_newton_solution",&last_time_solution);
249  GetIntegrator().AddDomainData("last_time_solution",&last_time_solution);
250  pde.SetStepPart("Old_for_1st_cycle");
251  GetIntegrator().ComputeNonlinearLhs(pde,residual);
252  GetIntegrator().ComputeNonlinearRhs(pde,tmp_residual);
253  tmp_residual *= -1;
254  residual += tmp_residual;
255 
256  GetIntegrator().DeleteDomainData("last_newton_solution");
257  GetIntegrator().DeleteDomainData("last_time_solution");
258  }
259  /*******************************************************************************************/
260 
261  template <typename INTEGRATOR, typename LINEARSOLVER, typename VECTOR>
262  template<typename PROBLEM>
263  bool FractionalStepThetaStepNewtonSolver<INTEGRATOR,LINEARSOLVER, VECTOR>
264  ::NonlinearSolve_Initial(PROBLEM& pde, VECTOR &solution, bool apply_boundary_values,
265  bool force_matrix_build, int priority, std::string algo_level)
266  {
267  bool build_matrix = force_matrix_build;
268  VECTOR residual;
269  VECTOR du;
270  std::stringstream out;
271  pde.GetOutputHandler()->InitNewtonOut(out);
272 
273  du.reinit(solution);
274  residual.reinit(solution);
275 
276  if(apply_boundary_values)
277  {
278  GetIntegrator().ApplyInitialBoundaryValues(pde,solution);
279  }
280  pde.GetDoFConstraints().distribute(solution);
281 
282  GetIntegrator().AddDomainData("last_newton_solution",&solution);
283 
284  GetIntegrator().ComputeNonlinearResidual(pde,residual);
285  residual *= -1.;
286 
287  pde.GetOutputHandler()->SetIterationNumber(0,"PDENewton");
288  pde.GetOutputHandler()->Write(residual,"Residual"+pde.GetType(),pde.GetDoFType());
289 
290  double res = residual.linfty_norm();
291  double firstres = res;
292  double lastres = res;
293 
294 
295  out<< algo_level << "Newton step: " <<0<<"\t Residual (abs.): "
296  <<pde.GetOutputHandler()->ZeroTolerance(res, 1.0)
297  <<"\n";
298 
299  out<< algo_level << "Newton step: " <<0<<"\t Residual (rel.): " << std::scientific << firstres/firstres;
300 
301 
302  pde.GetOutputHandler()->Write(out,priority);
303 
304  int iter=0;
305  while(res > nonlinear_global_tol_ && res > firstres * nonlinear_tol_)
306  {
307  iter++;
308 
309  if(iter > nonlinear_maxiter_)
310  {
311  GetIntegrator().DeleteDomainData("last_newton_solution");
312  throw DOpEIterationException("Iteration count exceeded bounds!","InstatStepNewtonSolver::NonlinearSolve_Initial");
313  }
314 
315  pde.GetOutputHandler()->SetIterationNumber(iter,"PDENewton");
316 
317  LINEARSOLVER::Solve(pde,GetIntegrator(),residual,du,build_matrix);
318  build_matrix = false;
319 
320  //Linesearch
321  {
322  solution += du;
323  GetIntegrator().ComputeNonlinearResidual(pde,residual);
324  residual *= -1.;
325 
326  pde.GetOutputHandler()->Write(residual,"Residual"+pde.GetType(),pde.GetDoFType());
327  pde.GetOutputHandler()->Write(du,"Update"+pde.GetType(),pde.GetDoFType());
328 
329  double newres = residual.linfty_norm();
330  int lineiter=0;
331  double rho = linesearch_rho_;
332  double alpha=1;
333 
334  while(newres > res)
335  {
336  out<< algo_level << "Newton step: " <<iter<<"\t Residual (rel.): "
337  <<pde.GetOutputHandler()->ZeroTolerance(newres/firstres, 1.0)
338  << "\t LineSearch {"<<lineiter<<"} ";
339 
340  pde.GetOutputHandler()->Write(out,priority+1);
341 
342  lineiter++;
343  if(lineiter > line_maxiter_)
344  {
345  GetIntegrator().DeleteDomainData("last_newton_solution");
346  throw DOpEIterationException("Line-Iteration count exceeded bounds!","InstatStepNewtonSolver::NonlinearSolve_Initial");
347  }
348  solution.add(alpha*(rho-1.),du);
349  alpha*= rho;
350 
351  GetIntegrator().ComputeNonlinearResidual(pde,residual);
352  residual *= -1.;
353  pde.GetOutputHandler()->Write(residual,"Residual"+pde.GetType(),pde.GetDoFType());
354 
355  newres = residual.linfty_norm();
356 
357  }
358  if(res/lastres > nonlinear_rho_)
359  {
360  build_matrix=true;
361  }
362  lastres=res;
363  res=newres;
364 
365  out << algo_level
366  << "Newton step: "
367  <<iter
368  <<"\t Residual (rel.): "
369  << pde.GetOutputHandler()->ZeroTolerance(res/firstres, 1.0)
370  << "\t LineSearch {"
371  <<lineiter
372  <<"} ";
373 
374 
375  pde.GetOutputHandler()->Write(out,priority);
376 
377  }//End of Linesearch
378  }
379  GetIntegrator().DeleteDomainData("last_newton_solution");
380 
381  return build_matrix;
382  }
383 
384  /*******************************************************************************************/
385  template <typename INTEGRATOR, typename LINEARSOLVER, typename VECTOR>
386  template<typename PROBLEM>
387  bool FractionalStepThetaStepNewtonSolver<INTEGRATOR,LINEARSOLVER, VECTOR>
388  ::NonlinearSolve(PROBLEM& pde,
389  const VECTOR &last_time_solution,
390  VECTOR &solution,
391  bool apply_boundary_values,
392  bool force_matrix_build,
393  int priority,
394  std::string /*algo_level*/)
395  {
396 
397  bool build_matrix = force_matrix_build;
398  VECTOR residual, time_residual, tmp_residual;
399  VECTOR du, tmp_last_time_solution;
400  std::stringstream out;
401  pde.GetOutputHandler()->InitNewtonOut(out);
402 
403  double res = 0.0;
404  double firstres = 0.0;
405  double lastres = 0.0;
406 
407  du.reinit(solution);
408  residual.reinit(solution);
409  time_residual.reinit(solution);
410  tmp_residual.reinit(solution);
411  tmp_last_time_solution.reinit(solution);
412 
413  //Transfer from previous timestep
414  residual +=solution;
415  // last_time_solution is very good starting value
416  solution = last_time_solution;
417 
418 
419  // First cycle of FS-scheme
420  if(apply_boundary_values)
421  {
422  GetIntegrator().ApplyInitialBoundaryValues(pde,solution);
423  }
424 
425  // Righthandside for the current timestep f^{n+1}
426  GetIntegrator().AddDomainData("last_time_solution",&last_time_solution);
427  GetIntegrator().AddDomainData("last_newton_solution",&solution);
428  pde.SetStepPart("New_for_1st_and_3rd_cycle");
429  GetIntegrator().ComputeNonlinearRhs(pde,tmp_residual);
430  tmp_residual *= -1;
431  residual += tmp_residual;
432 
433  // Save the part of the residual which is independent of u^{n+1}
434  time_residual = residual;
435  time_residual *=-1;
436 
437  // Calculate the "real" residual in the current timestep
438  GetIntegrator().ComputeNonlinearLhs(pde,tmp_residual); // modi, new
439  residual += tmp_residual;
440 
441  residual *=-1.; // due to A(U)(\psi) = - A(U)(du,\psi)
442 
443  pde.GetOutputHandler()->SetIterationNumber(0,"PDENewton");
444  pde.GetOutputHandler()->Write(residual,"Residual"+pde.GetType(),pde.GetDoFType());
445 
446  res = residual.linfty_norm();
447  firstres = res;
448  lastres = res;
449  int iter=0;
450 
451  out<<"\t\t Newton step: " <<0<<"\t Residual (abs.): "
452  <<pde.GetOutputHandler()->ZeroTolerance(res, 1.0)
453  <<"\n";
454  out<<"\t\t Newton step: " <<0<<"\t Residual (rel.): "<< std::scientific << firstres/firstres;
455 
456 
457 
458 
459  pde.GetOutputHandler()->Write(out,priority);
460  while(res > nonlinear_global_tol_ && res > firstres * nonlinear_tol_)
461  {
462  iter++;
463  if(iter > nonlinear_maxiter_)
464  {
465  throw DOpEIterationException("Iteration count exceeded bounds!","StatSolver::NonlinearSolve");
466  }
467 
468  pde.GetOutputHandler()->SetIterationNumber(iter,"PDENewton");
469  LINEARSOLVER::Solve(pde,GetIntegrator(),residual,du,build_matrix);
470  build_matrix = false;
471  //Linesearch
472  {
473  solution += du;
474  GetIntegrator().ComputeNonlinearLhs(pde,residual);
475  residual -= time_residual;
476  residual *= -1.;
477  pde.GetOutputHandler()->Write(residual,"Residual"+pde.GetType(),pde.GetDoFType());
478 
479  double newres = residual.linfty_norm();
480  int lineiter=0;
481  double rho = linesearch_rho_;
482  double alpha=1;
483 
484  while(newres > res)
485  {
486  out<<"\t\t\t Linesearch step: " <<lineiter<<"\t Residual (rel.): "
487  <<pde.GetOutputHandler()->ZeroTolerance(newres/firstres, 1.0);
488 
489  pde.GetOutputHandler()->Write(out,priority+1);
490  lineiter++;
491  if(lineiter > line_maxiter_)
492  {
493  throw DOpEIterationException("Line-Iteration count exceeded bounds!","StatSolver::NonlinearSolve");
494  }
495  solution.add(alpha*(rho-1.),du);
496  alpha*= rho;
497 
498  GetIntegrator().ComputeNonlinearLhs(pde,residual);
499  residual -= time_residual;
500  residual *= -1.;
501  pde.GetOutputHandler()->Write(residual,"Residual"+pde.GetType(),pde.GetDoFType());
502 
503  newres = residual.linfty_norm();
504 
505  }
506  if(res/lastres > nonlinear_rho_)
507  {
508  build_matrix=true;
509  }
510  lastres=res;
511  res=newres;
512 
513  out<<"\t\t Newton step: " <<iter<<"\t Residual (rel.): "
514  << pde.GetOutputHandler()->ZeroTolerance(res/firstres, 1.0)
515  << "\t LineSearch {"<<lineiter<<"} ";
516  pde.GetOutputHandler()->Write(out,priority);
517 
518  }//End of Linesearch
519  }
520  GetIntegrator().DeleteDomainData("last_time_solution");
521  GetIntegrator().DeleteDomainData("last_newton_solution");
522 
523 
524  // 2nd cycle of FS-scheme
525  tmp_last_time_solution = 0;
526  tmp_last_time_solution = solution;
527 
528  // needs to be set here, otherwise the
529  // values in element equation are not set!
530  GetIntegrator().AddDomainData("last_time_solution",&tmp_last_time_solution);
531 
532  // Calculate residual parts corresponding to the last time-step
533  GetIntegrator().AddDomainData("last_newton_solution",&tmp_last_time_solution);
534  pde.SetStepPart("Old_for_2nd_cycle");
535  GetIntegrator().ComputeNonlinearLhs(pde,residual);
536 
537  GetIntegrator().ComputeNonlinearRhs(pde,tmp_residual);
538  tmp_residual *= -1;
539  residual += tmp_residual;
540 
541  GetIntegrator().DeleteDomainData("last_newton_solution");
542  // Old time solution finished
543 
544  if(apply_boundary_values)
545  {
546  GetIntegrator().ApplyInitialBoundaryValues(pde,solution);
547  }
548 
549 
550  // Righthandside for the current timestep f^{n+1}
551  GetIntegrator().AddDomainData("last_newton_solution",&solution);
552  pde.SetStepPart("New_for_2nd_cycle");
553  GetIntegrator().ComputeNonlinearRhs(pde,tmp_residual);
554  tmp_residual *= -1;
555  residual += tmp_residual;
556 
557  // Save the part of the residual which is independent of u^{n+1}
558  time_residual = residual;
559  time_residual *=-1;
560 
561  // Calculate the "real" residual in the current timestep
562  GetIntegrator().ComputeNonlinearLhs(pde,tmp_residual); // modi, new
563  residual += tmp_residual;
564 
565  residual *=-1.; // due to A(U)(\psi) = - A(U)(du,\psi)
566 
567  pde.GetOutputHandler()->SetIterationNumber(0,"PDENewton");
568  pde.GetOutputHandler()->Write(residual,"Residual"+pde.GetType(),pde.GetDoFType());
569 
570  res = residual.linfty_norm();
571  firstres = res;
572  lastres = res;
573  iter=0;
574 
575  out<<"\t\t Newton step: " <<0<<"\t Residual (abs.): "
576  <<pde.GetOutputHandler()->ZeroTolerance(res, 1.0)
577  <<"\n";
578 
579  out<<"\t\t Newton step: " <<0<<"\t Residual (rel.): "<< std::scientific << firstres/firstres;
580 
581 
582  pde.GetOutputHandler()->Write(out,priority);
583  while(res > nonlinear_global_tol_ && res > firstres * nonlinear_tol_)
584  {
585  iter++;
586  if(iter > nonlinear_maxiter_)
587  {
588  throw DOpEIterationException("Iteration count exceeded bounds!","StatSolver::NonlinearSolve");
589  }
590 
591  pde.GetOutputHandler()->SetIterationNumber(iter,"PDENewton");
592  LINEARSOLVER::Solve(pde,GetIntegrator(),residual,du,build_matrix);
593  build_matrix = false;
594  //Linesearch
595  {
596  solution += du;
597  GetIntegrator().ComputeNonlinearLhs(pde,residual);
598  residual -= time_residual;
599  residual *= -1.;
600  pde.GetOutputHandler()->Write(residual,"Residual"+pde.GetType(),pde.GetDoFType());
601 
602  double newres = residual.linfty_norm();
603  int lineiter=0;
604  double rho = linesearch_rho_;
605  double alpha=1;
606 
607  while(newres > res)
608  {
609  out<<"\t\t\t Linesearch step: " <<lineiter<<"\t Residual (rel.): "
610  <<pde.GetOutputHandler()->ZeroTolerance(newres/firstres, 1.0);
611 
612  pde.GetOutputHandler()->Write(out,priority+1);
613  lineiter++;
614  if(lineiter > line_maxiter_)
615  {
616  throw DOpEIterationException("Line-Iteration count exceeded bounds!","StatSolver::NonlinearSolve");
617  }
618  solution.add(alpha*(rho-1.),du);
619  alpha*= rho;
620 
621  GetIntegrator().ComputeNonlinearLhs(pde,residual);
622  residual -= time_residual;
623  residual *= -1.;
624  pde.GetOutputHandler()->Write(residual,"Residual"+pde.GetType(),pde.GetDoFType());
625 
626  newres = residual.linfty_norm();
627 
628  }
629  if(res/lastres > nonlinear_rho_)
630  {
631  build_matrix=true;
632  }
633  lastres=res;
634  res=newres;
635 
636  out<<"\t\t Newton step: " <<iter<<"\t Residual (rel.): "
637  << pde.GetOutputHandler()->ZeroTolerance(res/firstres, 1.0)
638  << "\t LineSearch {"<<lineiter<<"} ";
639 
640  pde.GetOutputHandler()->Write(out,priority);
641 
642  }//End of Linesearch
643  }
644  GetIntegrator().DeleteDomainData("last_time_solution");
645  GetIntegrator().DeleteDomainData("last_newton_solution");
646 
647 
648  // 3rd cycle of FS-scheme
649  tmp_last_time_solution = 0;
650  tmp_last_time_solution = solution;
651 
652  // needs to be set here, otherwise the
653  // values in element equation are not set!
654  GetIntegrator().AddDomainData("last_time_solution",&tmp_last_time_solution);
655 
656  // Calculate residual parts corresponding to the last time-step
657  GetIntegrator().AddDomainData("last_newton_solution",&tmp_last_time_solution);
658  pde.SetStepPart("Old_for_3rd_cycle");
659  GetIntegrator().ComputeNonlinearLhs(pde,residual);
660 
661  GetIntegrator().ComputeNonlinearRhs(pde,tmp_residual);
662  tmp_residual *= -1;
663  residual += tmp_residual;
664 
665  GetIntegrator().DeleteDomainData("last_newton_solution");
666  // Old time solution finished
667 
668  if(apply_boundary_values)
669  {
670  GetIntegrator().ApplyInitialBoundaryValues(pde,solution);
671  }
672 
673  // Righthandside for the current timestep f^{n+1}
674  GetIntegrator().AddDomainData("last_newton_solution",&solution);
675  pde.SetStepPart("New_for_1st_and_3rd_cycle");
676  GetIntegrator().ComputeNonlinearRhs(pde,tmp_residual);
677  tmp_residual *= -1;
678  residual += tmp_residual;
679 
680  // Save the part of the residual which is independent of u^{n+1}
681  time_residual = residual;
682  time_residual *=-1;
683 
684  // Calculate the "real" residual in the current timestep
685  GetIntegrator().ComputeNonlinearLhs(pde,tmp_residual); // modi, new
686  residual += tmp_residual;
687 
688  residual *=-1.; // wg. A(U)(\psi) = - A(U)(du,\psi)
689 
690  pde.GetOutputHandler()->SetIterationNumber(0,"PDENewton");
691  pde.GetOutputHandler()->Write(residual,"Residual"+pde.GetType(),pde.GetDoFType());
692 
693  res = residual.linfty_norm();
694  firstres = res;
695  lastres = res;
696  iter=0;
697 
698  out<<"\t\t Newton step: " <<0<<"\t Residual (abs.): "
699  <<pde.GetOutputHandler()->ZeroTolerance(res, 1.0)
700  <<"\n";
701 
702  out<<"\t\t Newton step: " <<0<<"\t Residual (rel.): "<< std::scientific << firstres/firstres;
703 
704 
705  pde.GetOutputHandler()->Write(out,priority);
706  while(res > nonlinear_global_tol_ && res > firstres * nonlinear_tol_)
707  {
708  iter++;
709  if(iter > nonlinear_maxiter_)
710  {
711  throw DOpEIterationException("Iteration count exceeded bounds!","StatSolver::NonlinearSolve");
712  }
713 
714  pde.GetOutputHandler()->SetIterationNumber(iter,"PDENewton");
715  LINEARSOLVER::Solve(pde,GetIntegrator(),residual,du,build_matrix);
716  build_matrix = false;
717  //Linesearch
718  {
719  solution += du;
720  GetIntegrator().ComputeNonlinearLhs(pde,residual);
721  residual -= time_residual;
722  residual *= -1.;
723  pde.GetOutputHandler()->Write(residual,"Residual"+pde.GetType(),pde.GetDoFType());
724 
725  double newres = residual.linfty_norm();
726  int lineiter=0;
727  double rho = linesearch_rho_;
728  double alpha=1;
729 
730  while(newres > res)
731  {
732  out<<"\t\t\t Linesearch step: " <<lineiter<<"\t Residual (rel.): "
733  <<pde.GetOutputHandler()->ZeroTolerance(newres/firstres, 1.0);
734 
735  pde.GetOutputHandler()->Write(out,priority+1);
736  lineiter++;
737  if(lineiter > line_maxiter_)
738  {
739  throw DOpEIterationException("Line-Iteration count exceeded bounds!","StatSolver::NonlinearSolve");
740  }
741  solution.add(alpha*(rho-1.),du);
742  alpha*= rho;
743 
744  GetIntegrator().ComputeNonlinearLhs(pde,residual);
745  residual -= time_residual;
746  residual *= -1.;
747  pde.GetOutputHandler()->Write(residual,"Residual"+pde.GetType(),pde.GetDoFType());
748 
749  newres = residual.linfty_norm();
750 
751  }
752  if(res/lastres > nonlinear_rho_)
753  {
754  build_matrix=true;
755  }
756  lastres=res;
757  res=newres;
758 
759  out<<"\t\t Newton step: " <<iter<<"\t Residual (rel.): "
760  << pde.GetOutputHandler()->ZeroTolerance(res/firstres, 1.0)
761  << "\t LineSearch {"<<lineiter<<"} ";
762 
763  pde.GetOutputHandler()->Write(out,priority);
764 
765  }//End of Linesearch
766  }
767  GetIntegrator().DeleteDomainData("last_time_solution");
768  GetIntegrator().DeleteDomainData("last_newton_solution");
769 
770 
771 
772  return build_matrix;
773  }
774 
775 
776 /*******************************************************************************************/
777  template <typename INTEGRATOR, typename LINEARSOLVER, typename VECTOR>
778  INTEGRATOR& FractionalStepThetaStepNewtonSolver<INTEGRATOR,LINEARSOLVER, VECTOR>
779  ::GetIntegrator()
780  {
781  return integrator_;
782  }
783 
784  /*******************************************************************************************/
785 
786 }
787 #endif
DOpE::DOpEIterationException
Definition: dopeexception.h:76
DOpE::FractionalStepThetaStepNewtonSolver::GetIntegrator
INTEGRATOR & GetIntegrator()
Definition: fractional_step_theta_step_newtonsolver.h:779
DOpE::ParameterReader::get_double
double get_double(const std::string &entry_name)
Definition: parameterreader.h:115
DOpE::ParameterReader
Definition: parameterreader.h:36
DOpE::FractionalStepThetaStepNewtonSolver::~FractionalStepThetaStepNewtonSolver
~FractionalStepThetaStepNewtonSolver()
Definition: fractional_step_theta_step_newtonsolver.h:225
DOpE::ParameterReader::declare_entry
void declare_entry(const std::string &entry, const std::string &default_value, const Patterns::PatternBase &pattern=Patterns::Anything(), const std::string &documentation=std::string())
Definition: parameterreader.h:98
DOpE::FractionalStepThetaStepNewtonSolver::FractionalStepThetaStepNewtonSolver
FractionalStepThetaStepNewtonSolver(INTEGRATOR &integrator, ParameterReader &param_reader)
Definition: fractional_step_theta_step_newtonsolver.h:207
DOpE::FractionalStepThetaStepNewtonSolver::NonlinearSolve_Initial
bool NonlinearSolve_Initial(PROBLEM &pde, VECTOR &solution, bool apply_boundary_values=true, bool force_matrix_build=false, int priority=5, std::string algo_level="\t\t ")
Definition: fractional_step_theta_step_newtonsolver.h:264
DOpE::FractionalStepThetaStepNewtonSolver::ReInit
void ReInit(PROBLEM &pde)
Definition: fractional_step_theta_step_newtonsolver.h:233
DOpE::FractionalStepThetaStepNewtonSolver
Definition: fractional_step_theta_step_newtonsolver.h:55
DOpE::FractionalStepThetaStepNewtonSolver::NonlinearLastTimeEvals
void NonlinearLastTimeEvals(PROBLEM &pde, const VECTOR &last_time_solution, VECTOR &residual)
Definition: fractional_step_theta_step_newtonsolver.h:242
DOpE::ParameterReader::get_integer
int get_integer(const std::string &entry_name)
Definition: parameterreader.h:126
DOpE::ParameterReader::SetSubsection
void SetSubsection(const std::string subsection)
Definition: parameterreader.h:93
DOpE::FractionalStepThetaStepNewtonSolver::declare_params
static void declare_params(ParameterReader &param_reader)
Definition: fractional_step_theta_step_newtonsolver.h:189
DOpE::FractionalStepThetaStepNewtonSolver::NonlinearSolve
bool NonlinearSolve(PROBLEM &pde, const VECTOR &last_time_solution, VECTOR &solution, bool apply_boundary_values=true, bool force_matrix_build=false, int priority=5, std::string algo_level="\t\t ")
Definition: fractional_step_theta_step_newtonsolver.h:388
Generated on Mon Jan 25 2016 11:32:19 for DOpE by   doxygen 1.8.6