Newton.cpp 3.81 KB
 1 ``````#include "Newton.hpp" `````` niklas.baumgarten committed Oct 12, 2020 2 3 4 `````` void Newton::operator()(const IAssemble &A, Vector &u) { `````` jonathan.wunderlich committed Oct 14, 2020 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); `````` niklas.baumgarten committed Mar 05, 2021 13 `````` double eps = defaultEpsilon + defaultReduction * d; `````` jonathan.wunderlich committed Oct 14, 2020 14 15 `````` double d_previous = d; int LS_cnt = 0; `````` niklas.baumgarten committed Oct 12, 2020 16 `````` `````` jonathan.wunderlich committed Oct 14, 2020 17 18 `````` int JU_cnt = 0; // counts iteration steps without Jacobian update Vector c(r); `````` niklas.baumgarten committed Oct 12, 2020 19 `````` `````` jonathan.wunderlich committed Oct 14, 2020 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; `````` niklas.baumgarten committed Oct 12, 2020 24 `````` `````` jonathan.wunderlich committed Oct 14, 2020 25 26 27 28 `````` // Determines whether Jacobi Matrix is updated if (((iter - JU_cnt) >= JacobiUpdate) || (iter == 0)) { A.Jacobi(u, J); JU_cnt = iter; `````` niklas.baumgarten committed Feb 25, 2021 29 30 `````` c = (*solver)(J) * r; } else c = (*solver)(J, 0) * r; `````` niklas.baumgarten committed Oct 12, 2020 31 `````` `````` jonathan.wunderlich committed Oct 14, 2020 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; `````` niklas.baumgarten committed Oct 12, 2020 47 48 `````` d = A.Residual(u, r); E = A.Energy(u); `````` jonathan.wunderlich committed Oct 14, 2020 49 50 `````` if (d < d_previous) break; } `````` niklas.baumgarten committed Oct 12, 2020 51 `````` } `````` jonathan.wunderlich committed Oct 14, 2020 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; } `````` niklas.baumgarten committed Mar 15, 2021 62 `````` endOut(1); `````` niklas.baumgarten committed Oct 12, 2020 63 64 ``````} `````` niklas.baumgarten committed Feb 25, 2021 65 `````` `````` niklas.baumgarten committed Oct 12, 2020 66 ``````void ENewton::operator()(const IAssemble &A, Vector &u) { `````` jonathan.wunderlich committed Oct 14, 2020 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); `````` niklas.baumgarten committed Mar 05, 2021 74 `````` double eps = defaultEpsilon + defaultReduction * d; `````` jonathan.wunderlich committed Oct 14, 2020 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; `````` niklas.baumgarten committed Feb 25, 2021 91 `````` c = (*solver)(J) * r; `````` jonathan.wunderlich committed Mar 13, 2021 92 93 ``````// double rc = std::real(r * c); double rc = r * c; `````` jonathan.wunderlich committed Oct 14, 2020 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; `````` niklas.baumgarten committed Oct 12, 2020 115 116 `````` d = A.Residual(u, r); E = A.Energy(u); `````` jonathan.wunderlich committed Oct 14, 2020 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; } `````` niklas.baumgarten committed Oct 12, 2020 127 `````` } `````` jonathan.wunderlich committed Oct 14, 2020 128 129 130 131 `````` } vout(0) << "d(" << iter << ")= " << d << "\t E(" << iter << ")= " << E << " rate " << rate() << endl; mout.EndBlock(); `````` niklas.baumgarten committed Oct 12, 2020 132 ``````} `````` niklas.baumgarten committed Mar 01, 2021 133 `````` `````` jonathan.wunderlich committed Apr 12, 2021 134 ``````bool NewtonMethod(IAssemble &assemble, Vector &u) { `````` niklas.baumgarten committed Mar 05, 2021 135 136 `````` Newton newton; newton.PrintInfo(); `````` jonathan.wunderlich committed Apr 12, 2021 137 138 `````` assemble.PrintInfo(); newton.operator()(assemble, u); `````` niklas.baumgarten committed Mar 05, 2021 139 `````` return newton.converged(); `````` niklas.baumgarten committed Mar 01, 2021 140 ``````} `````` jonathan.wunderlich committed Apr 12, 2021 141 142 143 144 `````` bool NewtonMethod(IAssemble *assemble, Vector &u) { return NewtonMethod(*assemble, u); }``````