Reference documentation for deal.II version 8.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
newton.templates.h
1 // ---------------------------------------------------------------------
2 // @f$Id: newton.templates.h 31349 2013-10-20 19:07:06Z maier @f$
3 //
4 // Copyright (C) 2006 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 
18 #include <deal.II/algorithms/newton.h>
19 
20 #include <deal.II/base/parameter_handler.h>
21 #include <deal.II/base/logstream.h>
22 #include <deal.II/lac/vector_memory.h>
23 
24 #include <iomanip>
25 
26 
27 DEAL_II_NAMESPACE_OPEN
28 
29 namespace Algorithms
30 {
31  template <class VECTOR>
33  :
34  residual(&residual), inverse_derivative(&inverse_derivative),
35  assemble_now(false),
36  n_stepsize_iterations(21),
37  assemble_threshold(0.),
38  debug_vectors(false),
39  debug(0)
40  {}
41 
42 
43  template <class VECTOR>
44  void
46  {
47  param.enter_subsection("Newton");
49  param.declare_entry("Assemble threshold", "0.", Patterns::Double());
50  param.declare_entry("Stepsize iterations", "21", Patterns::Integer());
51  param.declare_entry("Debug level", "0", Patterns::Integer());
52  param.declare_entry("Debug vectors", "false", Patterns::Bool());
53  param.leave_subsection();
54  }
55 
56  template <class VECTOR>
57  void
59  {
60  param.enter_subsection("Newton");
61  control.parse_parameters (param);
62  assemble_threshold = param.get_double("Assemble threshold");
63  n_stepsize_iterations = param.get_integer("Stepsize iterations");
64  debug_vectors = param.get_bool("Debug vectors");
65  param.leave_subsection ();
66  }
67 
68  template <class VECTOR>
69  void
71  {
72  data_out = &output;
73  }
74 
75  template <class VECTOR>
76  void
78  {
79  residual->notify(e);
80  inverse_derivative->notify(e);
81  }
82 
83 
84  template <class VECTOR>
85  void
87  {
88  Assert (out.size() == 1, ExcNotImplemented());
89  deallog.push ("Newton");
90 
91  VECTOR &u = *out(0);
92 
93  if (debug>2)
94  deallog << "u: " << u.l2_norm() << std::endl;
95 
97  typename VectorMemory<VECTOR>::Pointer Du(mem);
98  typename VectorMemory<VECTOR>::Pointer res(mem);
99 
100  res->reinit(u);
101  NamedData<VECTOR *> src1;
102  NamedData<VECTOR *> src2;
103  VECTOR *p = &u;
104  src1.add(p, "Newton iterate");
105  src1.merge(in);
106  p = res;
107  src2.add(p, "Newton residual");
108  src2.merge(src1);
109  NamedData<VECTOR *> out1;
110  out1.add(p, "Residual");
111  p = Du;
112  NamedData<VECTOR *> out2;
113  out2.add(p, "Update");
114 
115  unsigned int step = 0;
116  // fill res with (f(u), v)
117  (*residual)(out1, src1);
118  double resnorm = res->l2_norm();
119  double old_residual = resnorm / assemble_threshold + 1;
120 
121  while (control.check(step++, resnorm) == SolverControl::iterate)
122  {
123  // assemble (Df(u), v)
124  if (resnorm/old_residual >= assemble_threshold)
125  inverse_derivative->notify (Events::bad_derivative);
126 
127  Du->reinit(u);
128  try
129  {
130  (*inverse_derivative)(out2, src2);
131  }
133  {
134  deallog << "Inner iteration failed after "
135  << e.last_step << " steps with residual "
136  << e.last_residual << std::endl;
137  }
138 
139  if (debug_vectors)
140  {
142  VECTOR *p = &u;
143  out.add(p, "solution");
144  p = Du;
145  out.add(p, "update");
146  p = res;
147  out.add(p, "residual");
148  *data_out << step;
149  *data_out << out;
150  }
151 
152  u.add(-1., *Du);
153  old_residual = resnorm;
154  (*residual)(out1, src1);
155  resnorm = res->l2_norm();
156 
157  // Step size control
158  unsigned int step_size = 0;
159  while (resnorm >= old_residual)
160  {
161  ++step_size;
162  if (step_size > n_stepsize_iterations)
163  {
164  deallog << "No smaller stepsize allowed!";
165  break;
166  }
167  if (control.log_history())
168  deallog << "Trying step size: 1/" << (1<<step_size)
169  << " since residual was " << resnorm << std::endl;
170  u.add(1./(1<<step_size), *Du);
171  (*residual)(out1, src1);
172  resnorm = res->l2_norm();
173  }
174  }
175  deallog.pop();
176 
177  // in case of failure: throw exception
178  if (control.last_check() != SolverControl::success)
179  AssertThrow(false, SolverControl::NoConvergence (control.last_step(),
180  control.last_value()));
181  // otherwise exit as normal
182  }
183 }
184 
185 
186 DEAL_II_NAMESPACE_CLOSE
void pop()
static void declare_parameters(ParameterHandler &param)
void add(DATA &v, const std::string &name)
Definition: named_data.h:253
Continue iteration.
void initialize(ParameterHandler &param)
#define AssertThrow(cond, exc)
Definition: exceptions.h:362
virtual void operator()(NamedData< VECTOR * > &out, const NamedData< VECTOR * > &in)
void enter_subsection(const std::string &subsection)
Newton(Operator< VECTOR > &residual, Operator< VECTOR > &inverse_derivative)
const unsigned int last_step
Stop iteration, goal reached.
#define Assert(cond, exc)
Definition: exceptions.h:299
const Event bad_derivative
unsigned int size() const
Number of stored data objects.
Definition: named_data.h:325
void declare_parameters(ParameterHandler &param)
void merge(NamedData< DATA2 > &)
Definition: named_data.h:278
bool get_bool(const std::string &entry_name) const
void push(const std::string &text)
virtual void notify(const Event &)
bool leave_subsection()
double get_double(const std::string &entry_name) const
void declare_entry(const std::string &entry, const std::string &default_value, const Patterns::PatternBase &pattern=Patterns::Anything(), const std::string &documentation=std::string())
::ExceptionBase & ExcNotImplemented()
long int get_integer(const std::string &entry_string) const