Newton.cpp 3.81 KB
Newer Older
1
#include "Newton.hpp"
2
3
4


void Newton::operator()(const IAssemble &A, Vector &u) {
5
6
7
8
9
10
11
12
  mout.StartBlock("Newton");
  double E, E_0;
  A.Initialize(u);
  vout(9) << "u(0)= " << endl << u << endl;
  Vector r(u);
  Matrix J(u);
  d = d_0 = A.Residual(u, r);
  E = E_0 = A.Energy(u);
13
  double eps = defaultEpsilon + defaultReduction * d;
14
15
  double d_previous = d;
  int LS_cnt = 0;
16

17
18
  int JU_cnt = 0; // counts iteration steps without Jacobian update
  Vector c(r);
19

20
21
22
23
  for (iter = 0; iter < max_iter; ++iter) {
    vout(3) << "r(" << iter << ")= " << endl << r << endl;
    if (d < eps) { break; }
    vout(1) << "d(" << iter << ")= " << d << endl;
24

25
26
27
28
    // Determines whether Jacobi Matrix is updated
    if (((iter - JU_cnt) >= JacobiUpdate) || (iter == 0)) {
      A.Jacobi(u, J);
      JU_cnt = iter;
29
30
      c = (*solver)(J) * r;
    } else c = (*solver)(J, 0) * r;
31

32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
    vout(3) << "c(" << iter << ")= " << endl << c << endl;
    u -= c;
    vout(5) << "u-c " << endl << u << endl;
    d = A.Residual(u, r);
    E = A.Energy(u);
    if (d > d_previous) {
      for (int l = 1; l <= LS_iter; ++l) {
        if (iter == 0)
          if (suppressLS) {
            vout(2) << "line search suppressed" << endl;
            break;
          }
        vout(1) << "line search " << l << ": d(" << iter << ")= " << d << endl;
        c *= 0.5;
        u += c;
47
48
        d = A.Residual(u, r);
        E = A.Energy(u);
49
50
        if (d < d_previous) break;
      }
51
    }
52
53
54
55
56
57
58
59
60
61
    if (d > d_previous) {
      vout(5) << "line search unsuccessful." << endl;
      ++LS_cnt;
      if (LS_cnt == 3) {
        vout(1) << "too many line searches unsuccessful." << endl;
        iter = max_iter;
      }
    }
    d_previous = d;
  }
62
  endOut(1);
63
64
}

65

66
void ENewton::operator()(const IAssemble &A, Vector &u) {
67
68
69
70
71
72
73
  mout.StartBlock("ENewton");
  double E, E_0;
  A.Initialize(u);
  Vector r(u);
  Matrix J(u);
  d = d_0 = A.Residual(u, r);
  E = E_0 = A.Energy(u);
74
  double eps = defaultEpsilon + defaultReduction * d;
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
  int LS_cnt = 0;
  int JU_cnt = 0;
  Vector c(r);
  for (iter = 0; iter < max_iter; ++iter) {
    vout(3) << "r(" << iter << ")= " << r << endl;
    if (d < eps) { break; }
    if (d > VeryLarge) { break; }
    if (!(E < VeryLarge)) {
      if (!(E > VeryLarge)) {
        iter = max_iter;
        break;
      }
    }
    vout(1) << "d(" << iter << ")= " << d << "\t E(" << iter << ")= " << E << endl;
    A.Jacobi(u, J);
    JU_cnt = iter;
91
    c = (*solver)(J) * r;
92
93
//    double rc = std::real(r * c);
    double rc = r * c;
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
    vout(3) << "c(" << iter << ")= " << c << endl;
    u -= c;
    vout(5) << "u-c " << u << endl;
    double d_old = d;
    d = A.Residual(u, r);
    double E_old = E;
    E = A.Energy(u);
    if (d < d_old) continue;
    if (E > E_old + sigma * rc) {
      for (int l = 1; l <= LS_iter; ++l) {
        if (iter == 0)
          if (suppressLS) {
            vout(2) << "line search suppressed" << endl;
            break;
          }
        vout(1) << "line search " << l
                << ": d(" << iter << ")= " << d
                << "\t E(" << iter << ")= " << E << endl;
        c *= 0.5;
        rc *= 0.5;
        u += c;
115
116
        d = A.Residual(u, r);
        E = A.Energy(u);
117
118
119
120
121
122
123
124
125
126
        if (E < E_old + sigma * rc) break;
      }
    }
    if (E > E_old + sigma * rc) {
      vout(5) << "line search unsuccessful." << endl;
      ++LS_cnt;
      if (LS_cnt == 3) {
        vout(1) << "too many line searches unsuccessful." << endl;
        iter = max_iter;
      }
127
    }
128
129
130
131
  }
  vout(0) << "d(" << iter << ")= " << d << "\t E(" << iter << ")= " << E
          << " rate " << rate() << endl;
  mout.EndBlock();
132
}
133

134
bool NewtonMethod(IAssemble &assemble, Vector &u) {
135
136
  Newton newton;
  newton.PrintInfo();
137
138
  assemble.PrintInfo();
  newton.operator()(assemble, u);
139
  return newton.converged();
140
}
141
142
143
144

bool NewtonMethod(IAssemble *assemble, Vector &u) {
  return NewtonMethod(*assemble, u);
}