Newton.hpp 4.33 KB
Newer Older
1
2
3
#ifndef _NEWTON_H_
#define _NEWTON_H_

4
#include "LinearSolver.hpp"
5
#include "Assemble.hpp"
6

7

8
class NonLinearSolver {
9
10
11
12
13
14
15
16
17
protected:
  int verbose = 1;

  int LS_iter = 3;

  int min_iter = 0;

  int max_iter = 10;

18
  double defaultEpsilon = 1e-10;
19

20
  double defaultReduction = 1e-5;
21
22
23

  std::string prefix = "NonLinear";

24
  int iter; // todo remove
25

26
  double d, d_0; // todo remove
27

28
  std::unique_ptr<LinearSolver> solver;
29

30
public:
31
  NonLinearSolver(std::unique_ptr<LinearSolver> &&solver, const string &prefix = "NonLinear") :
32
      solver(std::move(solver)) {
33
34
35
    config.get(prefix + "Verbose", verbose);
    config.get(prefix + "Steps", max_iter);
    config.get(prefix + "LineSearchSteps", LS_iter);
36
37
    config.get(prefix + "Epsilon", defaultEpsilon);
    config.get(prefix + "Reduction", defaultReduction);
38
39
40
    config.get(prefix + "MinimalStepNumber", min_iter);
  }

41
42
  void PrintInfo() const {
    mout.PrintInfo("Solver", verbose,
43
44
45
46
                   PrintInfoEntry(prefix + " Solver", this->Name(), 0),
                   PrintInfoEntry(prefix + " Steps", max_iter, 1),
                   PrintInfoEntry(prefix + " Epsilon", defaultEpsilon, 1),
                   PrintInfoEntry(prefix + " Reduction", defaultReduction, 1),
47
                   PrintInfoEntry(this->solver->GetLinearPrefix() + " Solver",
48
                                  this->solver->Name(), 0),
49
                   PrintInfoEntry(this->solver->GetLinearPrefix() + " Preconditioner",
50
                                  this->solver->GetPreconditioner().Name(), 0),
51
                   PrintInfoEntry(this->solver->GetLinearPrefix() + " Steps",
52
                                  this->solver->GetLinearMaxIter(), 1),
53
                   PrintInfoEntry(this->solver->GetLinearPrefix() + " Epsilon",
54
                                  this->solver->GetLinearDefaultEps(), 1),
55
                   PrintInfoEntry(this->solver->GetLinearPrefix() + " Reduction",
56
                                  this->solver->GetLinearDefaultReduction(), 1)
57
58
59
    );
  }

60
61
62
63
  virtual void operator()(const IAssemble &, Vector &) = 0;

  virtual std::string Name() const = 0;

64
  ~NonLinearSolver() {}
65

66
  // Todo make converged return value of operator
67
68
69
70
71
72
73
74
75
  bool converged() const { return (iter < max_iter); }

  double rate() const {
    if (abs(d) > VeryLarge) Exit("no convergence in iteration");
    if (iter) if (d_0 != 0) return pow(d / d_0, 1.0 / iter);
    return 0;
  }

  int Iter() const { return iter; }
76

77
  double defect() const { return d; }
78

79
80
81
82
83
  void endOut(int v) {
    vout(v) << "d(" << iter << ")= " << d << " rate " << rate() << endl;
    mout.EndBlock(verbose, v);
  }

84
  friend Logging &operator<<(Logging &s, const NonLinearSolver &S) {
85
86
    return s << "d(" << S.iter << ")= " << S.d << " rate " << S.rate() << endl;
  }
87
88
};

89
class Newton : public NonLinearSolver {
90
protected:
91
  int suppressLS = 0;
92

93
  int JacobiUpdate = 1;
94

95
public:
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
  /*
   *  Gives Newton's method with GMRES as linear solver
   *  preconditioned with an incomplete LU decomposition,
   *  if "LinearSolver" or "Preconditioner" are not found in config.
   */
  Newton() :
      NonLinearSolver(std::unique_ptr<LinearSolver>(GetLinearSolver()), "Newton") {
    config.get("NewtonSuppressFirstLineSearch", suppressLS);
    config.get("NewtonJacobiUpdate", JacobiUpdate);
  }

  [[deprecated]] Newton(LinearSolver &linearSolver) :
      NonLinearSolver(std::unique_ptr<LinearSolver>(std::move(&linearSolver)), "Newton") {
    config.get("NewtonSuppressFirstLineSearch", suppressLS);
    config.get("NewtonJacobiUpdate", JacobiUpdate);
  }

  Newton(std::unique_ptr<LinearSolver> &&linearSolver) :
      NonLinearSolver(std::move(linearSolver), "Newton") {
    config.get("NewtonSuppressFirstLineSearch", suppressLS);
    config.get("NewtonJacobiUpdate", JacobiUpdate);
  }

  void operator()(const IAssemble &, Vector &) override;

  std::string Name() const override { return "Newton"; }
122
123
};

124
125
126
bool NewtonMethod(IAssemble &assemble, Vector &u);

/// Note that the assemble is not deleted here!
127
bool NewtonMethod(IAssemble *assemble, Vector &u);
128

129
class ENewton : public Newton {
130
  double sigma = 0.25;
131
public:
132
133
134
135
136
137
  ENewton(LinearSolver &linearSolver) :
      Newton(std::unique_ptr<LinearSolver>(std::move(&linearSolver))) {
    config.get("NewtonSigma", sigma);
  }

  void operator()(const IAssemble &, Vector &) override;
138

139
  std::string Name() const override { return "ENewton"; }
140
141
142
};

#endif // of #ifndef _NEWTON_H_