24 #ifndef NEWTON_MIXED_SOLVER_H_
25 #define NEWTON_MIXED_SOLVER_H_
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>
55 template <
typename INTEGRATOR,
typename LINEARSOLVER,
typename VECTOR>
68 template<
typename PROBLEM>
99 template<
typename PROBLEM>
100 bool NonlinearSolve(PROBLEM& pde, VECTOR &solution,
bool apply_boundary_values=
true,
101 bool force_matrix_build=
false,
102 int priority = 5, std::string algo_level =
"\t\t ");
109 INTEGRATOR &integrator_;
113 double nonlinear_global_tol_, nonlinear_tol_, nonlinear_rho_;
114 double linesearch_rho_;
115 int nonlinear_maxiter_, line_maxiter_;
120 template <
typename INTEGRATOR,
typename LINEARSOLVER,
typename VECTOR>
125 param_reader.
declare_entry(
"nonlinear_global_tol",
"1.e-12",Patterns::Double(0),
"global tolerance for the newton iteration");
126 param_reader.
declare_entry(
"nonlinear_tol",
"1.e-10",Patterns::Double(0),
"relative tolerance for the newton iteration");
127 param_reader.
declare_entry(
"nonlinear_maxiter",
"10",Patterns::Integer(0),
"maximal number of newton iterations");
128 param_reader.
declare_entry(
"nonlinear_rho",
"0.1",Patterns::Double(0),
"minimal newton reduction, if actual reduction is less, matrix is rebuild ");
130 param_reader.
declare_entry(
"line_maxiter",
"4",Patterns::Integer(0),
"maximal number of linesearch steps");
131 param_reader.
declare_entry(
"linesearch_rho",
"0.9",Patterns::Double(0),
"reduction rate for the linesearch damping paramete");
133 LINEARSOLVER::declare_params(param_reader);
138 template <
typename INTEGRATOR,
typename LINEARSOLVER,
typename VECTOR>
141 : LINEARSOLVER(param_reader), integrator_(integrator)
144 nonlinear_global_tol_ = param_reader.
get_double (
"nonlinear_global_tol");
145 nonlinear_tol_ = param_reader.
get_double (
"nonlinear_tol");
146 nonlinear_maxiter_ = param_reader.
get_integer (
"nonlinear_maxiter");
147 nonlinear_rho_ = param_reader.
get_double (
"nonlinear_rho");
149 line_maxiter_ = param_reader.
get_integer (
"line_maxiter");
150 linesearch_rho_ = param_reader.
get_double (
"linesearch_rho");
156 template <
typename INTEGRATOR,
typename LINEARSOLVER,
typename VECTOR>
163 template <
typename INTEGRATOR,
typename LINEARSOLVER,
typename VECTOR>
164 template<
typename PROBLEM>
168 LINEARSOLVER::ReInit(pde);
172 template <
typename INTEGRATOR,
typename LINEARSOLVER,
typename VECTOR>
173 template<
typename PROBLEM>
177 bool apply_boundary_values,
178 bool force_matrix_build,
180 std::string algo_level)
182 bool build_matrix = force_matrix_build;
185 std::stringstream out;
186 pde.GetOutputHandler()->InitNewtonOut(out);
189 residual.reinit(solution);
191 if(apply_boundary_values)
193 GetIntegrator().ApplyInitialBoundaryValues(pde,solution);
196 dealii::Vector<double> copy_solution;
197 copy_solution= solution;
199 GetIntegrator().AddParamData(
"last_newton_solution",©_solution);
201 GetIntegrator().ComputeNonlinearResidual(pde,residual,
true);
204 pde.GetOutputHandler()->SetIterationNumber(0,
"PDENewton");
205 pde.GetOutputHandler()->Write(residual,
"Residual"+pde.GetType(),pde.GetDoFType());
207 double res = residual.linfty_norm();
208 double firstres = res;
209 double lastres = res;
211 out<<algo_level +
" Newton step: " <<0<<
"\t Residual (abs.): " <<pde.GetOutputHandler()->ZeroTolerance(res, 1.0)
214 out<<algo_level +
" Newton step: " <<0<<
"\t Residual (rel.): "
215 <<pde.GetOutputHandler()->ZeroTolerance(res/firstres, 1.0);
217 pde.GetOutputHandler()->Write(out,priority);
220 while(res > nonlinear_global_tol_ && res > firstres * nonlinear_tol_)
224 if(iter > nonlinear_maxiter_)
229 pde.GetOutputHandler()->SetIterationNumber(iter,
"PDENewton");
231 LINEARSOLVER::Solve(pde,GetIntegrator(),residual,du,build_matrix);
232 build_matrix =
false;
236 copy_solution = solution;
238 GetIntegrator().ComputeNonlinearResidual(pde,residual,
true);
240 pde.GetOutputHandler()->Write(residual,
"Residual"+pde.GetType(),pde.GetDoFType());
242 double newres = residual.linfty_norm();
244 double rho = linesearch_rho_;
249 out<<algo_level +
" Newton step: " <<iter<<
"\t Residual (rel.): "
250 << pde.GetOutputHandler()->ZeroTolerance(newres/firstres, 1.0)
251 <<
"\t LineSearch {"<<lineiter<<
"} ";
252 pde.GetOutputHandler()->Write(out,priority+1);
255 if(lineiter > line_maxiter_)
257 throw DOpEIterationException(
"Line-Iteration count exceeded bounds!",
"NewtonSolverMixedDims::NonlinearSolve");
259 solution.add(alpha*(rho-1.),du);
262 GetIntegrator().ComputeNonlinearResidual(pde,residual,
true);
264 pde.GetOutputHandler()->Write(residual,
"Residual"+pde.GetType(),pde.GetDoFType());
266 newres = residual.linfty_norm();
269 if(res/lastres > nonlinear_rho_)
276 out<<algo_level +
" Newton step: " <<iter<<
"\t Residual (rel.): "
277 << pde.GetOutputHandler()->ZeroTolerance(res/firstres, 1.0)
278 <<
"\t LineSearch {"<<lineiter<<
"} ";
279 pde.GetOutputHandler()->Write(out,priority);
282 copy_solution = solution;
284 GetIntegrator().DeleteParamData(
"last_newton_solution");
291 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: newtonsolvermixeddims.h:166
Definition: parameterreader.h:36
Definition: newtonsolvermixeddims.h:56
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
NewtonSolverMixedDimensions(INTEGRATOR &integrator, ParameterReader ¶m_reader)
Definition: newtonsolvermixeddims.h:140
INTEGRATOR & GetIntegrator()
Definition: newtonsolvermixeddims.h:293
int get_integer(const std::string &entry_name)
Definition: parameterreader.h:126
void SetSubsection(const std::string subsection)
Definition: parameterreader.h:93
static void declare_params(ParameterReader ¶m_reader)
Definition: newtonsolvermixeddims.h:122
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: newtonsolvermixeddims.h:175
~NewtonSolverMixedDimensions()
Definition: newtonsolvermixeddims.h:158