1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
|
// ATC headers
#include "NonLinearSolver.h"
#include "LinearSolver.h"
#include "ATC_Error.h"
#include "LammpsInterface.h"
using std::stringstream;
namespace ATC {
//===================================================================
// TangentOperator
//===================================================================
//===================================================================
// NonLinearSolver
//===================================================================
NonLinearSolver::NonLinearSolver(TangentOperator * f,
const BC_SET * bcs, const int dof,
bool parallel):
f_(f),
bcs_(bcs),
dof_(dof),
rNorm0P_(1.0),
tol_(1.e-10),
tolx_(1.e-8),
tol0_(1.e-6),
maxIterations_(20),
parallel_(parallel)
{
}
//--------------------------------------------------------------------
double NonLinearSolver::residual_norm(VECTOR & r)
{
if (bcs_) {
DENS_VEC R = r;
BC_SET::const_iterator itr;
for (itr = bcs_->begin(); itr != bcs_->end(); itr++) {
int i = itr->first;
R(i) = 0;
}
return R.norm();
}
else { return r.norm(); }
}
//--------------------------------------------------------------------
bool NonLinearSolver::solve(VECTOR & x)
{
f_->function(x, r_);
rNorm0_ = residual_norm(r_);
if (rNorm0_ < tol_*rNorm0P_) { // if a "solution" does pass here rNorm0_ will be too small to allow for convergence
return true; // note abs vs rel tol
}
if (rNorm0_ == 0.0) rNorm0_ = 1.0;
if (rNorm0_ < tol0_ ) rNorm0_ = rNorm0P_;
if (rNorm0P_ == 1.0) rNorm0P_ = rNorm0_;
rNormP_ = rNorm0_;
dx_.reset(r_.nRows()); // needs to be sized for linear solver
// newton's method
for (int iter = 0; iter < maxIterations_ ; iter++ ) {
// compute tangent
f_->tangent(x, r_, A_);
rNorm_ = residual_norm(r_);
rNorm_ /= rNorm0_;
if (rNorm_ < tol_) {
return true;
}
SPAR_MAT Asparse(A_);
LinearSolver linearSolver(Asparse, LinearSolver::AUTO_SOLVE, parallel_);
if (bcs_) {
linearSolver.allow_reinitialization();
linearSolver.initialize(bcs_);
if (iter > 0) linearSolver.set_homogeneous_bcs();
else { x.zero(); } // linear solve w/ bcs will replace guess
}
r_ *= -1;
linearSolver.solve(dx_,r_);
if (iter > 0 && rNorm_ > rNormP_) {
bool descent = line_search(x);
if (! descent ) {
// return false;
}
}
rNormP_ = rNorm_;
x += dx_;
}
stringstream ss;
ss << "WARNING NonLinearSolver: did not converge, iterations="<< maxIterations_ <<" error= " << rNorm_;
ATC::LammpsInterface::instance()->print_msg_once(ss.str());
return false;
}
//--------------------------------------------------------------------
bool NonLinearSolver::line_search(VECTOR & x)
{
double rNormP = rNormP_;
double dxnorm = dx_.norm();
while ( dxnorm > tolx_) {
dx_ *= 0.5; // bisection
dxnorm = dx_.norm();
f_->function(x+dx_,r_);
rNorm_ = residual_norm(r_)/rNorm0_;
if (rNorm_ < rNormP) return true;
}
return false; // no descent
}
} // end namespace ATC
|
