24 #ifndef _NEWTON_SOLVER_H_
25 #define _NEWTON_SOLVER_H_
27 #include <lac/vector.h>
28 #include <lac/block_sparsity_pattern.h>
29 #include <lac/block_sparse_matrix.h>
31 #include <numerics/vector_tools.h>
53 template <
typename INTEGRATOR,
typename LINEARSOLVER,
typename VECTOR>
66 template<
typename PROBLEM>
96 template<
typename PROBLEM>
97 bool NonlinearSolve(PROBLEM& pde, VECTOR &solution,
bool apply_boundary_values=
true,
98 bool force_matrix_build=
false,
99 int priority = 5, std::string algo_level =
"\t\t ");
106 INTEGRATOR &_integrator;
110 double _nonlinear_global_tol, _nonlinear_tol, _nonlinear_rho;
111 double _linesearch_rho;
112 int _nonlinear_maxiter, _line_maxiter;
117 template <
typename INTEGRATOR,
typename LINEARSOLVER,
typename VECTOR>
122 param_reader.
declare_entry(
"nonlinear_global_tol",
"1.e-12",Patterns::Double(0),
"global tolerance for the newton iteration");
123 param_reader.
declare_entry(
"nonlinear_tol",
"1.e-10",Patterns::Double(0),
"relative tolerance for the newton iteration");
124 param_reader.
declare_entry(
"nonlinear_maxiter",
"10",Patterns::Integer(0),
"maximal number of newton iterations");
125 param_reader.
declare_entry(
"nonlinear_rho",
"0.1",Patterns::Double(0),
"minimal newton reduction, if actual reduction is less, matrix is rebuild ");
127 param_reader.
declare_entry(
"line_maxiter",
"4",Patterns::Integer(0),
"maximal number of linesearch steps");
128 param_reader.
declare_entry(
"linesearch_rho",
"0.9",Patterns::Double(0),
"reduction rate for the linesearch damping paramete");
130 LINEARSOLVER::declare_params(param_reader);
135 template <
typename INTEGRATOR,
typename LINEARSOLVER,
typename VECTOR>
138 : LINEARSOLVER(param_reader), _integrator(integrator)
141 _nonlinear_global_tol = param_reader.
get_double (
"nonlinear_global_tol");
142 _nonlinear_tol = param_reader.
get_double (
"nonlinear_tol");
143 _nonlinear_maxiter = param_reader.
get_integer (
"nonlinear_maxiter");
144 _nonlinear_rho = param_reader.
get_double (
"nonlinear_rho");
146 _line_maxiter = param_reader.
get_integer (
"line_maxiter");
147 _linesearch_rho = param_reader.
get_double (
"linesearch_rho");
153 template <
typename INTEGRATOR,
typename LINEARSOLVER,
typename VECTOR>
160 template <
typename INTEGRATOR,
typename LINEARSOLVER,
typename VECTOR>
161 template<
typename PROBLEM>
165 LINEARSOLVER::ReInit(pde);
169 template <
typename INTEGRATOR,
typename LINEARSOLVER,
typename VECTOR>
170 template<
typename PROBLEM>
174 bool apply_boundary_values,
175 bool force_matrix_build,
177 std::string algo_level)
179 bool build_matrix = force_matrix_build;
182 std::stringstream out;
183 pde.GetOutputHandler()->InitNewtonOut(out);
186 residual.reinit(solution);
188 if(apply_boundary_values)
190 GetIntegrator().ApplyInitialBoundaryValues(pde,solution);
192 pde.GetDoFConstraints().distribute(solution);
194 GetIntegrator().AddDomainData(
"last_newton_solution",&solution);
196 GetIntegrator().ComputeNonlinearResidual(pde,residual);
199 pde.GetOutputHandler()->SetIterationNumber(0,
"PDENewton");
200 pde.GetOutputHandler()->Write(residual,
"Residual"+pde.GetType(),pde.GetDoFType());
202 double res = residual.linfty_norm();
203 double firstres = res;
204 double lastres = res;
207 out<< algo_level <<
"Newton step: " <<0<<
"\t Residual (abs.): "
208 <<pde.GetOutputHandler()->ZeroTolerance(res, 1.0)
211 out<< algo_level <<
"Newton step: " <<0<<
"\t Residual (rel.): " << std::scientific << firstres/firstres;
214 pde.GetOutputHandler()->Write(out,priority);
217 while(res > _nonlinear_global_tol && res > firstres * _nonlinear_tol)
221 if(iter > _nonlinear_maxiter)
223 GetIntegrator().DeleteDomainData(
"last_newton_solution");
227 pde.GetOutputHandler()->SetIterationNumber(iter,
"PDENewton");
229 LINEARSOLVER::Solve(pde,GetIntegrator(),residual,du,build_matrix);
230 build_matrix =
false;
235 GetIntegrator().ComputeNonlinearResidual(pde,residual);
238 pde.GetOutputHandler()->Write(residual,
"Residual"+pde.GetType(),pde.GetDoFType());
239 pde.GetOutputHandler()->Write(du,
"Update"+pde.GetType(),pde.GetDoFType());
241 double newres = residual.linfty_norm();
243 double rho = _linesearch_rho;
248 out<< algo_level <<
"Newton step: " <<iter<<
"\t Residual (rel.): "
249 <<pde.GetOutputHandler()->ZeroTolerance(newres/firstres, 1.0)
250 <<
"\t LineSearch {"<<lineiter<<
"} ";
252 pde.GetOutputHandler()->Write(out,priority+1);
255 if(lineiter > _line_maxiter)
257 GetIntegrator().DeleteDomainData(
"last_newton_solution");
260 solution.add(alpha*(rho-1.),du);
263 GetIntegrator().ComputeNonlinearResidual(pde,residual);
265 pde.GetOutputHandler()->Write(residual,
"Residual"+pde.GetType(),pde.GetDoFType());
267 newres = residual.linfty_norm();
270 if(res/lastres > _nonlinear_rho)
280 <<
"\t Residual (rel.): "
281 << pde.GetOutputHandler()->ZeroTolerance(res/firstres, 1.0)
287 pde.GetOutputHandler()->Write(out,priority);
291 GetIntegrator().DeleteDomainData(
"last_newton_solution");
297 template <
typename INTEGRATOR,
typename LINEARSOLVER,
typename VECTOR>
Definition: dopeexception.h:76
double get_double(const std::string &entry_name)
Definition: parameterreader.h:115
void ReInit(PROBLEM &pde)
Definition: newtonsolver.h:163
Definition: parameterreader.h:36
INTEGRATOR & GetIntegrator()
Definition: newtonsolver.h:299
~NewtonSolver()
Definition: newtonsolver.h:155
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
Definition: newtonsolver.h:54
int get_integer(const std::string &entry_name)
Definition: parameterreader.h:126
NewtonSolver(INTEGRATOR &integrator, ParameterReader ¶m_reader)
Definition: newtonsolver.h:137
static void declare_params(ParameterReader ¶m_reader)
Definition: newtonsolver.h:119
void SetSubsection(const std::string subsection)
Definition: parameterreader.h:93
bool NonlinearSolve(PROBLEM &pde, VECTOR &solution, bool apply_boundary_values=true, bool force_matrix_build=false, int priority=5, std::string algo_level="\t\t ")
Definition: newtonsolver.h:172