MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
nr.cc
Go to the documentation of this file.
1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/base/nr.cc,v 1.66 2017/01/12 14:46:10 masarati Exp $ */
2 /*
3  * MBDyn (C) is a multibody analysis code.
4  * http://www.mbdyn.org
5  *
6  * Copyright (C) 1996-2017
7  *
8  * Pierangelo Masarati <masarati@aero.polimi.it>
9  * Paolo Mantegazza <mantegazza@aero.polimi.it>
10  *
11  * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano
12  * via La Masa, 34 - 20156 Milano, Italy
13  * http://www.aero.polimi.it
14  *
15  * Changing this copyright notice is forbidden.
16  *
17  * This program is free software; you can redistribute it and/or modify
18  * it under the terms of the GNU General Public License as published by
19  * the Free Software Foundation (version 2 of the License).
20  *
21  *
22  * This program is distributed in the hope that it will be useful,
23  * but WITHOUT ANY WARRANTY; without even the implied warranty of
24  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25  * GNU General Public License for more details.
26  *
27  * You should have received a copy of the GNU General Public License
28  * along with this program; if not, write to the Free Software
29  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
30  */
31 
32  /*
33  *
34  * Copyright (C) 2003-2017
35  * Giuseppe Quaranta <quaranta@aero.polimi.it>
36  *
37  * classi che implementano la risoluzione del sistema nonlineare
38  */
39 
40 #include "mbconfig.h" /* This goes first in every *.c,*.cc file */
41 
42 #include <unistd.h>
43 #include "clock_time.h"
44 #include "solver.h"
45 #include "nr.h"
46 #ifdef USE_MPI
47 #include "mbcomm.h"
48 #include "schsolman.h"
49 #endif /* USE_MPI */
50 
51 #include "dofown.h"
52 #include "output.h"
53 
54 #include <cfloat>
55 #include <cmath>
56 #include <cstdlib>
57 
59  const bool bKJ,
60  const integer IterBfAss,
62 : NonlinearSolver(options), pRes(NULL),
63 pSol(NULL),
64 bTrueNewtonRaphson(bTNR),
65 IterationBeforeAssembly(IterBfAss),
66 bKeepJac(bKJ),
67 iPerformedIterations(0),
68 pPrevNLP(0)
69 {
70  NO_OP;
71 }
72 
74 {
75  NO_OP;
76 }
77 
78 void
80  Solver *pS,
81  const integer iMaxIter,
82  const doublereal& Tol,
83  integer& iIterCnt,
84  doublereal& dErr,
85  const doublereal& SolTol,
86  doublereal& dSolErr)
87 {
88  ASSERT(pS != NULL);
90 
91  iIterCnt = 0;
92  if ((!bKeepJac) || (pNLP != pPrevNLP)) {
94  }
95  pPrevNLP = pNLP;
96  dSolErr = 0.;
97 
98  doublereal dOldErr;
99  doublereal dErrFactor = 1.;
100  doublereal dErrDiff = 0.;
101  bool bJacBuilt = false;
102  doublereal dStartTimeCPU, dEndTimeCPU, dJacobianCPU, dResidualCPU = 0., dLinSolveCPU;
103 
104  while (true) {
105  if (outputCPUTime()) {
106  dStartTimeCPU = mbdyn_clock_time();
107  }
108 
109  pRes = pSM->pResHdl();
110  pSol = pSM->pSolHdl();
111  Size = pRes->iGetSize();
112 
113 #ifdef USE_EXTERNAL
114  SendExternal();
115 #endif /* USE_EXTERNAL */
116 
117  pRes->Reset();
118  bool forceJacobian(false);
119  try {
120  pNLP->Residual(pRes);
121  }
123  if (bHonorJacRequest) {
124  forceJacobian = true;
125  }
126  }
127 
128  /* FIXME: if Tol == 0., no convergence on residual
129  * is required, so we could simply don't compute
130  * the test; I'm leaving it in place so it appears
131  * in the output (maybe we could conditionally disable
132  * it?) */
133 
134  bool bTest = MakeResTest(pS, pNLP, *pRes, Tol, dErr, dErrDiff);
135 
136  if (outputCPUTime()) {
137  dEndTimeCPU = mbdyn_clock_time();
138  dResidualCPU += dEndTimeCPU - dStartTimeCPU;
139  AddTimeCPU(dResidualCPU, CPU_RESIDUAL);
140  }
141 
142  if (outputRes()) {
143  pS->PrintResidual(*pRes, iIterCnt);
144  }
145 
146  if (iIterCnt > 0) {
147  dErrFactor *= dErr/dOldErr;
148  }
149  dOldErr = dErr;
150 
151 #ifdef USE_MPI
152  if (!bParallel || MBDynComm.Get_rank() == 0)
153 #endif /* USE_MPI */
154  {
156  if (outputIters()) {
157  silent_cout("\tIteration(" << iIterCnt << ") " << dErr);
158 
159  if (bJacBuilt) {
160  silent_cout(" J");
161  }
162  }
163 
165  silent_cout(" cond=");
166  doublereal dCond;
167  if (pSM->bGetConditionNumber(dCond)) {
168  silent_cout(dCond);
170  AddCond(dCond);
171  silent_cout(" " << dGetCondMin() << " " << dGetCondMax() << " " << dGetCondAvg());
172  }
173  } else {
174  silent_cout("NA");
175  }
176  }
177 
178  if (outputCPUTime()) {
179  silent_cout(" CPU:" << dResidualCPU << "/" << dGetTimeCPU(CPU_RESIDUAL)
180  << "+" << dJacobianCPU << "/" << dGetTimeCPU(CPU_JACOBIAN)
181  << "+" << dLinSolveCPU << "/" << dGetTimeCPU(CPU_LINEAR_SOLVER));
182  }
183 
184  silent_cout(std::endl);
185  }
186  }
187 
188  pS->CheckTimeStepLimit(dErr, dErrDiff);
189 
190  if (bTest) {
191  return;
192  }
193 
194  if (iIterCnt >= std::abs(iMaxIter)) {
195  if (iMaxIter < 0 && dErrFactor < 1.) {
196  return;
197  }
198  if (outputBailout()) {
199  pS->PrintResidual(*pRes, iIterCnt);
200  }
202  }
203 
204  iIterCnt++;
205  bJacBuilt = false;
206 
207  if (outputCPUTime()) {
208  dStartTimeCPU = mbdyn_clock_time();
209  }
210 
213  || forceJacobian)
214  {
215  pSM->MatrReset();
216 rebuild_matrix:;
217  try {
218  pNLP->Jacobian(pSM->pMatHdl());
220  silent_cout("NewtonRaphsonSolver: "
221  "rebuilding matrix..."
222  << std::endl);
223 
224  /* need to rebuild the matrix... */
225  pSM->MatrInitialize();
226  /* FIXME: could loop forever! */
227  goto rebuild_matrix;
228 
229  } catch (...) {
230  throw;
231  }
232 
233  TotJac++;
234  bJacBuilt = true;
235  }
236 
237  if (outputCPUTime()) {
238  dEndTimeCPU = mbdyn_clock_time();
239  dJacobianCPU = dEndTimeCPU - dStartTimeCPU;
240  AddTimeCPU(dJacobianCPU, CPU_JACOBIAN);
241  }
242 
244 
245 #ifdef USE_MPI
246  if (!bParallel || MBDynComm.Get_rank() == 0)
247 #endif /* USE_MPI */
248  {
249  if (outputJac()) {
250  if (bJacBuilt) {
251  silent_cout("Jacobian:" << std::endl
252  << *(pSM->pMatHdl()));
253  } else {
254  silent_cout("Jacobian: unchanged" << std::endl);
255  }
256  }
257 
258  if (bJacBuilt && outputMatrixConditionNumber()) {
259  silent_cout("cond=" << pSM->pMatHdl()->ConditionNumber(GetCondMatNorm()) << std::endl);
260  }
261  }
262 
263  if (outputCPUTime()) {
264  dStartTimeCPU = mbdyn_clock_time();
265  }
266 
267  pSM->Solve();
268 
269  if (outputCPUTime()) {
270  dEndTimeCPU = mbdyn_clock_time();
271  dLinSolveCPU = dEndTimeCPU - dStartTimeCPU;
272  AddTimeCPU(dLinSolveCPU, CPU_LINEAR_SOLVER);
273  }
274 
275  if (outputSol()) {
276  pS->PrintSolution(*pSol, iIterCnt);
277  }
278 
279  if (outputCPUTime()) {
280  dStartTimeCPU = mbdyn_clock_time();
281  }
282 
283  pNLP->Update(pSol);
284 
285  bTest = MakeSolTest(pS, *pSol, SolTol, dSolErr);
286 
287  if (outputCPUTime()) {
288  dEndTimeCPU = mbdyn_clock_time();
289  dResidualCPU = dEndTimeCPU - dStartTimeCPU;
290  }
291 
292  if (outputIters()) {
293 #ifdef USE_MPI
294  if (!bParallel || MBDynComm.Get_rank() == 0)
295 #endif /* USE_MPI */
296  {
297  silent_cout("\t\tSolErr "
298  << dSolErr << std::endl);
299  }
300  }
301 
302  if (bTest) {
304  }
305 
306  // allow to bail out in case of multiple CTRL^C
309  }
310  }
311 }
virtual void Reset(void)=0
virtual VectorHandler * pResHdl(void) const =0
void Solve(const NonlinearProblem *pNLP, Solver *pS, const integer iMaxIter, const doublereal &Tol, integer &iIterCnt, doublereal &dErr, const doublereal &SolTol, doublereal &dSolErr)
Definition: nr.cc:79
bool outputMatrixConditionNumber(void) const
virtual SolutionManager * pGetSolutionManager(void) const
Definition: solver.h:398
void AddCond(doublereal dCond)
Definition: nonlin.h:331
virtual void Jacobian(MatrixHandler *pJac) const =0
doublereal dGetCondMax() const
Definition: nonlin.h:275
bool outputIters(void) const
integer TotJac
Definition: nonlin.h:252
#define MBDYN_EXCEPT_ARGS
Definition: except.h:63
bool outputBailout(void) const
NewtonRaphsonSolver(const bool bTNR, const bool bKJ, const integer IterBfAss, const NonlinearSolverOptions &options)
Definition: nr.cc:58
MatrixHandler::Norm_t GetCondMatNorm(void) const
virtual bool MakeSolTest(Solver *pS, const VectorHandler &Vec, const doublereal &dTol, doublereal &dTest)
Definition: nonlin.cc:497
virtual void PrintResidual(const VectorHandler &Res, integer iIterCnt) const
Definition: solver.cc:5194
virtual bool MakeResTest(Solver *pS, const NonlinearProblem *pNLP, const VectorHandler &Vec, const doublereal &dTol, doublereal &dTest, doublereal &dTestDiff)
Definition: nonlin.cc:474
double mbdyn_clock_time()
Definition: clock_time.cc:51
virtual void CheckTimeStepLimit(doublereal dErr, doublereal dErrDiff) const
Definition: solver.cc:5205
doublereal dGetCondAvg() const
Definition: nonlin.h:277
bool outputRes(void) const
integer Size
Definition: nonlin.h:251
virtual void Residual(VectorHandler *pRes) const =0
bool outputJac(void) const
#define NO_OP
Definition: myassert.h:74
doublereal dGetCondMin() const
Definition: nonlin.h:276
virtual integer iGetSize(void) const =0
doublereal dGetTimeCPU(CPUTimeType eType) const
Definition: nonlin.h:346
~NewtonRaphsonSolver(void)
Definition: nr.cc:73
bool outputSolverConditionNumber(void) const
virtual doublereal ConditionNumber(enum Norm_t eNorm=NORM_1) const
Definition: mh.cc:451
bool bKeepJac
Definition: nr.h:52
int mbdyn_stop_at_end_of_iteration(void)
Definition: solver.cc:172
bool outputSol(void) const
VectorHandler * pRes
Definition: nr.h:48
virtual MatrixHandler * pMatHdl(void) const =0
integer iPerformedIterations
Definition: nr.h:53
const NonlinearProblem * pPrevNLP
Definition: nr.h:54
virtual void MatrReset(void)=0
void AddTimeCPU(doublereal dTime, CPUTimeType eType)
Definition: nonlin.h:355
integer IterationBeforeAssembly
Definition: nr.h:51
virtual void Solve(void)=0
#define ASSERT(expression)
Definition: colamd.c:977
virtual void MatrInitialize(void)
Definition: solman.cc:72
struct option options[]
Definition: ann_in.c:46
Definition: solver.h:78
bool bTrueNewtonRaphson
Definition: nr.h:50
bool bGetConditionNumber(doublereal &dCond) const
Definition: solman.cc:107
virtual void PrintSolution(const VectorHandler &Sol, integer iIterCnt) const
Definition: solver.cc:5200
bool outputCPUTime(void) const
VectorHandler * pSol
Definition: nr.h:49
virtual VectorHandler * pSolHdl(void) const =0
virtual void Update(const VectorHandler *pSol) const =0
double doublereal
Definition: colamd.c:52
long int integer
Definition: colamd.c:51
bool outputSolverConditionStat(void) const