1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/base/solver.cc,v 1.301 2017/10/14 23:58:06 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 * Copyright 1999-2017 Giuseppe Quaranta <quaranta@aero.polimi.it>
34 * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano
35 *
36 * This copyright statement applies to the MPI related code, which was
37 * merged from files schur.h/schur.cc
38 */
39
40 /*
41 *
42 * Copyright (C) 2003-2017
43 * Giuseppe Quaranta <quaranta@aero.polimi.it>
44 *
45 */
46
47 /* metodo per la soluzione del modello */
48
49 #include "mbconfig.h" /* This goes first in every *.c,*.cc file */
50
51 /* required for configure time macros with paths */
52 #include "mbdefs.h"
53
54 #include <cstdlib>
55 #include <cstring>
56 #include <limits>
57 #include <ac/unistd.h>
58 #include <cerrno>
59 #include <csignal>
60 #include <cfloat>
61 #include <cmath>
62 #include <vector>
63 #include <algorithm>
64 extern "C" { int get_nprocs_conf(void); int get_nprocs(void); }
65
66 #include "solver.h"
67 #include "dataman.h"
68 #include "mtdataman.h"
69 #include "thirdorderstepsol.h"
70 #include "nr.h"
71 #include "linesearch.h"
72 #include "bicg.h"
73 #include "gmres.h"
74 #include "solman.h"
75 #include "readlinsol.h"
76 #include "ls.h"
77 #include "naivewrap.h"
78 #include "Rot.hh"
79 #include "cleanup.h"
80 #include "drive_.h"
81 #include "TimeStepControl.h"
82 #include "solver_impl.h"
83
84 #include "ac/lapack.h"
85 #include "ac/arpack.h"
86 #include "eigjdqz.h"
87
88 extern "C" {int get_nprocs(void);}
89
90 #ifdef HAVE_SIGNAL
91 /*
92 * MBDyn starts with mbdyn_keep_going == MBDYN_KEEP_GOING.
93 *
94 * A single CTRL^C sets it to MBDYN_STOP_AT_END_OF_TIME_STEP,
95 * which results in exiting at the end of the time step,
96 * after the output in case of success.
97 *
98 * A second CTRL^C sets it to MBDYN_STOP_AT_END_OF_ITERATION,
99 * which results in exiting at the end of the current iteration,
100 * after printing debug output if required.
101 *
102 * A further CTRL^C sets it to MBDYN_STOP_NOW and in throwing
103 * an exception.
104 */
105 enum {
106 MBDYN_KEEP_GOING = 0,
107 MBDYN_STOP_AT_END_OF_TIME_STEP = 1,
108 MBDYN_STOP_AT_END_OF_ITERATION = 2,
109 MBDYN_STOP_NOW = 3
110
111 };
112
113 volatile sig_atomic_t mbdyn_keep_going = MBDYN_KEEP_GOING;
114 #if defined(__FreeBSD__) || defined(__DragonFly__)
115 __sighandler_t *mbdyn_sh_term = SIG_DFL;
116 __sighandler_t *mbdyn_sh_int = SIG_DFL;
117 __sighandler_t *mbdyn_sh_hup = SIG_DFL;
118 __sighandler_t *mbdyn_sh_pipe = SIG_DFL;
119 #else // !defined(__FreeBSD__)
120 __sighandler_t mbdyn_sh_term = SIG_DFL;
121 __sighandler_t mbdyn_sh_int = SIG_DFL;
122 __sighandler_t mbdyn_sh_hup = SIG_DFL;
123 __sighandler_t mbdyn_sh_pipe = SIG_DFL;
124 #endif // !defined(__FreeBSD__)
125
126 extern "C" void
mbdyn_really_exit_handler(int signum)127 mbdyn_really_exit_handler(int signum)
128 {
129 ::mbdyn_keep_going = MBDYN_STOP_NOW;
130 switch (signum) {
131 case SIGTERM:
132 signal(signum, ::mbdyn_sh_term);
133 break;
134
135 case SIGINT:
136 signal(signum, ::mbdyn_sh_int);
137 break;
138
139 #ifdef SIGHUP
140 case SIGHUP:
141 signal(signum, ::mbdyn_sh_hup);
142 break;
143 #endif // SIGHUP
144
145 #ifdef SIGPIPE
146 case SIGPIPE:
147 signal(signum, ::mbdyn_sh_pipe);
148 break;
149 #endif // SIGPIPE
150 }
151
152 mbdyn_cleanup();
153
154 throw ErrInterrupted(MBDYN_EXCEPT_ARGS);
155 }
156
157 extern "C" void
mbdyn_modify_last_iteration_handler(int signum)158 mbdyn_modify_last_iteration_handler(int signum)
159 {
160 ::mbdyn_keep_going = MBDYN_STOP_AT_END_OF_ITERATION;
161 signal(signum, mbdyn_really_exit_handler);
162 }
163
164 extern "C" void
mbdyn_modify_final_time_handler(int signum)165 mbdyn_modify_final_time_handler(int signum)
166 {
167 ::mbdyn_keep_going = MBDYN_STOP_AT_END_OF_TIME_STEP;
168 signal(signum, mbdyn_modify_last_iteration_handler);
169 }
170 #endif /* HAVE_SIGNAL */
171
172 extern "C" int
mbdyn_stop_at_end_of_iteration(void)173 mbdyn_stop_at_end_of_iteration(void)
174 {
175 #ifdef HAVE_SIGNAL
176 return ::mbdyn_keep_going >= MBDYN_STOP_AT_END_OF_ITERATION;
177 #else // ! HAVE_SIGNAL
178 return 0;
179 #endif // ! HAVE_SIGNAL
180 }
181
182 extern "C" int
mbdyn_stop_at_end_of_time_step(void)183 mbdyn_stop_at_end_of_time_step(void)
184 {
185 #ifdef HAVE_SIGNAL
186 return ::mbdyn_keep_going >= MBDYN_STOP_AT_END_OF_TIME_STEP;
187 #else // ! HAVE_SIGNAL
188 return 0;
189 #endif // ! HAVE_SIGNAL
190 }
191
192 extern "C" void
mbdyn_set_stop_at_end_of_iteration(void)193 mbdyn_set_stop_at_end_of_iteration(void)
194 {
195 #ifdef HAVE_SIGNAL
196 ::mbdyn_keep_going = MBDYN_STOP_AT_END_OF_ITERATION;
197 #endif // HAVE_SIGNAL
198 }
199
200 extern "C" void
mbdyn_set_stop_at_end_of_time_step(void)201 mbdyn_set_stop_at_end_of_time_step(void)
202 {
203 #ifdef HAVE_SIGNAL
204 ::mbdyn_keep_going = MBDYN_STOP_AT_END_OF_TIME_STEP;
205 #endif // HAVE_SIGNAL
206 }
207
208 extern "C" void
mbdyn_signal_init(int pre)209 mbdyn_signal_init(int pre)
210 {
211 #ifdef HAVE_SIGNAL
212 #if defined(__FreeBSD__) || defined(__DragonFly__)
213 __sighandler_t *hdl;
214 #else // ! defined(__FreeBSD__)
215 __sighandler_t hdl;
216 #endif // ! defined(__FreeBSD__)
217
218 if (pre) {
219 hdl = mbdyn_really_exit_handler;
220
221 } else {
222 hdl = mbdyn_modify_final_time_handler;
223 }
224 /*
225 * FIXME: don't do this if compiling with USE_RTAI
226 * Re FIXME: use sigaction() ...
227 */
228 ::mbdyn_sh_term = signal(SIGTERM, hdl);
229 if (::mbdyn_sh_term == SIG_IGN) {
230 signal(SIGTERM, SIG_IGN);
231 }
232
233 ::mbdyn_sh_int = signal(SIGINT, hdl);
234 if (::mbdyn_sh_int == SIG_IGN) {
235 signal(SIGINT, SIG_IGN);
236 }
237
238 #ifdef SIGHUP
239 ::mbdyn_sh_hup = signal(SIGHUP, hdl);
240 if (::mbdyn_sh_hup == SIG_IGN) {
241 signal(SIGHUP, SIG_IGN);
242 }
243 #endif // SIGHUP
244
245 #ifdef SIGPIPE
246 ::mbdyn_sh_pipe = signal(SIGPIPE, hdl);
247 if (::mbdyn_sh_pipe == SIG_IGN) {
248 signal(SIGPIPE, SIG_IGN);
249 }
250 #endif // SIGPIPE
251 #endif /* HAVE_SIGNAL */
252 }
253
254 int
mbdyn_reserve_stack(unsigned long size)255 mbdyn_reserve_stack(unsigned long size)
256 {
257 int buf[size];
258
259 #ifdef HAVE_MEMSET
260 memset(buf, 0, size*sizeof(int));
261 #else /* !HAVE_MEMSET */
262 for (unsigned long i = 0; i < size; i++) {
263 buf[i] = 0;
264 }
265 #endif /* !HAVE_MEMSET */
266
267 #ifdef HAVE_MLOCKALL
268 return mlockall(MCL_CURRENT | MCL_FUTURE);
269 #else /* !HAVE_MLOCKALL */
270 return 0;
271 #endif /* !HAVE_MLOCKALL */
272 }
273
274 /* Costruttore: esegue la simulazione */
Solver(MBDynParser & HPar,const std::string & sInFName,const std::string & sOutFName,unsigned int nThreads,bool bPar)275 Solver::Solver(MBDynParser& HPar,
276 const std::string& sInFName,
277 const std::string& sOutFName,
278 unsigned int nThreads,
279 bool bPar)
280 :
281 #ifdef USE_MULTITHREAD
282 nThreads(nThreads),
283 #endif /* USE_MULTITHREAD */
284 pTSC(0),
285 dCurrTimeStep(0.),
286 iStIter(0),
287 dTime(0.),
288 MaxTimeStep(new ConstDriveCaller(::dDefaultMaxTimeStep)),
289 dMinTimeStep(::dDefaultMinTimeStep),
290 CurrStep(StepIntegrator::NEWSTEP),
291 sInputFileName(sInFName),
292 sOutputFileName(sOutFName),
293 HP(HPar),
294 iMaxIterations(::iDefaultMaxIterations),
295 EigAn(),
296 pRTSolver(0),
297 iNumPreviousVectors(2),
298 iUnkStates(1),
299 pdWorkSpace(0),
300 qX(),
301 qXPrime(),
302 pX(0),
303 pXPrime(0),
304 dInitialTime(0.),
305 dFinalTime(0.),
306 dRefTimeStep(0.),
307 dInitialTimeStep(1.),
308 dMaxResidual(std::numeric_limits<doublereal>::max()),
309 dMaxResidualDiff(std::numeric_limits<doublereal>::max()),
310 eTimeStepLimit(TS_SOFT_LIMIT),
311 iDummyStepsNumber(::iDefaultDummyStepsNumber),
312 dDummyStepsRatio(::dDefaultDummyStepsRatio),
313 eAbortAfter(AFTER_UNKNOWN),
314 RegularType(INT_UNKNOWN),
315 DummyType(INT_UNKNOWN),
316 pDerivativeSteps(0),
317 pFirstDummyStep(0),
318 pDummySteps(0),
319 pFirstRegularStep(0),
320 pRegularSteps(0),
321 pCurrStepIntegrator(0),
322 pRhoRegular(0),
323 pRhoAlgebraicRegular(0),
324 pRhoDummy(0),
325 pRhoAlgebraicDummy(0),
326 dDerivativesCoef(::dDefaultDerivativesCoefficient),
327 CurrLinearSolver(),
328 ResTest(NonlinearSolverTest::NORM),
329 SolTest(NonlinearSolverTest::NONE),
330 bScale(false),
331 bTrueNewtonRaphson(true),
332 bKeepJac(false),
333 iIterationsBeforeAssembly(0),
334 NonlinearSolverType(NonlinearSolver::UNKNOWN),
335 /* for matrix-free solvers */
336 MFSolverType(MatrixFreeSolver::UNKNOWN),
337 dIterTol(::dDefaultTol),
338 PcType(Preconditioner::FULLJACOBIANMATRIX),
339 iPrecondSteps(::iDefaultPreconditionerSteps),
340 iIterativeMaxSteps(::iDefaultPreconditionerSteps),
341 dIterertiveEtaMax(defaultIterativeEtaMax),
342 dIterertiveTau(defaultIterativeTau),
343 /* end of matrix-free solvers */
344 /* for line search solver */
345 LineSearch(),
346 /* end of line search solver */
347 /* for parallel solvers */
348 bParallel(bPar),
349 pSDM(0),
350 iNumLocDofs(0),
351 iNumIntDofs(0),
352 pLocDofs(0),
353 pIntDofs(0),
354 pDofs(0),
355 pLocalSM(0),
356 /* end of parallel solvers */
357 pDM(0),
358 iNumDofs(0),
359 pSM(0),
360 pNLS(0),
361 eStatus(SOLVER_STATUS_UNINITIALIZED),
362 bOutputCounter(false),
363 outputCounterPrefix(),
364 outputCounterPostfix(),
365 iTotIter(0),
366 dTotErr(0.),
367 dTest(std::numeric_limits<double>::max()),
368 dSolTest(std::numeric_limits<double>::max()),
369 bSolConv(false),
370 bOut(false),
371 lStep(0)
372 {
373 DEBUGCOUTFNAME("Solver::Solver");
374 ::InitTimeStepData();
375 ASSERT(!sInFName.empty());
376 }
377
378 #ifdef USE_MULTITHREAD
379 void
ThreadPrepare(void)380 Solver::ThreadPrepare(void)
381 {
382 /* check for thread potential */
383 if (nThreads == 0) {
384 int n = get_nprocs();
385
386 if (n > 1) {
387 silent_cout("no multithread requested "
388 "with a potential of " << n
389 << " CPUs" << std::endl);
390 nThreads = n;
391
392 } else {
393 nThreads = 1;
394 }
395 }
396 }
397 #endif /* USE_MULTITHREAD */
398
399 bool
Prepare(void)400 Solver::Prepare(void)
401 {
402 DEBUGCOUTFNAME("Solver::Prepare");
403
404 // consistency check
405 if (eStatus != SOLVER_STATUS_UNINITIALIZED) {
406 silent_cerr("Prepare() must be called first" << std::endl);
407 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
408 }
409
410 mbdyn_signal_init(1);
411
412 /* Legge i dati relativi al metodo di integrazione */
413 ReadData(HP);
414
415 #ifdef USE_MULTITHREAD
416 ThreadPrepare();
417 #endif /* USE_MULTITHREAD */
418
419 if (pRTSolver) {
420 pRTSolver->Setup();
421 }
422
423 #ifdef USE_MPI
424 int mpi_finalize = 0;
425 #endif
426
427 #ifdef USE_SCHUR
428 if (bParallel) {
429 DEBUGLCOUT(MYDEBUG_MEM, "creating parallel SchurDataManager"
430 << std::endl);
431
432 SAFENEWWITHCONSTRUCTOR(pSDM,
433 SchurDataManager,
434 SchurDataManager(HP,
435 OutputFlags,
436 this,
437 dInitialTime,
438 sOutputFileName.c_str(),
439 sInputFileName.c_str(),
440 eAbortAfter == AFTER_INPUT));
441
442 pDM = pSDM;
443
444 } else
445 #endif // USE_SCHUR
446 {
447 /* chiama il gestore dei dati generali della simulazione */
448 #ifdef USE_MULTITHREAD
449 if (nThreads > 1) {
450 if (!(CurrLinearSolver.GetSolverFlags() & LinSol::SOLVER_FLAGS_ALLOWS_MT_ASS)) {
451 /* conservative: dir may use too much memory */
452 if (!CurrLinearSolver.AddSolverFlags(LinSol::SOLVER_FLAGS_ALLOWS_MT_ASS)) {
453 bool b;
454
455 #if defined(USE_UMFPACK)
456 b = CurrLinearSolver.SetSolver(LinSol::UMFPACK_SOLVER,
457 LinSol::SOLVER_FLAGS_ALLOWS_MT_ASS);
458 #elif defined(USE_Y12)
459 b = CurrLinearSolver.SetSolver(LinSol::Y12_SOLVER,
460 LinSol::SOLVER_FLAGS_ALLOWS_MT_ASS);
461 #else
462 b = false;
463 #endif
464 if (!b) {
465 silent_cerr("unable to select a CC-capable solver"
466 << std::endl);
467 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
468 }
469 }
470 }
471
472 silent_cout("Creating multithread solver "
473 "with " << nThreads << " threads "
474 "and "
475 << CurrLinearSolver.GetSolverName()
476 << " linear solver"
477 << std::endl);
478
479 SAFENEWWITHCONSTRUCTOR(pDM,
480 MultiThreadDataManager,
481 MultiThreadDataManager(HP,
482 OutputFlags,
483 this,
484 dInitialTime,
485 sOutputFileName.c_str(),
486 sInputFileName.c_str(),
487 eAbortAfter == AFTER_INPUT,
488 nThreads));
489
490 } else
491 #endif /* USE_MULTITHREAD */
492 {
493 DEBUGLCOUT(MYDEBUG_MEM, "creating DataManager"
494 << std::endl);
495
496 silent_cout("Creating scalar solver "
497 "with "
498 << CurrLinearSolver.GetSolverName()
499 << " linear solver"
500 << std::endl);
501
502 SAFENEWWITHCONSTRUCTOR(pDM,
503 DataManager,
504 DataManager(HP,
505 OutputFlags,
506 this,
507 dInitialTime,
508 sOutputFileName.c_str(),
509 sInputFileName.c_str(),
510 eAbortAfter == AFTER_INPUT));
511 }
512 }
513
514 // log symbol table
515 std::ostream& log = pDM->GetLogFile();
516 log << "Symbol table:" << std::endl;
517 log << HP.GetMathParser().GetSymbolTable();
518
519 #ifdef HAVE_ENVIRON
520 // log environment
521 log << "Environment:" << std::endl;
522 for (int i = 0; environ[i] != NULL; i++) {
523 log << " " << environ[i] << std::endl;
524 }
525 #endif // HAVE_ENVIRON
526
527 // close input stream
528 HP.Close();
529
530 /* Si fa dare il DriveHandler e linka i drivers di rho ecc. */
531 const DriveHandler* pDH = pDM->pGetDrvHdl();
532 SetOutputDriveHandler(pDH);
533
534 bOutputCounter = outputCounter() && isatty(fileno(stderr));
535 outputCounterPrefix = bOutputCounter ? "\n" : "";
536 outputCounterPostfix = outputStep() ? "\n" : "\r";
537
538 /* Si fa dare l'std::ostream al file di output per il log */
539 std::ostream& Out = pDM->GetOutFile();
540
541 if (eAbortAfter == AFTER_INPUT) {
542 /* Esce */
543 pDM->Output(0, dTime, 0., true);
544 Out << "End of Input; no simulation or assembly is required."
545 << std::endl;
546 return false;
547
548 } else if (eAbortAfter == AFTER_ASSEMBLY) {
549 /* Fa l'output dell'assemblaggio iniziale e poi esce */
550 pDM->Output(0, dTime, 0., true);
551 Out << "End of Initial Assembly; no simulation is required."
552 << std::endl;
553 return false;
554 }
555
556 #ifdef USE_SCHUR
557 /* Qui crea le partizioni: principale fra i processi, se parallelo */
558 if (bParallel) {
559 pSDM->CreatePartition();
560 }
561 #endif // USE_SCHUR
562
563 pRegularSteps->SetDriveHandler(pDH);
564 if (iDummyStepsNumber) {
565 pDummySteps->SetDriveHandler(pDH);
566 }
567
568 /* Costruisce i vettori della soluzione ai vari passi */
569 DEBUGLCOUT(MYDEBUG_MEM, "creating solution vectors" << std::endl);
570
571 #ifdef USE_SCHUR
572 if (bParallel) {
573 iNumDofs = pSDM->HowManyDofs(SchurDataManager::TOTAL);
574 pDofs = pSDM->pGetDofsList();
575
576 iNumLocDofs = pSDM->HowManyDofs(SchurDataManager::LOCAL);
577 pLocDofs = pSDM->GetDofsList(SchurDataManager::LOCAL);
578
579 iNumIntDofs = pSDM->HowManyDofs(SchurDataManager::INTERNAL);
580 pIntDofs = pSDM->GetDofsList(SchurDataManager::INTERNAL);
581
582 } else
583 #endif // USE_SCHUR
584 {
585 iNumDofs = pDM->iGetNumDofs();
586 }
587
588 /* relink those known drive callers that might need
589 * the data manager, but were verated ahead of it */
590 pTSC->SetDriveHandler(pDM->pGetDrvHdl());
591
592 ASSERT(iNumDofs > 0);
593
594 integer iRSteps = pRegularSteps->GetIntegratorNumPreviousStates();
595 integer iFSteps = 0;
596 if (iDummyStepsNumber) {
597 iFSteps = pDummySteps->GetIntegratorNumPreviousStates();
598 }
599 iNumPreviousVectors = (iRSteps < iFSteps) ? iFSteps : iRSteps;
600
601 integer iRUnkStates = pRegularSteps->GetIntegratorNumUnknownStates();
602 integer iFUnkStates = 0;
603 if (iDummyStepsNumber) {
604 iFUnkStates = pDummySteps->GetIntegratorNumUnknownStates();
605 }
606 iUnkStates = (iRUnkStates < iFUnkStates) ? iFUnkStates : iRUnkStates;
607
608 /* allocate workspace for previous time steps */
609 SAFENEWARR(pdWorkSpace, doublereal,
610 2*(iNumPreviousVectors)*iNumDofs);
611 /* allocate MyVectorHandlers for previous time steps: use workspace */
612 for (int ivec = 0; ivec < iNumPreviousVectors; ivec++) {
613 SAFENEWWITHCONSTRUCTOR(pX,
614 MyVectorHandler,
615 MyVectorHandler(iNumDofs, pdWorkSpace+ivec*iNumDofs));
616 qX.push_back(pX);
617 SAFENEWWITHCONSTRUCTOR(pXPrime,
618 MyVectorHandler,
619 MyVectorHandler(iNumDofs,
620 pdWorkSpace+((iNumPreviousVectors)+ivec)*iNumDofs));
621 qXPrime.push_back(pXPrime);
622 pX = 0;
623 pXPrime = 0;
624 }
625 /* allocate MyVectorHandlers for unknown time step(s): own memory */
626 SAFENEWWITHCONSTRUCTOR(pX,
627 MyVectorHandler,
628 MyVectorHandler(iUnkStates*iNumDofs));
629 SAFENEWWITHCONSTRUCTOR(pXPrime,
630 MyVectorHandler,
631 MyVectorHandler(iUnkStates*iNumDofs));
632
633
634 /* Resetta i vettori */
635 for (int ivec = 0; ivec < iNumPreviousVectors; ivec++) {
636 qX[ivec]->Reset();
637 qXPrime[ivec]->Reset();
638 }
639 pX->Reset();
640 pXPrime->Reset();
641
642 /*
643 * Immediately link DataManager to current solution
644 *
645 * this should work as long as the last unknown time step is put
646 * at the beginning of pX, pXPrime
647 */
648 pDM->LinkToSolution(*pX, *pXPrime);
649
650 /* a questo punto si costruisce il nonlinear solver */
651 pNLS = AllocateNonlinearSolver();
652
653 Scale.ResizeReset(iNumDofs);
654 if (bScale) {
655 /* collects scale factors from data manager */
656 pDM->SetScale(Scale);
657 }
658
659 /*
660 * prepare tests for nonlinear solver;
661 *
662 * test on residual may allow pre-scaling;
663 * test on solution (difference between two iterations) does not
664 */
665 NonlinearSolverTest *pResTest = 0;
666 if (bScale) {
667 NonlinearSolverTestScale *pResTestScale = 0;
668
669 switch (ResTest) {
670 case NonlinearSolverTest::NORM:
671 SAFENEW(pResTestScale, NonlinearSolverTestScaleNorm);
672 break;
673
674 case NonlinearSolverTest::MINMAX:
675 SAFENEW(pResTestScale, NonlinearSolverTestScaleMinMax);
676 break;
677
678 default:
679 ASSERT(0);
680 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
681 }
682
683 /* registers scale factors at nonlinear solver */
684 pResTestScale->SetScale(&Scale);
685
686 pResTest = pResTestScale;
687
688
689 } else {
690 switch (ResTest) {
691 case NonlinearSolverTest::NONE:
692 SAFENEW(pResTest, NonlinearSolverTestNone);
693 break;
694
695 case NonlinearSolverTest::NORM:
696 SAFENEW(pResTest, NonlinearSolverTestNorm);
697 break;
698
699 case NonlinearSolverTest::MINMAX:
700 SAFENEW(pResTest, NonlinearSolverTestMinMax);
701 break;
702
703 default:
704 ASSERT(0);
705 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
706 }
707 }
708
709 NonlinearSolverTest *pSolTest = 0;
710 switch (SolTest) {
711 case NonlinearSolverTest::NONE:
712 SAFENEW(pSolTest, NonlinearSolverTestNone);
713 break;
714
715 case NonlinearSolverTest::NORM:
716 SAFENEW(pSolTest, NonlinearSolverTestNorm);
717 break;
718
719 case NonlinearSolverTest::MINMAX:
720 SAFENEW(pSolTest, NonlinearSolverTestMinMax);
721 break;
722
723 default:
724 ASSERT(0);
725 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
726 }
727
728 /* registers tests in nonlinear solver */
729 pNLS->SetTest(pResTest, pSolTest);
730
731 /*
732 * Dell'assemblaggio iniziale dei vincoli se ne occupa il DataManager
733 * in quanto e' lui il responsabile dei dati della simulazione,
734 * e quindi anche della loro coerenza. Inoltre e' lui a sapere
735 * quali equazioni sono di vincolo o meno.
736 */
737
738 pDM->SetValue(*pX, *pXPrime);
739
740
741 /*
742 * Prepare output
743 */
744 pDM->OutputPrepare();
745
746 /*
747 * If eigenanalysis is requested, prepare output for it
748 */
749 if (EigAn.bAnalysis) {
750 pDM->OutputEigPrepare(EigAn.Analyses.size(), iNumDofs);
751 if (EigAn.iFNameWidth == EigenAnalysis::EIGAN_WIDTH_COMPUTE) {
752 EigAn.iFNameWidth = int(std::log10(double(EigAn.Analyses.size()))) + 1;
753 }
754 }
755
756 /*
757 * Dialoga con il DataManager per dargli il tempo iniziale
758 * e per farsi inizializzare i vettori di soluzione e derivata
759 */
760 /* FIXME: the time is already set by DataManager, but FileDrivers
761 * have not been ServePending'd
762 */
763 dTime = dInitialTime;
764 pDM->SetTime(dTime, dInitialTimeStep, 0);
765
766 EigAn.currAnalysis = std::find_if(EigAn.Analyses.begin(), EigAn.Analyses.end(),
767 bind2nd(std::greater<doublereal>(), dTime));
768 if (EigAn.currAnalysis != EigAn.Analyses.end() && EigAn.currAnalysis != EigAn.Analyses.begin()) {
769 --EigAn.currAnalysis;
770 }
771
772 // if eigenanalysis is requested and currAnalysis points
773 // past the end of the array, the analysis was requested
774 // at Time < initial time; perform *before* derivatives
775 if (EigAn.bAnalysis
776 && ((EigAn.currAnalysis == EigAn.Analyses.end()
777 && EigAn.Analyses.back() < dTime)
778 || (EigAn.currAnalysis != EigAn.Analyses.end()
779 && *EigAn.currAnalysis < dTime)))
780 {
781 Eig();
782 if (EigAn.currAnalysis != EigAn.Analyses.end()) {
783 ++EigAn.currAnalysis;
784 }
785 }
786
787 /* calcolo delle derivate */
788 DEBUGLCOUT(MYDEBUG_DERIVATIVES, "derivatives solution step"
789 << std::endl);
790
791 mbdyn_signal_init(0);
792
793 /* settaggio degli output Types */
794 unsigned OF = OutputFlags;
795 if ( DEBUG_LEVEL_MATCH(MYDEBUG_RESIDUAL) ) {
796 OF |= OUTPUT_RES;
797 }
798 if ( DEBUG_LEVEL_MATCH(MYDEBUG_JAC) ) {
799 OF |= OUTPUT_JAC;
800 }
801 if ( DEBUG_LEVEL_MATCH(MYDEBUG_SOL) ) {
802 OF |= OUTPUT_SOL;
803 }
804 pNLS->SetOutputFlags(OF);
805 if (pOutputMeter) {
806 pOutputMeter->SetDrvHdl(pDM->pGetDrvHdl());
807 pNLS->SetOutputMeter(pOutputMeter->pCopy());
808 }
809
810 pDerivativeSteps->SetDataManager(pDM);
811 pDerivativeSteps->OutputTypes(DEBUG_LEVEL_MATCH(MYDEBUG_PRED));
812 if (iDummyStepsNumber) {
813 pFirstDummyStep->SetDataManager(pDM);
814 pFirstDummyStep->OutputTypes(DEBUG_LEVEL_MATCH(MYDEBUG_PRED));
815 pDummySteps->SetDataManager(pDM);
816 pDummySteps->OutputTypes(DEBUG_LEVEL_MATCH(MYDEBUG_PRED));
817 }
818 pFirstRegularStep->SetDataManager(pDM);
819 pFirstRegularStep->OutputTypes(DEBUG_LEVEL_MATCH(MYDEBUG_PRED));
820 pRegularSteps->SetDataManager(pDM);
821 pRegularSteps->OutputTypes(DEBUG_LEVEL_MATCH(MYDEBUG_PRED));
822
823 #ifdef USE_EXTERNAL
824 pNLS->SetExternal(External::EMPTY);
825 #endif /* USE_EXTERNAL */
826 /* Setup SolutionManager(s) */
827 SetupSolmans(pDerivativeSteps->GetIntegratorNumUnknownStates());
828
829 /* Derivative steps */
830 pCurrStepIntegrator = pDerivativeSteps;
831 try {
832 if (outputStep()) {
833 if (outputCounter()) {
834 silent_cout(std::endl);
835 }
836 silent_cout("Derivatives t=" << dTime << " coef=" << dDerivativesCoef << std::endl);
837 }
838 dTest = pDerivativeSteps->Advance(this,
839 dInitialTimeStep, 1., StepIntegrator::NEWSTEP,
840 qX, qXPrime, pX, pXPrime,
841 iStIter, dTest, dSolTest);
842 }
843 catch (NonlinearSolver::NoConvergence) {
844 silent_cerr("Initial derivatives calculation " << iStIter
845 << " does not converge; aborting..." << std::endl
846 << "(hint: try playing with the \"derivatives coefficient\" value)" << std::endl);
847 pDM->Output(0, dTime, 0., true);
848 throw ErrMaxIterations(MBDYN_EXCEPT_ARGS);
849 }
850 catch (NonlinearSolver::ErrSimulationDiverged) {
851 /*
852 * Mettere qui eventuali azioni speciali
853 * da intraprendere in caso di errore ...
854 */
855 throw SimulationDiverged(MBDYN_EXCEPT_ARGS);
856 }
857 catch (LinearSolver::ErrFactor& err) {
858 /*
859 * Mettere qui eventuali azioni speciali
860 * da intraprendere in caso di errore ...
861 */
862 silent_cerr("Initial derivatives failed because no pivot element "
863 "could be found for column " << err.iCol
864 << " (" << pDM->GetDofDescription(err.iCol) << "); "
865 "aborting..." << std::endl);
866 throw SimulationDiverged(MBDYN_EXCEPT_ARGS);
867 }
868 catch (NonlinearSolver::ConvergenceOnSolution) {
869 bSolConv = true;
870 }
871 catch (EndOfSimulation& eos) {
872 silent_cerr("Simulation ended during the derivatives steps:\n" << eos.what() << "\n");
873 return false;
874 }
875
876 SAFEDELETE(pDerivativeSteps);
877 pDerivativeSteps = 0;
878
879 #if 0
880 /* don't sum up the derivatives error */
881 dTotErr += dTest;
882 #endif
883 iTotIter += iStIter;
884
885 if (outputMsg()) {
886 Out << "# Derivatives solution step at time " << dInitialTime
887 << " performed in " << iStIter
888 << " iterations with " << dTest
889 << " error" << std::endl;
890 }
891
892 DEBUGCOUT("Derivatives solution step has been performed successfully"
893 " in " << iStIter << " iterations" << std::endl);
894
895 #ifdef USE_EXTERNAL
896 /* comunica che gli ultimi dati inviati sono la condizione iniziale */
897 if (iDummyStepsNumber == 0) {
898 External::SendInitial();
899 }
900 #endif /* USE_EXTERNAL */
901
902 if (eAbortAfter == AFTER_DERIVATIVES) {
903 /*
904 * Fa l'output della soluzione delle derivate iniziali ed esce
905 */
906
907
908 pDM->Output(0, dTime, 0., true);
909 Out << "End of derivatives; no simulation is required."
910 << std::endl;
911 return false;
912
913 } else if (mbdyn_stop_at_end_of_time_step()) {
914 /*
915 * Fa l'output della soluzione delle derivate iniziali ed esce
916 */
917 pDM->Output(0, dTime, 0., true);
918 Out << "Interrupted during derivatives computation." << std::endl;
919 throw ErrInterrupted(MBDYN_EXCEPT_ARGS);
920 }
921
922 // if eigenanalysis is requested and currAnalysis points
923 // past the end of the array, the analysis was requested
924 // at Time == initial time; perform *after* derivatives
925 if (EigAn.bAnalysis) {
926 ASSERT(EigAn.Analyses.size() > 0);
927
928 if ((EigAn.currAnalysis == EigAn.Analyses.end()
929 && EigAn.Analyses.back() == dTime)
930 || (EigAn.currAnalysis != EigAn.Analyses.end()
931 && *EigAn.currAnalysis == dTime))
932 {
933 Eig();
934 if (EigAn.currAnalysis != EigAn.Analyses.end()) {
935 ++EigAn.currAnalysis;
936 }
937 }
938 }
939
940 /* Dati comuni a passi fittizi e normali */
941 lStep = 1;
942
943 if (iDummyStepsNumber > 0) {
944 /* passi fittizi */
945
946 /*
947 * inizio integrazione: primo passo a predizione lineare
948 * con sottopassi di correzione delle accelerazioni
949 * e delle reazioni vincolari
950 */
951 pDM->BeforePredict(*pX, *pXPrime, *qX[0], *qXPrime[0]);
952 Flip();
953
954 dRefTimeStep = dInitialTimeStep*dDummyStepsRatio;
955 dCurrTimeStep = dRefTimeStep;
956 /* FIXME: do we need to serve pending drives in dummy steps? */
957 pDM->SetTime(dTime + dCurrTimeStep, dCurrTimeStep, 0);
958
959 DEBUGLCOUT(MYDEBUG_FSTEPS, "Current time step: "
960 << dCurrTimeStep << std::endl);
961
962 if (outputStep()) {
963 silent_cout("Dummy Step(" << lStep << ") t=" << dTime + dCurrTimeStep << " dt=" << dCurrTimeStep << std::endl);
964 }
965
966 ASSERT(pFirstDummyStep != 0);
967
968 /* Setup SolutionManager(s) */
969 SetupSolmans(pFirstDummyStep->GetIntegratorNumUnknownStates());
970 /* pFirstDummyStep */
971 pCurrStepIntegrator = pFirstDummyStep;
972 try {
973 dTest = pFirstDummyStep->Advance(this,
974 dRefTimeStep, 1.,
975 StepIntegrator::NEWSTEP,
976 qX, qXPrime, pX, pXPrime,
977 iStIter, dTest, dSolTest);
978 }
979 catch (NonlinearSolver::NoConvergence) {
980 silent_cerr("First dummy step does not converge; "
981 "TimeStep=" << dCurrTimeStep
982 << " cannot be reduced further; "
983 "aborting..." << std::endl);
984 pDM->Output(0, dTime, dCurrTimeStep, true);
985 throw ErrMaxIterations(MBDYN_EXCEPT_ARGS);
986 }
987 catch (NonlinearSolver::ErrSimulationDiverged) {
988 /*
989 * Mettere qui eventuali azioni speciali
990 * da intraprendere in caso di errore ...
991 */
992 throw SimulationDiverged(MBDYN_EXCEPT_ARGS);
993 }
994 catch (LinearSolver::ErrFactor& err) {
995 /*
996 * Mettere qui eventuali azioni speciali
997 * da intraprendere in caso di errore ...
998 */
999 silent_cerr("First dummy step failed because no pivot element "
1000 "could be found for column " << err.iCol
1001 << " (" << pDM->GetDofDescription(err.iCol) << "); "
1002 "aborting..." << std::endl);
1003 throw SimulationDiverged(MBDYN_EXCEPT_ARGS);
1004 }
1005 catch (NonlinearSolver::ConvergenceOnSolution) {
1006 bSolConv = true;
1007 }
1008 catch (EndOfSimulation& eos) {
1009 silent_cerr("Simulation ended during the first dummy step:\n"
1010 << eos.what() << "\n");
1011 return false;
1012 }
1013
1014 SAFEDELETE(pFirstDummyStep);
1015 pFirstDummyStep = 0;
1016
1017 dRefTimeStep = dCurrTimeStep;
1018 dTime += dRefTimeStep;
1019
1020 #if 0
1021 /* don't sum up the derivatives error */
1022 dTotErr += dTest;
1023 #endif
1024 iTotIter += iStIter;
1025
1026 if (mbdyn_stop_at_end_of_time_step()) {
1027 /*
1028 * Fa l'output della soluzione delle derivate iniziali
1029 * ed esce
1030 */
1031 #ifdef DEBUG_FICTITIOUS
1032 pDM->Output(0, dTime, dCurrTimeStep, true);
1033 #endif /* DEBUG_FICTITIOUS */
1034 Out << "Interrupted during first dummy step." << std::endl;
1035 throw ErrInterrupted(MBDYN_EXCEPT_ARGS);
1036 }
1037
1038 #ifdef DEBUG_FICTITIOUS
1039 pDM->Output(0, dTime, dCurrTimeStep, true);
1040 #endif /* DEBUG_FICTITIOUS */
1041
1042 /* Passi fittizi successivi */
1043 if (iDummyStepsNumber > 1) {
1044 /* Setup SolutionManager(s) */
1045 SetupSolmans(pDummySteps->GetIntegratorNumUnknownStates());
1046 }
1047
1048 for (int iSubStep = 2;
1049 iSubStep <= iDummyStepsNumber;
1050 iSubStep++)
1051 {
1052 pDM->BeforePredict(*pX, *pXPrime,
1053 *qX[0], *qXPrime[0]);
1054 Flip();
1055
1056 DEBUGLCOUT(MYDEBUG_FSTEPS, "Dummy step "
1057 << iSubStep
1058 << "; current time step: " << dCurrTimeStep
1059 << std::endl);
1060
1061 if (outputStep()) {
1062 silent_cout("Dummy Step(" << iSubStep << ") t=" << dTime + dCurrTimeStep << " dt=" << dCurrTimeStep << std::endl);
1063 }
1064
1065 pCurrStepIntegrator = pDummySteps;
1066 ASSERT(pDummySteps!= 0);
1067 try {
1068 pDM->SetTime(dTime + dCurrTimeStep, dCurrTimeStep, 0);
1069 dTest = pDummySteps->Advance(this,
1070 dRefTimeStep,
1071 dCurrTimeStep/dRefTimeStep,
1072 StepIntegrator::NEWSTEP,
1073 qX, qXPrime, pX, pXPrime,
1074 iStIter, dTest, dSolTest);
1075 }
1076 catch (NonlinearSolver::NoConvergence) {
1077 silent_cerr("Dummy step " << iSubStep
1078 << " does not converge; "
1079 "TimeStep=" << dCurrTimeStep
1080 << " cannot be reduced further; "
1081 "aborting..." << std::endl);
1082 pDM->Output(0, dTime, dCurrTimeStep, true);
1083 throw ErrMaxIterations(MBDYN_EXCEPT_ARGS);
1084 }
1085
1086 catch (NonlinearSolver::ErrSimulationDiverged) {
1087 /*
1088 * Mettere qui eventuali azioni speciali
1089 * da intraprendere in caso di errore ...
1090 */
1091 throw SimulationDiverged(MBDYN_EXCEPT_ARGS);
1092 }
1093 catch (LinearSolver::ErrFactor& err) {
1094 /*
1095 * Mettere qui eventuali azioni speciali
1096 * da intraprendere in caso di errore ...
1097 */
1098 silent_cerr("Dummy step " << iSubStep
1099 << " failed because no pivot element "
1100 "could be found for column " << err.iCol
1101 << " (" << pDM->GetDofDescription(err.iCol) << "); "
1102 "aborting..." << std::endl);
1103 throw SimulationDiverged(MBDYN_EXCEPT_ARGS);
1104 }
1105 catch (NonlinearSolver::ConvergenceOnSolution) {
1106 bSolConv = true;
1107 }
1108 catch (EndOfSimulation& eos) {
1109 silent_cerr("Simulation ended during the dummy steps:\n"
1110 << eos.what() << "\n");
1111 return false;
1112 }
1113
1114 #if 0
1115 /* don't sum up the derivatives error */
1116 dTotErr += dTest;
1117 #endif
1118 iTotIter += iStIter;
1119
1120 #ifdef DEBUG
1121 if (DEBUG_LEVEL(MYDEBUG_FSTEPS)) {
1122 Out << "Step " << lStep
1123 << " time " << dTime + dCurrTimeStep
1124 << " step " << dCurrTimeStep
1125 << " iterations " << iStIter
1126 << " error " << dTest << std::endl;
1127 }
1128 #endif /* DEBUG */
1129
1130 DEBUGLCOUT(MYDEBUG_FSTEPS, "Substep " << iSubStep
1131 << " of step " << lStep
1132 << " has been successfully completed "
1133 "in " << iStIter << " iterations"
1134 << std::endl);
1135
1136 if (mbdyn_stop_at_end_of_time_step()) {
1137 /* */
1138 #ifdef DEBUG_FICTITIOUS
1139 pDM->Output(0, dTime, dCurrTimeStep);
1140 #endif /* DEBUG_FICTITIOUS */
1141 Out << "Interrupted during dummy steps."
1142 << std::endl;
1143 throw ErrInterrupted(MBDYN_EXCEPT_ARGS);
1144 }
1145
1146 dTime += dRefTimeStep;
1147 }
1148 if (outputMsg()) {
1149 Out << "# Initial solution after dummy steps "
1150 "at time " << dTime
1151 << " performed in " << iStIter
1152 << " iterations with " << dTest
1153 << " error" << std::endl;
1154 }
1155
1156 DEBUGLCOUT(MYDEBUG_FSTEPS,
1157 "Dummy steps have been successfully completed "
1158 "in " << iStIter << " iterations" << std::endl);
1159 #ifdef USE_EXTERNAL
1160 /* comunica che gli ultimi dati inviati sono la condizione iniziale */
1161 External::SendInitial();
1162 #endif /* USE_EXTERNAL */
1163 } /* Fine dei passi fittizi */
1164
1165 /* Output delle "condizioni iniziali" */
1166 bOut = pDM->Output(0, dTime, dCurrTimeStep);
1167
1168 if (outputMsg()) {
1169 Out
1170 << "# Key for lines starting with \"Step\":"
1171 << std::endl
1172 << "# Step Time TStep NIter ResErr SolErr SolConv Out"
1173 << std::endl
1174 << "Step " << 0
1175 << " " << dTime + dCurrTimeStep
1176 << " " << dCurrTimeStep
1177 << " " << iStIter
1178 << " " << dTest
1179 << " " << dSolTest
1180 << " " << bSolConv
1181 << " " << bOut
1182 << std::endl;
1183 }
1184
1185
1186 if (eAbortAfter == AFTER_DUMMY_STEPS) {
1187 Out << "End of dummy steps; no simulation is required."
1188 << std::endl;
1189 return false;
1190
1191 } else if (mbdyn_stop_at_end_of_time_step()) {
1192 /* Fa l'output della soluzione ed esce */
1193 Out << "Interrupted during dummy steps." << std::endl;
1194 throw ErrInterrupted(MBDYN_EXCEPT_ARGS);
1195 }
1196
1197 eStatus = SOLVER_STATUS_PREPARED;
1198
1199 return true;
1200 }
1201
1202 bool
Start(void)1203 Solver::Start(void)
1204 {
1205 DEBUGCOUTFNAME("Solver::Start");
1206
1207 // consistency check
1208 if (eStatus != SOLVER_STATUS_PREPARED) {
1209 silent_cerr("Start() must be called after Prepare()" << std::endl);
1210 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
1211 }
1212
1213 /* primo passo regolare */
1214
1215 #ifdef USE_EXTERNAL
1216 /* il prossimo passo e' un regular */
1217 pNLS->SetExternal(External::REGULAR);
1218 #endif /* USE_EXTERNAL */
1219
1220 lStep = 1; /* Resetto di nuovo lStep */
1221
1222 DEBUGCOUT("Step " << lStep << " has been successfully completed "
1223 "in " << iStIter << " iterations" << std::endl);
1224
1225
1226 DEBUGCOUT("Current time step: " << dCurrTimeStep << std::endl);
1227
1228 pDM->BeforePredict(*pX, *pXPrime, *qX[0], *qXPrime[0]);
1229
1230 Flip();
1231 dRefTimeStep = dInitialTimeStep;
1232 dCurrTimeStep = dRefTimeStep;
1233
1234 CurrStep = StepIntegrator::NEWSTEP;
1235 pTSC->Init(iMaxIterations, dMinTimeStep, MaxTimeStep, dInitialTimeStep);
1236
1237 /* Setup SolutionManager(s) */
1238 ASSERT(pFirstRegularStep!= 0);
1239 SetupSolmans(pFirstRegularStep->GetIntegratorNumUnknownStates(), true);
1240 pCurrStepIntegrator = pFirstRegularStep;
1241
1242 IfFirstStepIsToBeRepeated:
1243 try {
1244 pDM->SetTime(dTime + dCurrTimeStep, dCurrTimeStep, 1);
1245 if (outputStep()) {
1246 if (outputCounter()) {
1247 silent_cout(std::endl);
1248 }
1249 silent_cout("Step(" << 1 << ':' << 0 << ") t=" << dTime + dCurrTimeStep << " dt=" << dCurrTimeStep << std::endl);
1250 }
1251 dTest = pFirstRegularStep->Advance(this, dRefTimeStep,
1252 dCurrTimeStep/dRefTimeStep, CurrStep,
1253 qX, qXPrime, pX, pXPrime,
1254 iStIter, dTest, dSolTest);
1255 }
1256 catch (NonlinearSolver::NoConvergence) {
1257 if (dCurrTimeStep > dMinTimeStep) {
1258 /* Riduce il passo */
1259 CurrStep = StepIntegrator::REPEATSTEP;
1260 doublereal dOldCurrTimeStep = dCurrTimeStep;
1261 dCurrTimeStep = pTSC->dGetNewStepTime(CurrStep, iStIter);
1262 if (dCurrTimeStep < dOldCurrTimeStep) {
1263 DEBUGCOUT("Changing time step"
1264 " from " << dOldCurrTimeStep
1265 << " to " << dCurrTimeStep
1266 << " during first step after "
1267 << iStIter << " iterations"
1268 << std::endl);
1269 goto IfFirstStepIsToBeRepeated;
1270 }
1271 }
1272
1273 silent_cerr("Max iterations number "
1274 << std::abs(pFirstRegularStep->GetIntegratorMaxIters())
1275 << " has been reached during "
1276 "first step, Time=" << dTime << "; "
1277 << "TimeStep=" << dCurrTimeStep
1278 << " cannot be reduced further; "
1279 "aborting..." << std::endl);
1280 pDM->Output(0, dTime, dCurrTimeStep, true);
1281
1282 throw Solver::ErrMaxIterations(MBDYN_EXCEPT_ARGS);
1283 }
1284 catch (NonlinearSolver::ErrSimulationDiverged) {
1285 /*
1286 * Mettere qui eventuali azioni speciali
1287 * da intraprendere in caso di errore ...
1288 */
1289
1290 throw SimulationDiverged(MBDYN_EXCEPT_ARGS);
1291 }
1292 catch (LinearSolver::ErrFactor& err) {
1293 /*
1294 * Mettere qui eventuali azioni speciali
1295 * da intraprendere in caso di errore ...
1296 */
1297 silent_cerr("First step failed because no pivot element "
1298 "could be found for column " << err.iCol
1299 << " (" << pDM->GetDofDescription(err.iCol) << "); "
1300 "aborting..." << std::endl);
1301 throw SimulationDiverged(MBDYN_EXCEPT_ARGS);
1302 }
1303 catch (NonlinearSolver::ConvergenceOnSolution) {
1304 bSolConv = true;
1305 }
1306 catch (EndOfSimulation& eos) {
1307 silent_cerr("Simulation ended during the first regular step:\n"
1308 << eos.what() << "\n");
1309 return false;
1310 }
1311
1312 SAFEDELETE(pFirstRegularStep);
1313 pFirstRegularStep = 0;
1314
1315 bOut = pDM->Output(lStep, dTime + dCurrTimeStep, dCurrTimeStep);
1316
1317 /* Si fa dare l'std::ostream al file di output per il log */
1318 std::ostream& Out = pDM->GetOutFile();
1319
1320 if (mbdyn_stop_at_end_of_time_step()) {
1321 /* Fa l'output della soluzione al primo passo ed esce */
1322 Out << "Interrupted during first step." << std::endl;
1323 throw ErrInterrupted(MBDYN_EXCEPT_ARGS);
1324 }
1325
1326 if (outputMsg()) {
1327 Out
1328 << "Step " << lStep
1329 << " " << dTime + dCurrTimeStep
1330 << " " << dCurrTimeStep
1331 << " " << iStIter
1332 << " " << dTest
1333 << " " << dSolTest
1334 << " " << bSolConv
1335 << " " << bOut
1336 << std::endl;
1337 }
1338
1339 bSolConv = false;
1340
1341 dRefTimeStep = dCurrTimeStep;
1342 dTime += dRefTimeStep;
1343
1344 dTotErr += dTest;
1345 iTotIter += iStIter;
1346
1347 if (EigAn.bAnalysis
1348 && EigAn.currAnalysis != EigAn.Analyses.end()
1349 && *EigAn.currAnalysis <= dTime)
1350 {
1351 std::vector<doublereal>::iterator i = std::find_if(EigAn.Analyses.begin(),
1352 EigAn.Analyses.end(), bind2nd(std::greater<doublereal>(), dTime));
1353 if (i != EigAn.Analyses.end()) {
1354 EigAn.currAnalysis = --i;
1355 }
1356 Eig();
1357 ++EigAn.currAnalysis;
1358 }
1359
1360 if (pRTSolver) {
1361 pRTSolver->Init();
1362 }
1363
1364 /* Altri passi regolari */
1365 ASSERT(pRegularSteps != 0);
1366
1367 /* Setup SolutionManager(s) */
1368 SetupSolmans(pRegularSteps->GetIntegratorNumUnknownStates(), true);
1369 pCurrStepIntegrator = pRegularSteps;
1370
1371 eStatus = SOLVER_STATUS_STARTED;
1372
1373 return true;
1374 }
1375
1376 bool
Advance(void)1377 Solver::Advance(void)
1378 {
1379 DEBUGCOUTFNAME("Solver::Advance");
1380
1381 // consistency check
1382 if (eStatus != SOLVER_STATUS_STARTED) {
1383 silent_cerr("Started() must be called first" << std::endl);
1384 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
1385 }
1386
1387 CurrStep = StepIntegrator::NEWSTEP;
1388
1389 if (pDM->EndOfSimulation() || dTime >= dFinalTime) {
1390 if (pRTSolver) {
1391 pRTSolver->StopCommanded();
1392 }
1393 silent_cout(outputCounterPrefix
1394 << "End of simulation at time "
1395 << dTime << " after "
1396 << lStep << " steps;" << std::endl
1397 << "output in file \"" << sOutputFileName << "\"" << std::endl
1398 << "total iterations: " << iTotIter << std::endl
1399 << "total Jacobian matrices: " << pNLS->TotalAssembledJacobian() << std::endl
1400 << "total error: " << dTotErr << std::endl);
1401
1402 if (pRTSolver) {
1403 pRTSolver->Log();
1404 }
1405
1406 return false;
1407
1408 } else if (pRTSolver && pRTSolver->IsStopCommanded()) {
1409 silent_cout(outputCounterPrefix
1410 << "Simulation is stopped by RTAI task" << std::endl
1411 << "Simulation ended at time "
1412 << dTime << " after "
1413 << lStep << " steps;" << std::endl
1414 << "total iterations: " << iTotIter << std::endl
1415 << "total Jacobian matrices: " << pNLS->TotalAssembledJacobian() << std::endl
1416 << "total error: " << dTotErr << std::endl);
1417 pRTSolver->Log();
1418 return false;
1419
1420 } else if (mbdyn_stop_at_end_of_time_step()
1421 #ifdef USE_MPI
1422 || (MPI_Finalized(&mpi_finalize), mpi_finalize)
1423 #endif /* USE_MPI */
1424 )
1425 {
1426 if (pRTSolver) {
1427 pRTSolver->StopCommanded();
1428 }
1429
1430 silent_cout(outputCounterPrefix
1431 << "Interrupted!" << std::endl
1432 << "Simulation ended at time "
1433 << dTime << " after "
1434 << lStep << " steps;" << std::endl
1435 << "total iterations: " << iTotIter << std::endl
1436 << "total Jacobian matrices: " << pNLS->TotalAssembledJacobian() << std::endl
1437 << "total error: " << dTotErr << std::endl);
1438
1439 if (pRTSolver) {
1440 pRTSolver->Log();
1441 }
1442
1443 throw ErrInterrupted(MBDYN_EXCEPT_ARGS);
1444 }
1445
1446 lStep++;
1447 pDM->BeforePredict(*pX, *pXPrime, *qX[0], *qXPrime[0]);
1448
1449 Flip();
1450
1451 if (pRTSolver) {
1452 pRTSolver->Wait();
1453 }
1454
1455 int retries = -1;
1456 IfStepIsToBeRepeated:
1457 try {
1458 retries++;
1459 pDM->SetTime(dTime + dCurrTimeStep, dCurrTimeStep, lStep);
1460 if (outputStep()) {
1461 if (outputCounter()) {
1462 silent_cout(std::endl);
1463 }
1464 silent_cout("Step(" << lStep << ':' << retries << ") t=" << dTime + dCurrTimeStep << " dt=" << dCurrTimeStep << std::endl);
1465 }
1466 dTest = pRegularSteps->Advance(this, dRefTimeStep,
1467 dCurrTimeStep/dRefTimeStep, CurrStep,
1468 qX, qXPrime, pX, pXPrime, iStIter,
1469 dTest, dSolTest);
1470 }
1471 catch (NonlinearSolver::NoConvergence) {
1472 if (dCurrTimeStep > dMinTimeStep) {
1473 /* Riduce il passo */
1474 CurrStep = StepIntegrator::REPEATSTEP;
1475 doublereal dOldCurrTimeStep = dCurrTimeStep;
1476 dCurrTimeStep = pTSC->dGetNewStepTime(CurrStep, iStIter);
1477 if (dCurrTimeStep < dOldCurrTimeStep) {
1478 DEBUGCOUT("Changing time step"
1479 " from " << dOldCurrTimeStep
1480 << " to " << dCurrTimeStep
1481 << " during step "
1482 << lStep << " after "
1483 << iStIter << " iterations"
1484 << std::endl);
1485 goto IfStepIsToBeRepeated;
1486 }
1487 }
1488
1489 silent_cerr(outputCounterPrefix
1490 << "Max iterations number "
1491 << std::abs(pRegularSteps->GetIntegratorMaxIters())
1492 << " has been reached during "
1493 "Step=" << lStep << ", "
1494 "Time=" << dTime + dCurrTimeStep << "; "
1495 "TimeStep=" << dCurrTimeStep
1496 << " cannot be reduced further; "
1497 "aborting..." << std::endl);
1498 throw ErrMaxIterations(MBDYN_EXCEPT_ARGS);
1499 }
1500 catch (NonlinearSolver::ErrSimulationDiverged) {
1501 if (dCurrTimeStep > dMinTimeStep) {
1502 /* Riduce il passo */
1503 CurrStep = StepIntegrator::REPEATSTEP;
1504 doublereal dOldCurrTimeStep = dCurrTimeStep;
1505 dCurrTimeStep = pTSC->dGetNewStepTime(CurrStep, iStIter);
1506 if (dCurrTimeStep < dOldCurrTimeStep) {
1507 DEBUGCOUT("Changing time step"
1508 " from " << dOldCurrTimeStep
1509 << " to " << dCurrTimeStep
1510 << " during step "
1511 << lStep << " after "
1512 << iStIter << " iterations"
1513 << std::endl);
1514 goto IfStepIsToBeRepeated;
1515 }
1516 }
1517
1518 silent_cerr(outputCounterPrefix
1519 << "Simulation diverged after "
1520 << iStIter << " iterations, before "
1521 "reaching max iteration number "
1522 << std::abs(pRegularSteps->GetIntegratorMaxIters())
1523 << " during Step=" << lStep << ", "
1524 "Time=" << dTime + dCurrTimeStep << "; "
1525 "TimeStep=" << dCurrTimeStep
1526 << " cannot be reduced further; "
1527 "aborting..." << std::endl);
1528 throw SimulationDiverged(MBDYN_EXCEPT_ARGS);
1529 }
1530 catch (LinearSolver::ErrFactor& err) {
1531 /*
1532 * Mettere qui eventuali azioni speciali
1533 * da intraprendere in caso di errore ...
1534 */
1535 silent_cerr(outputCounterPrefix
1536 << "Simulation failed because no pivot element "
1537 "could be found for column " << err.iCol
1538 << " (" << pDM->GetDofDescription(err.iCol) << ") "
1539 "after " << iStIter << " iterations "
1540 "during step " << lStep << "; "
1541 "aborting..." << std::endl);
1542 throw SimulationDiverged(MBDYN_EXCEPT_ARGS);
1543 }
1544 catch (NonlinearSolver::ConvergenceOnSolution) {
1545 bSolConv = true;
1546 }
1547 catch (EndOfSimulation& eos) {
1548 silent_cerr(outputCounterPrefix
1549 << "Simulation ended during a regular step:\n"
1550 << eos.what() << "\n");
1551 #ifdef USE_MPI
1552 MBDynComm.Abort(0);
1553 #endif /* USE_MPI */
1554 if (pRTSolver) {
1555 pRTSolver->StopCommanded();
1556 }
1557
1558 silent_cout("Simulation ended at time "
1559 << dTime << " after "
1560 << lStep << " steps;" << std::endl
1561 << "total iterations: " << iTotIter << std::endl
1562 << "total Jacobian matrices: " << pNLS->TotalAssembledJacobian() << std::endl
1563 << "total error: " << dTotErr << std::endl);
1564
1565 if (pRTSolver) {
1566 pRTSolver->Log();
1567 }
1568
1569 return false;
1570 }
1571
1572 dTotErr += dTest;
1573 iTotIter += iStIter;
1574
1575 bOut = pDM->Output(lStep, dTime + dCurrTimeStep, dCurrTimeStep);
1576
1577 /* Si fa dare l'std::ostream al file di output per il log */
1578 std::ostream& Out = pDM->GetOutFile();
1579
1580 if (outputMsg()) {
1581 Out << "Step " << lStep
1582 << " " << dTime + dCurrTimeStep
1583 << " " << dCurrTimeStep
1584 << " " << iStIter
1585 << " " << dTest
1586 << " " << dSolTest
1587 << " " << bSolConv
1588 << " " << bOut
1589 << std::endl;
1590 }
1591
1592 if (bOutputCounter) {
1593 silent_cout("Step " << std::setw(5) << lStep
1594 << " " << std::setw(13) << dTime + dCurrTimeStep
1595 << " " << std::setw(13) << dCurrTimeStep
1596 << " " << std::setw(4) << iStIter
1597 << " " << std::setw(13) << dTest
1598 << " " << std::setw(13) << dSolTest
1599 << " " << bSolConv
1600 << " " << bOut
1601 << outputCounterPostfix);
1602 }
1603
1604 DEBUGCOUT("Step " << lStep
1605 << " has been successfully completed "
1606 "in " << iStIter << " iterations" << std::endl);
1607
1608 dRefTimeStep = dCurrTimeStep;
1609 dTime += dRefTimeStep;
1610
1611 bSolConv = false;
1612
1613 if (EigAn.bAnalysis
1614 && EigAn.currAnalysis != EigAn.Analyses.end()
1615 && *EigAn.currAnalysis <= dTime)
1616 {
1617 std::vector<doublereal>::iterator i = std::find_if(EigAn.Analyses.begin(),
1618 EigAn.Analyses.end(), bind2nd(std::greater<doublereal>(), dTime));
1619 if (i != EigAn.Analyses.end()) {
1620 EigAn.currAnalysis = --i;
1621 }
1622 Eig(bOutputCounter);
1623 ++EigAn.currAnalysis;
1624 }
1625
1626 /* Calcola il nuovo timestep */
1627 dCurrTimeStep = pTSC->dGetNewStepTime(CurrStep, iStIter);
1628 DEBUGCOUT("Current time step: " << dCurrTimeStep << std::endl);
1629
1630 return true;
1631 }
1632
1633 void
Run(void)1634 Solver::Run(void)
1635 {
1636 DEBUGCOUTFNAME("Solver::Run");
1637
1638 if (Prepare()) {
1639 if (Start()) {
1640 while (Advance()) {
1641 NO_OP;
1642 }
1643 }
1644 }
1645 }
1646
1647 /* Distruttore */
~Solver(void)1648 Solver::~Solver(void)
1649 {
1650 DEBUGCOUTFNAME("Solver::~Solver");
1651
1652 if (!qX.empty()) {
1653 for (int ivec = 0; ivec < iNumPreviousVectors; ivec++) {
1654 if (qX[ivec] != 0) {
1655 SAFEDELETE(qX[ivec]);
1656 SAFEDELETE(qXPrime[ivec]);
1657 }
1658 }
1659 }
1660
1661 if (pX) {
1662 SAFEDELETE(pX);
1663 }
1664
1665 if (pXPrime) {
1666 SAFEDELETE(pXPrime);
1667 }
1668
1669 if (pdWorkSpace) {
1670 SAFEDELETEARR(pdWorkSpace);
1671 }
1672
1673 if (pDM) {
1674 SAFEDELETE(pDM);
1675 }
1676
1677 if (pRTSolver) {
1678 SAFEDELETE(pRTSolver);
1679 }
1680
1681 if (pDerivativeSteps) {
1682 SAFEDELETE(pDerivativeSteps);
1683 }
1684
1685 if (pFirstDummyStep) {
1686 SAFEDELETE(pFirstDummyStep);
1687 }
1688
1689 if (pDummySteps) {
1690 SAFEDELETE(pDummySteps);
1691 }
1692
1693 if (pFirstRegularStep) {
1694 SAFEDELETE(pFirstRegularStep);
1695 }
1696
1697 if (pRegularSteps) {
1698 SAFEDELETE(pRegularSteps);
1699 }
1700
1701 if (pSM) {
1702 SAFEDELETE(pSM);
1703 }
1704
1705 if (pNLS) {
1706 SAFEDELETE(pNLS);
1707 }
1708
1709 if (pTSC) {
1710 delete pTSC;
1711 }
1712
1713 DestroyTimeStepData();
1714 }
1715
1716 /*scrive il contributo al file di restart*/
1717 std::ostream &
Restart(std::ostream & out,DataManager::eRestart type) const1718 Solver::Restart(std::ostream& out,DataManager::eRestart type) const
1719 {
1720
1721 out << "begin: initial value;" << std::endl;
1722 switch(type) {
1723 case DataManager::ATEND:
1724 out << " # initial time: " << pDM->dGetTime() << ";"
1725 << std::endl
1726 << " # final time: " << dFinalTime << ";"
1727 << std::endl
1728 << " # time step: " << dInitialTimeStep << ";"
1729 << std::endl;
1730 break;
1731 case DataManager::ITERATIONS:
1732 case DataManager::TIME:
1733 case DataManager::TIMES:
1734 out << " initial time: " << pDM->dGetTime()<< ";" << std::endl
1735 << " final time: " << dFinalTime << ";" << std::endl
1736 << " time step: " << dInitialTimeStep << ";"
1737 << std::endl;
1738 break;
1739 default:
1740 ASSERT(0);
1741 }
1742
1743 out << " method: ";
1744 switch(RegularType) {
1745 case INT_CRANKNICOLSON:
1746 out << "Crank Nicolson; " << std::endl;
1747 break;
1748 case INT_MS2:
1749 out << "ms, ";
1750 pRhoRegular->Restart(out) << ", ";
1751 pRhoAlgebraicRegular->Restart(out) << ";" << std::endl;
1752 break;
1753 case INT_HOPE:
1754 out << "hope, ";
1755 pRhoRegular->Restart(out) << ", ";
1756 pRhoAlgebraicRegular->Restart(out) << ";" << std::endl;
1757 break;
1758
1759 case INT_THIRDORDER:
1760 out << "thirdorder, ";
1761 if (!pRhoRegular)
1762 out << "ad hoc;" << std::endl;
1763 else
1764 pRhoRegular->Restart(out) << ";" << std::endl;
1765 break;
1766 case INT_IMPLICITEULER:
1767 out << "implicit euler;" << std::endl;
1768 break;
1769 default:
1770 ASSERT(0);
1771 }
1772
1773 integer iMI = pRegularSteps->GetIntegratorMaxIters();
1774 out << " max iterations: " << std::abs(iMI);
1775 if (iMI < 0) out << ", at most";
1776 out
1777 << ";" << std::endl
1778 << " tolerance: " << pRegularSteps->GetIntegratorDTol();
1779 switch(ResTest) {
1780 case NonlinearSolverTest::NORM:
1781 out << ", test, norm" ;
1782 break;
1783 case NonlinearSolverTest::MINMAX:
1784 out << ", test, minmax" ;
1785 break;
1786 case NonlinearSolverTest::NONE:
1787 NO_OP;
1788 default:
1789 silent_cerr("unhandled nonlinear solver test type" << std::endl);
1790 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
1791 }
1792
1793 if (bScale) {
1794 out << ", scale"
1795 << ", " << pRegularSteps->GetIntegratorDSolTol();
1796 }
1797
1798 switch (SolTest) {
1799 case NonlinearSolverTest::NORM:
1800 out << ", test, norm" ;
1801 break;
1802 case NonlinearSolverTest::MINMAX:
1803 out << ", test, minmax" ;
1804 break;
1805 case NonlinearSolverTest::NONE:
1806 NO_OP;
1807 default:
1808 ASSERT(0);
1809 }
1810 out
1811 << ";" << std::endl;
1812
1813 if ( pDerivativeSteps != 0 ) {
1814 out
1815 << " derivatives max iterations: " << pDerivativeSteps->GetIntegratorMaxIters() << ";" << std::endl
1816 << " derivatives tolerance: " << pDerivativeSteps->GetIntegratorDTol() << ";" << std::endl
1817 << " derivatives coefficient: " << dDerivativesCoef << ";" << std::endl;
1818 }
1819
1820 if (iDummyStepsNumber) {
1821 out
1822 << " dummy steps max iterations: " << pDummySteps->GetIntegratorMaxIters() << ";" << std::endl
1823 << " dummy steps tolerance: " << pDummySteps->GetIntegratorDTol() << ";" << std::endl;
1824 }
1825 out
1826 << " dummy steps number: " << iDummyStepsNumber << ";" << std::endl
1827 << " dummy steps ratio: " << dDummyStepsRatio << ";" << std::endl;
1828 switch (NonlinearSolverType) {
1829 case NonlinearSolver::MATRIXFREE:
1830 out << " # nonlinear solver: matrix free;" << std::endl;
1831 break;
1832 case NonlinearSolver::NEWTONRAPHSON:
1833 default :
1834 out << " nonlinear solver: newton raphson";
1835 if (!bTrueNewtonRaphson) {
1836 out << ", modified, " << iIterationsBeforeAssembly;
1837 if (bKeepJac) {
1838 out << ", keep jacobian matrix";
1839 }
1840 if (bHonorJacRequest) {
1841 out << ", honor element requests";
1842 }
1843 }
1844 out << ";" << std::endl;
1845 }
1846 out << " solver: ";
1847 RestartLinSol(out, CurrLinearSolver);
1848 out << "end: initial value;" << std::endl << std::endl;
1849 return out;
1850 }
1851
1852 /* Dati dell'integratore */
1853 void
ReadData(MBDynParser & HP)1854 Solver::ReadData(MBDynParser& HP)
1855 {
1856 DEBUGCOUTFNAME("MultiStepIntegrator::ReadData");
1857
1858 /* parole chiave */
1859 static const char*const sKeyWords[] = {
1860 "begin",
1861 "initial" "value",
1862 "multistep", /* deprecated */
1863 "end",
1864
1865 "initial" "time",
1866 "final" "time",
1867 "time" "step",
1868 "min" "time" "step",
1869 "max" "time" "step",
1870 "tolerance",
1871 "max" "residual",
1872 "max" "iterations",
1873 "modify" "residual" "test",
1874 "enforce" "constraint" "equations",
1875 /* DEPRECATED */
1876 "fictitious" "steps" "number",
1877 "fictitious" "steps" "ratio",
1878 "fictitious" "steps" "tolerance",
1879 "fictitious" "steps" "max" "iterations",
1880 /* END OF DEPRECATED */
1881
1882 "dummy" "steps" "number",
1883 "dummy" "steps" "ratio",
1884 "dummy" "steps" "tolerance",
1885 "dummy" "steps" "max" "iterations",
1886
1887 "abort" "after",
1888 "input",
1889 "assembly",
1890 "derivatives",
1891
1892 /* DEPRECATED */ "fictitious" "steps" /* END OF DEPRECATED */ ,
1893 "dummy" "steps",
1894
1895 "output",
1896 "none",
1897 "iterations",
1898 "residual",
1899 "solution",
1900 /* DEPRECATED */ "jacobian" /* END OF DEPRECATED */ ,
1901 "jacobian" "matrix",
1902 "bailout",
1903 "messages",
1904 "counter",
1905 "matrix" "condition" "number",
1906 "solver" "condition" "number",
1907 "cpu" "time",
1908 "output" "meter",
1909
1910 "method",
1911 /* DEPRECATED */ "fictitious" "steps" "method" /* END OF DEPRECATED */ ,
1912 "dummy" "steps" "method",
1913
1914 "Crank" "Nicolson",
1915 /* DEPRECATED */ "Crank" "Nicholson" /* END OF DEPRECATED */ ,
1916 /* DEPRECATED */ "nostro" /* END OF DEPRECATED */ ,
1917 "ms",
1918 "hope",
1919 "bdf",
1920 "thirdorder",
1921 "implicit" "euler",
1922
1923 "derivatives" "coefficient",
1924 "derivatives" "tolerance",
1925 "derivatives" "max" "iterations",
1926
1927 /* DEPRECATED */
1928 "true",
1929 "modified",
1930 /* END OF DEPRECATED */
1931
1932 "strategy",
1933 "factor",
1934 "no" "change",
1935 "change",
1936
1937 "pod",
1938 "eigen" "analysis",
1939
1940 /* DEPRECATED */
1941 "solver",
1942 "interface" "solver",
1943 /* END OF DEPRECATED */
1944 "linear" "solver",
1945 "interface" "linear" "solver",
1946
1947 /* DEPRECATED */
1948 "preconditioner",
1949 /* END OF DEPRECATED */
1950
1951 "nonlinear" "solver",
1952 "default",
1953 "newton" "raphson",
1954 "line" "search",
1955 "matrix" "free",
1956 "bicgstab",
1957 "gmres",
1958 /* DEPRECATED */ "full" "jacobian" /* END OF DEPRECATED */ ,
1959 "full" "jacobian" "matrix",
1960
1961 /* RTAI stuff */
1962 "real" "time",
1963
1964 /* multithread stuff */
1965 "threads",
1966
1967 NULL
1968 };
1969
1970 /* enum delle parole chiave */
1971 enum KeyWords {
1972 UNKNOWN = -1,
1973 BEGIN = 0,
1974 INITIAL_VALUE,
1975 MULTISTEP,
1976 END,
1977
1978 INITIALTIME,
1979 FINALTIME,
1980 TIMESTEP,
1981 MINTIMESTEP,
1982 MAXTIMESTEP,
1983 TOLERANCE,
1984 MAXRESIDUAL,
1985 MAXITERATIONS,
1986 MODIFY_RES_TEST,
1987 ENFORCE_CONSTRAINT_EQUATIONS,
1988 FICTITIOUSSTEPSNUMBER,
1989 FICTITIOUSSTEPSRATIO,
1990 FICTITIOUSSTEPSTOLERANCE,
1991 FICTITIOUSSTEPSMAXITERATIONS,
1992
1993 DUMMYSTEPSNUMBER,
1994 DUMMYSTEPSRATIO,
1995 DUMMYSTEPSTOLERANCE,
1996 DUMMYSTEPSMAXITERATIONS,
1997
1998 ABORTAFTER,
1999 INPUT,
2000 ASSEMBLY,
2001 DERIVATIVES,
2002 FICTITIOUSSTEPS,
2003 DUMMYSTEPS,
2004
2005 OUTPUT,
2006 NONE,
2007 ITERATIONS,
2008 RESIDUAL,
2009 SOLUTION,
2010 JACOBIAN,
2011 JACOBIANMATRIX,
2012 BAILOUT,
2013 MESSAGES,
2014 COUNTER,
2015 MATRIX_COND_NUM,
2016 SOLVER_COND_NUM,
2017 CPU_TIME,
2018 OUTPUTMETER,
2019
2020 METHOD,
2021 FICTITIOUSSTEPSMETHOD,
2022 DUMMYSTEPSMETHOD,
2023 CRANKNICOLSON,
2024 CRANKNICHOLSON,
2025 NOSTRO,
2026 MS,
2027 HOPE,
2028 BDF,
2029 THIRDORDER,
2030 IMPLICITEULER,
2031
2032 DERIVATIVESCOEFFICIENT,
2033 DERIVATIVESTOLERANCE,
2034 DERIVATIVESMAXITERATIONS,
2035
2036 /* DEPRECATED */
2037 NR_TRUE,
2038 MODIFIED,
2039 /* END OF DEPRECATED */
2040
2041 STRATEGY,
2042 STRATEGYFACTOR,
2043 STRATEGYNOCHANGE,
2044 STRATEGYCHANGE,
2045
2046 POD,
2047 EIGENANALYSIS,
2048
2049 /* DEPRECATED */
2050 SOLVER,
2051 INTERFACESOLVER,
2052 /* END OF DEPRECATED */
2053 LINEARSOLVER,
2054 INTERFACELINEARSOLVER,
2055
2056 /* DEPRECATED */
2057 PRECONDITIONER,
2058 /* END OF DEPRECATED */
2059
2060 NONLINEARSOLVER,
2061 DEFAULT,
2062 NEWTONRAPHSON,
2063 LINESEARCH,
2064 MATRIXFREE,
2065 BICGSTAB,
2066 GMRES,
2067 FULLJACOBIAN,
2068 FULLJACOBIANMATRIX,
2069
2070 /* RTAI stuff */
2071 REALTIME,
2072
2073 THREADS,
2074
2075 LASTKEYWORD
2076 };
2077
2078 /* tabella delle parole chiave */
2079 KeyTable K(HP, sKeyWords);
2080
2081 /* legge i dati della simulazione */
2082 if (KeyWords(HP.GetDescription()) != BEGIN) {
2083 silent_cerr("Error: <begin> expected at line "
2084 << HP.GetLineData() << "; aborting..." << std::endl);
2085 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2086 }
2087
2088 switch (KeyWords(HP.GetWord())) {
2089 case MULTISTEP:
2090 pedantic_cout("warning: \"begin: multistep\" is deprecated; "
2091 "use \"begin: initial value;\" instead." << std::endl);
2092 case INITIAL_VALUE:
2093 break;
2094
2095 default:
2096 silent_cerr("Error: \"begin: initial value;\" expected at line "
2097 << HP.GetLineData() << "; aborting..." << std::endl);
2098 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2099 }
2100
2101 bool bMethod(false);
2102 bool bDummyStepsMethod(false);
2103
2104 /* dati letti qui ma da passare alle classi
2105 * StepIntegration e NonlinearSolver
2106 */
2107
2108 doublereal dTol = ::dDefaultTol;
2109 doublereal dSolutionTol = 0.;
2110 iMaxIterations = ::iDefaultMaxIterations;
2111 bool bModResTest = false;
2112 bool bSetScaleAlgebraic = false;
2113
2114 /* Dati dei passi fittizi di trimmaggio iniziale */
2115 doublereal dDummyStepsTolerance = ::dDefaultDummyStepsTolerance;
2116 integer iDummyStepsMaxIterations = ::iDefaultMaxIterations;
2117
2118 /* Dati del passo iniziale di calcolo delle derivate */
2119
2120 doublereal dDerivativesTol = ::dDefaultTol;
2121 integer iDerivativesMaxIterations = ::iDefaultMaxIterations;
2122 integer iDerivativesCoefMaxIter = 0;
2123 doublereal dDerivativesCoefFactor = 10.;
2124
2125 #ifdef USE_MULTITHREAD
2126 bool bSolverThreads(false);
2127 unsigned nSolverThreads = 0;
2128 #endif /* USE_MULTITHREAD */
2129
2130
2131 /* Ciclo infinito */
2132 while (true) {
2133 KeyWords CurrKeyWord = KeyWords(HP.GetDescription());
2134
2135 switch (CurrKeyWord) {
2136 case INITIALTIME:
2137 dInitialTime = HP.GetReal();
2138 DEBUGLCOUT(MYDEBUG_INPUT, "Initial time is "
2139 << dInitialTime << std::endl);
2140 break;
2141
2142 case FINALTIME:
2143 if (HP.IsKeyWord("forever")) {
2144 dFinalTime = std::numeric_limits<doublereal>::max();
2145 } else {
2146 dFinalTime = HP.GetReal();
2147 }
2148 DEBUGLCOUT(MYDEBUG_INPUT, "Final time is "
2149 << dFinalTime << std::endl);
2150
2151 if(dFinalTime <= dInitialTime) {
2152 silent_cerr("warning: final time " << dFinalTime
2153 << " is less than initial time "
2154 << dInitialTime << ';' << std::endl
2155 << "this will cause the simulation"
2156 " to abort" << std::endl);
2157 }
2158 break;
2159
2160 case TIMESTEP:
2161 dInitialTimeStep = HP.GetReal();
2162 DEBUGLCOUT(MYDEBUG_INPUT, "Initial time step is "
2163 << dInitialTimeStep << std::endl);
2164
2165 if (dInitialTimeStep == 0.) {
2166 silent_cerr("warning, null initial time step"
2167 " is not allowed" << std::endl);
2168 } else if (dInitialTimeStep < 0.) {
2169 dInitialTimeStep = -dInitialTimeStep;
2170 silent_cerr("warning, negative initial time step"
2171 " is not allowed;" << std::endl
2172 << "its modulus " << dInitialTimeStep
2173 << " will be considered" << std::endl);
2174 }
2175 break;
2176
2177 case MINTIMESTEP:
2178 try {
2179 dMinTimeStep = HP.GetReal(std::numeric_limits<doublereal>::min(), HighParser::range_gt<doublereal>(0.));
2180
2181 } catch (HighParser::ErrValueOutOfRange<doublereal> e) {
2182 silent_cerr("invalid min time step " << e.Get() << " (must be positive) [" << e.what() << "] at line " << HP.GetLineData() << std::endl);
2183 throw e;
2184 }
2185 break;
2186
2187 case MAXTIMESTEP:
2188 if (HP.IsKeyWord("unlimited")) {
2189 DriveCaller *pDC = 0;
2190 SAFENEWWITHCONSTRUCTOR(pDC, ConstDriveCaller, ConstDriveCaller(0.));
2191 MaxTimeStep.Set(pDC);
2192
2193 } else {
2194 MaxTimeStep.Set(HP.GetDriveCaller());
2195 }
2196
2197 #ifdef DEBUG
2198 if (typeid(*MaxTimeStep.pGetDriveCaller()) == typeid(PostponedDriveCaller)) {
2199 DEBUGLCOUT(MYDEBUG_INPUT, "Max time step is postponed" << std::endl);
2200
2201 } else {
2202 DEBUGLCOUT(MYDEBUG_INPUT, "Max time step is " << MaxTimeStep.dGet() << std::endl);
2203 }
2204 #endif // DEBUG
2205
2206 eTimeStepLimit = (HP.IsKeyWord("hard" "limit")
2207 && HP.GetYesNoOrBool())
2208 ? TS_HARD_LIMIT
2209 : TS_SOFT_LIMIT;
2210
2211 if (dGetInitialMaxTimeStep() == 0.) {
2212 silent_cout("no max time step limit will be"
2213 " considered" << std::endl);
2214
2215 } else if (dGetInitialMaxTimeStep() < 0.) {
2216 silent_cerr("negative max time step"
2217 " is not allowed" << std::endl);
2218 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2219 }
2220 break;
2221
2222 case FICTITIOUSSTEPSNUMBER:
2223 case DUMMYSTEPSNUMBER:
2224 iDummyStepsNumber = HP.GetInt();
2225 if (iDummyStepsNumber < 0) {
2226 iDummyStepsNumber =
2227 ::iDefaultDummyStepsNumber;
2228 silent_cerr("negative dummy steps number"
2229 " is illegal" << std::endl);
2230 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2231
2232 } else if (iDummyStepsNumber == 1) {
2233 silent_cerr("warning, a single dummy step"
2234 " may be useless" << std::endl);
2235 }
2236
2237 DEBUGLCOUT(MYDEBUG_INPUT, "Dummy steps number: "
2238 << iDummyStepsNumber << std::endl);
2239 break;
2240
2241 case FICTITIOUSSTEPSRATIO:
2242 case DUMMYSTEPSRATIO:
2243 dDummyStepsRatio = HP.GetReal();
2244 if (dDummyStepsRatio < 0.) {
2245 silent_cerr("negative dummy steps ratio"
2246 " is illegal" << std::endl);
2247 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2248 }
2249
2250 if (dDummyStepsRatio > 1.) {
2251 silent_cerr("warning, dummy steps ratio"
2252 " is larger than one." << std::endl
2253 << "Something like 1.e-3 should"
2254 " be safer ..." << std::endl);
2255 }
2256
2257 DEBUGLCOUT(MYDEBUG_INPUT, "Dummy steps ratio: "
2258 << dDummyStepsRatio << std::endl);
2259 break;
2260
2261 case FICTITIOUSSTEPSTOLERANCE:
2262 case DUMMYSTEPSTOLERANCE:
2263 dDummyStepsTolerance = HP.GetReal();
2264 if (dDummyStepsTolerance <= 0.) {
2265 silent_cerr("negative dummy steps"
2266 " tolerance is illegal" << std::endl);
2267 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2268 }
2269 DEBUGLCOUT(MYDEBUG_INPUT,
2270 "Dummy steps tolerance: "
2271 << dDummyStepsTolerance << std::endl);
2272 break;
2273
2274 case ABORTAFTER: {
2275 KeyWords WhenToAbort(KeyWords(HP.GetWord()));
2276 switch (WhenToAbort) {
2277 case INPUT:
2278 eAbortAfter = AFTER_INPUT;
2279 DEBUGLCOUT(MYDEBUG_INPUT,
2280 "Simulation will abort after"
2281 " data input" << std::endl);
2282 break;
2283
2284 case ASSEMBLY:
2285 eAbortAfter = AFTER_ASSEMBLY;
2286 DEBUGLCOUT(MYDEBUG_INPUT,
2287 "Simulation will abort after"
2288 " initial assembly" << std::endl);
2289 break;
2290
2291 case DERIVATIVES:
2292 eAbortAfter = AFTER_DERIVATIVES;
2293 DEBUGLCOUT(MYDEBUG_INPUT,
2294 "Simulation will abort after"
2295 " derivatives solution" << std::endl);
2296 break;
2297
2298 case FICTITIOUSSTEPS:
2299 case DUMMYSTEPS:
2300 eAbortAfter = AFTER_DUMMY_STEPS;
2301 DEBUGLCOUT(MYDEBUG_INPUT,
2302 "Simulation will abort after"
2303 " dummy steps solution" << std::endl);
2304 break;
2305
2306 default:
2307 silent_cerr("Don't know when to abort,"
2308 " so I'm going to abort now" << std::endl);
2309 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2310 }
2311 } break;
2312
2313 case OUTPUT: {
2314 unsigned OF = OUTPUT_DEFAULT;
2315 bool setOutput = false;
2316
2317 while (HP.IsArg()) {
2318 KeyWords OutputFlag(KeyWords(HP.GetWord()));
2319 switch (OutputFlag) {
2320 case NONE:
2321 OF = OUTPUT_NONE;
2322 setOutput = true;
2323 break;
2324
2325 case ITERATIONS:
2326 OF |= OUTPUT_ITERS;
2327 break;
2328
2329 case RESIDUAL:
2330 OF |= OUTPUT_RES;
2331 break;
2332
2333 case SOLUTION:
2334 OF |= OUTPUT_SOL;
2335 break;
2336
2337 case JACOBIAN:
2338 case JACOBIANMATRIX:
2339 OF |= OUTPUT_JAC;
2340 break;
2341
2342 case BAILOUT:
2343 OF |= OUTPUT_BAILOUT;
2344 break;
2345
2346 case MESSAGES:
2347 OF |= OUTPUT_MSG;
2348 break;
2349
2350 case COUNTER:
2351 OF |= OUTPUT_COUNTER;
2352 break;
2353
2354 case MATRIX_COND_NUM:
2355 DelOutputFlags(OUTPUT_MAT_COND_NUM);
2356 if (HP.IsKeyWord("norm")) {
2357 if (HP.IsKeyWord("inf")) {
2358 OF |= OUTPUT_MAT_COND_NUM_INF;
2359 } else {
2360 const double P = HP.GetReal();
2361 if (P != 1.) {
2362 silent_cerr("Only inf or 1 norm are supported for condition numbers at line " << HP.GetLineData() << std::endl);
2363 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2364 }
2365 OF |= OUTPUT_MAT_COND_NUM_1;
2366 }
2367 } else {
2368 OF |= OUTPUT_MAT_COND_NUM_1;
2369 }
2370 break;
2371
2372 case SOLVER_COND_NUM:
2373 OF |= OUTPUT_SOLVER_COND_NUM;
2374 if (HP.IsKeyWord("stat")) {
2375 if (HP.GetYesNoOrBool()) {
2376 OF |= OUTPUT_SOLVER_COND_STAT;
2377 } else {
2378 DelOutputFlags(OUTPUT_SOLVER_COND_STAT);
2379 }
2380 }
2381 break;
2382 case CPU_TIME:
2383 OF |= OUTPUT_CPU_TIME;
2384 break;
2385
2386 default:
2387 silent_cerr("Warning, unknown output flag "
2388 "at line " << HP.GetLineData()
2389 << "; ignored" << std::endl);
2390 break;
2391 }
2392 }
2393
2394 if (setOutput) {
2395 SetOutputFlags(OF);
2396
2397 } else {
2398 AddOutputFlags(OF);
2399 }
2400 } break;
2401
2402 case OUTPUTMETER:
2403 SetOutputMeter(HP.GetDriveCaller(true));
2404 break;
2405
2406 case METHOD: {
2407 if (bMethod) {
2408 silent_cerr("error: multiple definition"
2409 " of integration method at line "
2410 << HP.GetLineData() << std::endl);
2411 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2412 }
2413 bMethod = true;
2414
2415 KeyWords KMethod = KeyWords(HP.GetWord());
2416 switch (KMethod) {
2417 case CRANKNICHOLSON:
2418 silent_cout("warning: \"crank nicolson\" is the correct spelling "
2419 "at line " << HP.GetLineData() << std::endl);
2420 case CRANKNICOLSON:
2421 RegularType = INT_CRANKNICOLSON;
2422 break;
2423
2424 case BDF:
2425 /* default (order 2) */
2426 RegularType = INT_MS2;
2427
2428 if (HP.IsKeyWord("order")) {
2429 int iOrder = HP.GetInt();
2430
2431 switch (iOrder) {
2432 case 1:
2433 RegularType = INT_IMPLICITEULER;
2434 break;
2435
2436 case 2:
2437 break;
2438
2439 default:
2440 silent_cerr("unhandled BDF order " << iOrder << std::endl);
2441 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2442 }
2443 }
2444
2445 if (RegularType == INT_MS2) {
2446 SAFENEW(pRhoRegular, NullDriveCaller);
2447 SAFENEW(pRhoAlgebraicRegular, NullDriveCaller);
2448 }
2449 break;
2450
2451 case NOSTRO:
2452 silent_cerr("integration method \"nostro\" "
2453 "is deprecated; use \"ms\" "
2454 "instead at line "
2455 << HP.GetLineData()
2456 << std::endl);
2457 case MS:
2458 case HOPE:
2459 pRhoRegular = HP.GetDriveCaller(true);
2460
2461 pRhoAlgebraicRegular = 0;
2462 if (HP.IsArg()) {
2463 pRhoAlgebraicRegular = HP.GetDriveCaller(true);
2464 } else {
2465 pRhoAlgebraicRegular = pRhoRegular->pCopy();
2466 }
2467
2468 switch (KMethod) {
2469 case NOSTRO:
2470 case MS:
2471 RegularType = INT_MS2;
2472 break;
2473
2474 case HOPE:
2475 RegularType = INT_HOPE;
2476 break;
2477
2478 default:
2479 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2480 }
2481 break;
2482
2483 case THIRDORDER:
2484 if (HP.IsKeyWord("ad" "hoc")) {
2485 /* do nothing */ ;
2486 } else {
2487 pRhoRegular = HP.GetDriveCaller(true);
2488 }
2489 RegularType = INT_THIRDORDER;
2490 break;
2491
2492 case IMPLICITEULER:
2493 RegularType = INT_IMPLICITEULER;
2494 break;
2495
2496 default:
2497 silent_cerr("Unknown integration method at line "
2498 << HP.GetLineData() << std::endl);
2499 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2500 }
2501 break;
2502 }
2503
2504 case FICTITIOUSSTEPSMETHOD:
2505 case DUMMYSTEPSMETHOD: {
2506 if (bDummyStepsMethod) {
2507 silent_cerr("error: multiple definition "
2508 "of dummy steps integration method "
2509 "at line " << HP.GetLineData()
2510 << std::endl);
2511 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2512 }
2513 bDummyStepsMethod = true;
2514
2515 KeyWords KMethod = KeyWords(HP.GetWord());
2516 switch (KMethod) {
2517 case CRANKNICHOLSON:
2518 silent_cout("warning: \"crank nicolson\" is the correct spelling" << std::endl);
2519 case CRANKNICOLSON:
2520 DummyType = INT_CRANKNICOLSON;
2521 break;
2522
2523 case BDF:
2524 /* default (order 2) */
2525 DummyType = INT_MS2;
2526
2527 if (HP.IsKeyWord("order")) {
2528 int iOrder = HP.GetInt();
2529
2530 switch (iOrder) {
2531 case 1:
2532 DummyType = INT_IMPLICITEULER;
2533 break;
2534
2535 case 2:
2536 break;
2537
2538 default:
2539 silent_cerr("unhandled BDF order " << iOrder << std::endl);
2540 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2541 }
2542 }
2543
2544 if (DummyType == INT_MS2) {
2545 SAFENEW(pRhoDummy, NullDriveCaller);
2546 SAFENEW(pRhoAlgebraicDummy, NullDriveCaller);
2547 }
2548 break;
2549
2550 case NOSTRO:
2551 case MS:
2552 case HOPE:
2553 pRhoDummy = HP.GetDriveCaller(true);
2554
2555 if (HP.IsArg()) {
2556 pRhoAlgebraicDummy = HP.GetDriveCaller(true);
2557 } else {
2558 pRhoAlgebraicDummy = pRhoDummy->pCopy();
2559 }
2560
2561 switch (KMethod) {
2562 case NOSTRO:
2563 case MS:
2564 DummyType = INT_MS2;
2565 break;
2566
2567 case HOPE:
2568 DummyType = INT_HOPE;
2569 break;
2570
2571 default:
2572 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2573 }
2574 break;
2575
2576 case THIRDORDER:
2577 if (HP.IsKeyWord("ad" "hoc")) {
2578 /* do nothing */ ;
2579 } else {
2580 pRhoDummy = HP.GetDriveCaller(true);
2581 }
2582 DummyType = INT_THIRDORDER;
2583 break;
2584
2585 case IMPLICITEULER:
2586 DummyType = INT_IMPLICITEULER;
2587 break;
2588
2589 default:
2590 silent_cerr("Unknown integration method at line " << HP.GetLineData() << std::endl);
2591 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2592 }
2593 break;
2594 }
2595
2596 case TOLERANCE: {
2597 /*
2598 * residual tolerance; can be the keyword "null",
2599 * which means that the convergence test
2600 * will be computed on the solution, or a number
2601 */
2602 if (HP.IsKeyWord("null")) {
2603 dTol = 0.;
2604
2605 } else {
2606 dTol = HP.GetReal();
2607 if (dTol < 0.) {
2608 dTol = ::dDefaultTol;
2609 silent_cerr("warning, residual tolerance "
2610 "< 0. is illegal; "
2611 "using default value " << dTol
2612 << std::endl);
2613 }
2614 }
2615
2616 /* safe default */
2617 if (dTol == 0.) {
2618 ResTest = NonlinearSolverTest::NONE;
2619 }
2620
2621 if (HP.IsArg()) {
2622 if (HP.IsKeyWord("test")) {
2623 if (HP.IsKeyWord("norm")) {
2624 ResTest = NonlinearSolverTest::NORM;
2625 } else if (HP.IsKeyWord("minmax")) {
2626 ResTest = NonlinearSolverTest::MINMAX;
2627 } else if (HP.IsKeyWord("none")) {
2628 ResTest = NonlinearSolverTest::NONE;
2629 } else {
2630 silent_cerr("unknown test "
2631 "method at line "
2632 << HP.GetLineData()
2633 << std::endl);
2634 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2635 }
2636
2637 if (HP.IsKeyWord("scale")) {
2638 if (ResTest == NonlinearSolverTest::NONE) {
2639 silent_cerr("it's a nonsense "
2640 "to scale a disabled test; "
2641 "\"scale\" ignored"
2642 << std::endl);
2643 bScale = false;
2644 } else {
2645 bScale = true;
2646 }
2647 }
2648 }
2649 }
2650
2651 if (HP.IsArg()) {
2652 if (!HP.IsKeyWord("null")) {
2653 dSolutionTol = HP.GetReal();
2654 }
2655
2656 /* safe default */
2657 if (dSolutionTol != 0.) {
2658 SolTest = NonlinearSolverTest::NORM;
2659 }
2660
2661 if (HP.IsArg()) {
2662 if (HP.IsKeyWord("test")) {
2663 if (HP.IsKeyWord("norm")) {
2664 SolTest = NonlinearSolverTest::NORM;
2665 } else if (HP.IsKeyWord("minmax")) {
2666 SolTest = NonlinearSolverTest::MINMAX;
2667 } else if (HP.IsKeyWord("none")) {
2668 SolTest = NonlinearSolverTest::NONE;
2669 } else {
2670 silent_cerr("unknown test "
2671 "method at line "
2672 << HP.GetLineData()
2673 << std::endl);
2674 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2675 }
2676 }
2677 }
2678
2679 } else if (dTol == 0.) {
2680 silent_cerr("need solution tolerance "
2681 "with null residual tolerance"
2682 << std::endl);
2683 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2684 }
2685
2686 if (dSolutionTol < 0.) {
2687 dSolutionTol = 0.;
2688 silent_cerr("warning, solution tolerance "
2689 "< 0. is illegal; "
2690 "solution test is disabled"
2691 << std::endl);
2692 }
2693
2694 if (dTol == 0. && dSolutionTol == 0.) {
2695 silent_cerr("both residual and solution "
2696 "tolerances are zero" << std::endl);
2697 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2698 }
2699
2700 DEBUGLCOUT(MYDEBUG_INPUT, "tolerance = " << dTol
2701 << ", " << dSolutionTol << std::endl);
2702 break;
2703 }
2704
2705 case MAXRESIDUAL: {
2706 if (HP.IsKeyWord("differential" "equations")) {
2707 dMaxResidualDiff = HP.GetReal();
2708
2709 if (dMaxResidualDiff <= 0.) {
2710 silent_cerr("error: max residual for differential equations "
2711 "must be greater than zero at line "
2712 << HP.GetLineData() << std::endl);
2713 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2714 }
2715 }
2716
2717 if (HP.IsKeyWord("all" "equations") || HP.IsArg()) {
2718 dMaxResidual = HP.GetReal();
2719
2720 if (dMaxResidual <= 0.) {
2721 silent_cerr("error: max residual must be greater than zero at line "
2722 << HP.GetLineData() << std::endl);
2723 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2724 }
2725 }
2726 break;
2727 }
2728
2729 case DERIVATIVESTOLERANCE: {
2730 dDerivativesTol = HP.GetReal();
2731 if (dDerivativesTol <= 0.) {
2732 dDerivativesTol = ::dDefaultTol;
2733 silent_cerr("warning, derivatives "
2734 "tolerance <= 0.0 is illegal; "
2735 "using default value "
2736 << dDerivativesTol
2737 << std::endl);
2738 }
2739 DEBUGLCOUT(MYDEBUG_INPUT,
2740 "Derivatives tolerance = "
2741 << dDerivativesTol
2742 << std::endl);
2743 break;
2744 }
2745
2746 case MAXITERATIONS: {
2747 iMaxIterations = HP.GetInt();
2748 if (iMaxIterations < 1) {
2749 iMaxIterations = ::iDefaultMaxIterations;
2750 silent_cerr("warning, max iterations "
2751 "< 1 is illegal; using default value "
2752 << iMaxIterations
2753 << std::endl);
2754 }
2755 DEBUGLCOUT(MYDEBUG_INPUT,
2756 "Max iterations = "
2757 << iMaxIterations << std::endl);
2758 if (HP.IsKeyWord("at" "most")) {
2759 iMaxIterations = -iMaxIterations;
2760 }
2761 break;
2762 }
2763
2764 case MODIFY_RES_TEST:
2765 if (bParallel) {
2766 silent_cerr("\"modify residual test\" "
2767 "not supported by schur data manager "
2768 "at line " << HP.GetLineData()
2769 << "; ignored" << std::endl);
2770 } else {
2771 bModResTest = true;
2772 DEBUGLCOUT(MYDEBUG_INPUT,
2773 "Modify residual test" << std::endl);
2774 }
2775 break;
2776
2777 case ENFORCE_CONSTRAINT_EQUATIONS:
2778 if (HP.IsKeyWord("constraint" "violations")) {
2779 eScaleFlags = SCALE_ALGEBRAIC_EQUATIONS_YES;
2780
2781 bSetScaleAlgebraic = !HP.IsKeyWord("scale" "factor");
2782
2783 if (!bSetScaleAlgebraic) {
2784 dScaleAlgebraic = HP.GetReal();
2785 }
2786 } else if (HP.IsKeyWord("constraint" "violation" "rates")) {
2787 eScaleFlags = SCALE_ALGEBRAIC_EQUATIONS_NO;
2788 } else {
2789 silent_cerr("Keyword \"constraint violations\" or "
2790 "\"constraint violation rates\" expected at line "
2791 << HP.GetLineData() << std::endl);
2792
2793 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2794 }
2795 break;
2796
2797 case DERIVATIVESMAXITERATIONS: {
2798 iDerivativesMaxIterations = HP.GetInt();
2799 if (iDerivativesMaxIterations < 1) {
2800 iDerivativesMaxIterations = ::iDefaultMaxIterations;
2801 silent_cerr("warning, derivatives "
2802 "max iterations < 1 is illegal; "
2803 "using default value "
2804 << iDerivativesMaxIterations
2805 << std::endl);
2806 }
2807 DEBUGLCOUT(MYDEBUG_INPUT, "Derivatives "
2808 "max iterations = "
2809 << iDerivativesMaxIterations
2810 << std::endl);
2811 break;
2812 }
2813
2814 case FICTITIOUSSTEPSMAXITERATIONS:
2815 case DUMMYSTEPSMAXITERATIONS: {
2816 iDummyStepsMaxIterations = HP.GetInt();
2817 if (iDummyStepsMaxIterations < 1) {
2818 iDummyStepsMaxIterations = ::iDefaultMaxIterations;
2819 silent_cerr("warning, dummy steps "
2820 "max iterations < 1 is illegal; "
2821 "using default value "
2822 << iDummyStepsMaxIterations
2823 << std::endl);
2824 }
2825 DEBUGLCOUT(MYDEBUG_INPUT, "Dummy steps "
2826 "max iterations = "
2827 << iDummyStepsMaxIterations
2828 << std::endl);
2829 break;
2830 }
2831
2832 case DERIVATIVESCOEFFICIENT: {
2833 const bool bAutoDerivCoef = HP.IsKeyWord("auto");
2834
2835 if (!bAutoDerivCoef) {
2836 dDerivativesCoef = HP.GetReal();
2837 if (dDerivativesCoef <= 0.) {
2838 dDerivativesCoef = ::dDefaultDerivativesCoefficient;
2839 silent_cerr("warning, derivatives "
2840 "coefficient <= 0. is illegal; "
2841 "using default value "
2842 << dDerivativesCoef
2843 << std::endl);
2844 }
2845 DEBUGLCOUT(MYDEBUG_INPUT, "Derivatives coefficient = "
2846 << dDerivativesCoef << std::endl);
2847 }
2848
2849 if (bAutoDerivCoef || HP.IsKeyWord("auto")) {
2850 if (HP.IsKeyWord("max" "iterations")) {
2851 iDerivativesCoefMaxIter = HP.GetInt();
2852 if (iDerivativesCoefMaxIter < 0) {
2853 silent_cerr("number of iterations must be greater than or equal to zero at line " << HP.GetLineData() << std::endl);
2854 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2855 }
2856
2857 } else {
2858 iDerivativesCoefMaxIter = 6;
2859 }
2860
2861 if (HP.IsKeyWord("factor")) {
2862 dDerivativesCoefFactor = HP.GetReal();
2863
2864 if (dDerivativesCoefFactor <= 1.) {
2865 silent_cerr("factor must be greater than one at line " << HP.GetLineData() << std::endl);
2866 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2867 }
2868 }
2869 }
2870
2871 break;
2872 }
2873
2874 case NEWTONRAPHSON: {
2875 pedantic_cout("Newton Raphson is deprecated; use "
2876 "\"nonlinear solver: newton raphson "
2877 "[ , modified, <steps> ]\" instead"
2878 << std::endl);
2879 KeyWords NewRaph(KeyWords(HP.GetWord()));
2880 switch(NewRaph) {
2881 case MODIFIED:
2882 bTrueNewtonRaphson = 0;
2883 if (HP.IsArg()) {
2884 iIterationsBeforeAssembly = HP.GetInt();
2885 } else {
2886 iIterationsBeforeAssembly = ::iDefaultIterationsBeforeAssembly;
2887 }
2888 DEBUGLCOUT(MYDEBUG_INPUT, "Modified "
2889 "Newton-Raphson will be used; "
2890 "matrix will be assembled "
2891 "at most after "
2892 << iIterationsBeforeAssembly
2893 << " iterations" << std::endl);
2894 break;
2895
2896 default:
2897 silent_cerr("warning: unknown case; "
2898 "using default" << std::endl);
2899
2900 /* no break: fall-thru to next case */
2901 case NR_TRUE:
2902 bTrueNewtonRaphson = 1;
2903 iIterationsBeforeAssembly = 0;
2904 break;
2905 }
2906 break;
2907 }
2908
2909 case END:
2910 switch (KeyWords(HP.GetWord())) {
2911 case MULTISTEP:
2912 pedantic_cout("\"end: multistep;\" is deprecated; "
2913 "use \"end: initial value;\" instead." << std::endl);
2914 case INITIAL_VALUE:
2915 break;
2916
2917 default:
2918 silent_cerr("\"end: initial value;\" expected "
2919 "at line " << HP.GetLineData()
2920 << "; aborting..." << std::endl);
2921 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2922 }
2923 goto EndOfCycle;
2924
2925 case STRATEGY:
2926 pTSC = ReadTimeStepData(this, HP);
2927 break;
2928
2929 case POD:
2930 silent_cerr("line " << HP.GetLineData()
2931 << ": POD analysis not supported (ignored)"
2932 << std::endl);
2933 for (; HP.IsArg();) {
2934 (void)HP.GetReal();
2935 }
2936 break;
2937
2938 case EIGENANALYSIS:
2939 // initialize output precision: 0 means use default precision
2940 EigAn.iMatrixPrecision = 0;
2941 EigAn.iResultsPrecision = 0;
2942
2943 // read eigenanalysis time (to be changed)
2944 if (HP.IsKeyWord("list")) {
2945 int iNumTimes = HP.GetInt();
2946 if (iNumTimes <= 0) {
2947 silent_cerr("invalid number of eigenanalysis times "
2948 "at line " << HP.GetLineData()
2949 << std::endl);
2950 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2951 }
2952
2953 EigAn.Analyses.resize(iNumTimes);
2954 for (std::vector<doublereal>::iterator i = EigAn.Analyses.begin();
2955 i != EigAn.Analyses.end(); ++i)
2956 {
2957 *i = HP.GetReal();
2958 if (i > EigAn.Analyses.begin() && *i <= *(i-1)) {
2959 silent_cerr("eigenanalysis times must be in strict ascending order "
2960 "at line " << HP.GetLineData()
2961 << std::endl);
2962 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2963 }
2964 }
2965
2966 } else {
2967 EigAn.Analyses.resize(1);
2968 EigAn.Analyses[0] = HP.GetReal();
2969 }
2970
2971 ASSERT(EigAn.Analyses.size() > 0);
2972 // initialize EigAn
2973 EigAn.currAnalysis = EigAn.Analyses.begin();
2974 EigAn.bAnalysis = true;
2975
2976 // permute is the default; use "balance, no" to disable
2977 EigAn.uFlags = EigenAnalysis::EIG_PERMUTE;
2978
2979 while (HP.IsArg()) {
2980 if (HP.IsKeyWord("parameter")) {
2981 EigAn.dParam = HP.GetReal();
2982
2983 } else if (HP.IsKeyWord("output" "matrices")) {
2984 EigAn.uFlags |= EigenAnalysis::EIG_OUTPUT_MATRICES;
2985
2986 } else if (HP.IsKeyWord("output" "full" "matrices")) {
2987 #ifndef USE_EIG
2988 silent_cerr("\"output full matrices\" needs eigenanalysis support" << std::endl);
2989 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2990 #endif // !USE_EIG
2991 EigAn.uFlags |= EigenAnalysis::EIG_OUTPUT_FULL_MATRICES;
2992
2993 } else if (HP.IsKeyWord("output" "sparse" "matrices")) {
2994 EigAn.uFlags |= EigenAnalysis::EIG_OUTPUT_SPARSE_MATRICES;
2995
2996 } else if (HP.IsKeyWord("output" "eigenvectors")) {
2997 #ifndef USE_EIG
2998 silent_cerr("\"output eigenvectors\" needs eigenanalysis support" << std::endl);
2999 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3000 #endif // !USE_EIG
3001 EigAn.uFlags |= EigenAnalysis::EIG_OUTPUT_EIGENVECTORS;
3002
3003 } else if (HP.IsKeyWord("output" "geometry")) {
3004 EigAn.uFlags |= EigenAnalysis::EIG_OUTPUT_GEOMETRY;
3005
3006 } else if (HP.IsKeyWord("matrix" "output" "precision")) {
3007 EigAn.iMatrixPrecision = HP.GetInt();
3008 if (EigAn.iMatrixPrecision <= 0) {
3009 silent_cerr("negative or null \"matrix output precision\" "
3010 "parameter at line " << HP.GetLineData()
3011 << std::endl);
3012 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3013 }
3014
3015 } else if (HP.IsKeyWord("results" "output" "precision")) {
3016 EigAn.iResultsPrecision = HP.GetInt();
3017 if (EigAn.iResultsPrecision <= 0) {
3018 silent_cerr("negative or null \"results output precision\" "
3019 "parameter at line " << HP.GetLineData()
3020 << std::endl);
3021 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3022 }
3023
3024 } else if (HP.IsKeyWord("upper" "frequency" "limit")) {
3025 EigAn.dUpperFreq = HP.GetReal();
3026 if (EigAn.dUpperFreq < 0.) {
3027 silent_cerr("invalid \"upper frequency limit\" "
3028 "at line " << HP.GetLineData()
3029 << std::endl);
3030 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3031 }
3032
3033 } else if (HP.IsKeyWord("lower" "frequency" "limit")) {
3034 EigAn.dLowerFreq = HP.GetReal();
3035 if (EigAn.dLowerFreq < 0.) {
3036 silent_cerr("invalid \"lower frequency limit\" "
3037 "at line " << HP.GetLineData()
3038 << std::endl);
3039 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3040 }
3041
3042 } else if (HP.IsKeyWord("use" "lapack")) {
3043 if (EigAn.uFlags & EigenAnalysis::EIG_USE_MASK) {
3044 silent_cerr("eigenanalysis routine already selected "
3045 "at line " << HP.GetLineData()
3046 << std::endl);
3047 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3048 }
3049 #ifdef USE_LAPACK
3050 EigAn.uFlags |= EigenAnalysis::EIG_USE_LAPACK;
3051 #else // !USE_LAPACK
3052 silent_cerr("\"use lapack\" "
3053 "needs to configure --with-lapack "
3054 "at line " << HP.GetLineData()
3055 << std::endl);
3056 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3057 #endif // !USE_LAPACK
3058
3059 } else if (HP.IsKeyWord("use" "arpack")) {
3060 if (EigAn.uFlags & EigenAnalysis::EIG_USE_MASK) {
3061 silent_cerr("eigenanalysis routine already selected "
3062 "at line " << HP.GetLineData()
3063 << std::endl);
3064 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3065 }
3066 #ifdef USE_ARPACK
3067 EigAn.uFlags |= EigenAnalysis::EIG_USE_ARPACK;
3068
3069 EigAn.arpack.iNEV = HP.GetInt();
3070 if (EigAn.arpack.iNEV <= 0) {
3071 silent_cerr("invalid number of eigenvalues "
3072 "at line " << HP.GetLineData()
3073 << std::endl);
3074 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3075 }
3076
3077 EigAn.arpack.iNCV = HP.GetInt();
3078 if (EigAn.arpack.iNCV <= 0
3079 || EigAn.arpack.iNCV <= EigAn.arpack.iNEV + 2)
3080 {
3081 silent_cerr("invalid number of Arnoldi vectors "
3082 "(must be > NEV+2) "
3083 "at line " << HP.GetLineData()
3084 << std::endl);
3085 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3086 }
3087
3088 if (EigAn.arpack.iNCV <= 2*EigAn.arpack.iNEV) {
3089 silent_cerr("warning, possibly incorrect number of Arnoldi vectors "
3090 "(should be > 2*NEV) "
3091 "at line " << HP.GetLineData()
3092 << std::endl);
3093 }
3094
3095 EigAn.arpack.dTOL = HP.GetReal();
3096 if (EigAn.arpack.dTOL < 0.) {
3097 silent_cerr("tolerance must be non-negative "
3098 "at line " << HP.GetLineData()
3099 << std::endl);
3100 EigAn.arpack.dTOL = 0.;
3101 }
3102 #else // !USE_ARPACK
3103 silent_cerr("\"use arpack\" "
3104 "needs to configure --with-arpack "
3105 "at line " << HP.GetLineData()
3106 << std::endl);
3107 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3108 #endif // !USE_ARPACK
3109
3110 } else if (HP.IsKeyWord("use" "jdqz")) {
3111 if (EigAn.uFlags & EigenAnalysis::EIG_USE_MASK) {
3112 silent_cerr("eigenanalysis routine already selected "
3113 "at line " << HP.GetLineData()
3114 << std::endl);
3115 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3116 }
3117 #ifdef USE_JDQZ
3118 EigAn.uFlags |= EigenAnalysis::EIG_USE_JDQZ;
3119
3120 EigAn.jdqz.kmax = HP.GetInt();
3121 if (EigAn.jdqz.kmax <= 0) {
3122 silent_cerr("invalid number of eigenvalues "
3123 "at line " << HP.GetLineData()
3124 << std::endl);
3125 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3126 }
3127
3128 EigAn.jdqz.jmax = HP.GetInt();
3129 if (EigAn.jdqz.jmax < 20
3130 || EigAn.jdqz.jmax < 2*EigAn.jdqz.kmax)
3131 {
3132 silent_cerr("invalid size of the search space "
3133 "(must be >= 20 && >= 2*kmax) "
3134 "at line " << HP.GetLineData()
3135 << std::endl);
3136 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3137 }
3138
3139 EigAn.jdqz.jmin = 2*EigAn.jdqz.kmax;
3140
3141 EigAn.jdqz.eps = HP.GetReal();
3142 if (EigAn.jdqz.eps <= 0.) {
3143 silent_cerr("tolerance must be non-negative "
3144 "at line " << HP.GetLineData()
3145 << std::endl);
3146 EigAn.jdqz.eps = std::numeric_limits<doublereal>::epsilon();
3147 }
3148 #else // !USE_JDQZ
3149 silent_cerr("\"use jdqz\" "
3150 "needs to configure --with-jdqz "
3151 "at line " << HP.GetLineData()
3152 << std::endl);
3153 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3154 #endif // !USE_JDQZ
3155
3156 } else if (HP.IsKeyWord("balance")) {
3157 if (HP.IsKeyWord("no")) {
3158 EigAn.uFlags &= ~EigenAnalysis::EIG_BALANCE;
3159
3160 } else if (HP.IsKeyWord("permute")) {
3161 EigAn.uFlags |= EigenAnalysis::EIG_PERMUTE;
3162
3163 } else if (HP.IsKeyWord("scale")) {
3164 EigAn.uFlags |= EigenAnalysis::EIG_SCALE;
3165
3166 } else if (HP.IsKeyWord("all")) {
3167 EigAn.uFlags |= EigenAnalysis::EIG_BALANCE;
3168
3169 } else {
3170 silent_cerr("unknown balance option "
3171 "at line " << HP.GetLineData()
3172 << std::endl);
3173 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3174 }
3175
3176 } else if (HP.IsKeyWord("suffix" "width")) {
3177 if (HP.IsKeyWord("compute")) {
3178 EigAn.iFNameWidth = EigenAnalysis::EIGAN_WIDTH_COMPUTE;
3179
3180 } else {
3181 EigAn.iFNameWidth = HP.GetInt(0, HighParser::range_ge<integer>(0));
3182 }
3183
3184 } else if (HP.IsKeyWord("suffix" "format")) {
3185 EigAn.iFNameFormat = HP.GetStringWithDelims();
3186 // check?
3187
3188 } else {
3189 silent_cerr("unknown option "
3190 "at line " << HP.GetLineData()
3191 << std::endl);
3192 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3193 }
3194 }
3195
3196 // lower must be less than upper
3197 if (EigAn.dLowerFreq > EigAn.dUpperFreq) {
3198 silent_cerr("upper frequency limit " << EigAn.dUpperFreq
3199 << " less than lower frequency limit " << EigAn.dLowerFreq
3200 << std::endl);
3201 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3202 }
3203
3204 // if only upper is defined, make lower equal to -upper
3205 if (EigAn.dLowerFreq == -1.) {
3206 EigAn.dLowerFreq = -EigAn.dUpperFreq;
3207 }
3208
3209 switch (EigAn.uFlags & EigenAnalysis::EIG_USE_MASK) {
3210 case EigenAnalysis::EIG_USE_LAPACK:
3211 if (EigAn.uFlags & EigenAnalysis::EIG_OUTPUT_SPARSE_MATRICES) {
3212 silent_cerr("sparse matrices output "
3213 "incompatible with lapack "
3214 "at line " << HP.GetLineData()
3215 << std::endl);
3216 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3217 }
3218 if (EigAn.uFlags & EigenAnalysis::EIG_OUTPUT_MATRICES) {
3219 EigAn.uFlags |= EigenAnalysis::EIG_OUTPUT_FULL_MATRICES;
3220 }
3221 break;
3222
3223 case EigenAnalysis::EIG_USE_ARPACK:
3224 if (EigAn.uFlags & EigenAnalysis::EIG_OUTPUT_FULL_MATRICES) {
3225 silent_cerr("full matrices output "
3226 "incompatible with arpack "
3227 "at line " << HP.GetLineData()
3228 << std::endl);
3229 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3230 }
3231 if (EigAn.uFlags & EigenAnalysis::EIG_OUTPUT_MATRICES) {
3232 EigAn.uFlags |= EigenAnalysis::EIG_OUTPUT_SPARSE_MATRICES;
3233 }
3234 break;
3235
3236 case EigenAnalysis::EIG_USE_JDQZ:
3237 if (EigAn.uFlags & EigenAnalysis::EIG_OUTPUT_FULL_MATRICES) {
3238 silent_cerr("full matrices output "
3239 "incompatible with jdqz "
3240 "at line " << HP.GetLineData()
3241 << std::endl);
3242 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3243 }
3244 if (EigAn.uFlags & EigenAnalysis::EIG_OUTPUT_MATRICES) {
3245 EigAn.uFlags |= EigenAnalysis::EIG_OUTPUT_SPARSE_MATRICES;
3246 }
3247 break;
3248
3249 default:
3250 break;
3251 }
3252
3253 #ifdef USE_EIG
3254 // if an eigenanalysis routine is selected
3255 // or sparse matrix output is not requested,
3256 // force direct eigensolution
3257 if ((EigAn.uFlags & EigenAnalysis::EIG_USE_MASK)
3258 || !(EigAn.uFlags & EigenAnalysis::EIG_OUTPUT_SPARSE_MATRICES))
3259 {
3260 EigAn.uFlags |= EigenAnalysis::EIG_SOLVE;
3261 }
3262
3263 // if no eigenanalysis routine is selected,
3264 // force the use of LAPACK's
3265 if ((EigAn.uFlags & EigenAnalysis::EIG_SOLVE)
3266 && !(EigAn.uFlags & EigenAnalysis::EIG_USE_MASK))
3267 {
3268 EigAn.uFlags |= EigenAnalysis::EIG_USE_LAPACK;
3269 if (EigAn.uFlags & EigenAnalysis::EIG_OUTPUT_MATRICES) {
3270 EigAn.uFlags |= EigenAnalysis::EIG_OUTPUT_FULL_MATRICES;
3271 }
3272 }
3273 #else // !USE_EIG
3274 if (EigAn.uFlags & EigenAnalysis::EIG_OUTPUT_MATRICES) {
3275 EigAn.uFlags |= EigenAnalysis::EIG_OUTPUT_SPARSE_MATRICES;
3276 }
3277 silent_cerr("warning: \"eigenanalysis\" not supported; ignored" << std::endl);
3278 #endif // !USE_EIG
3279 break;
3280
3281 case SOLVER:
3282 silent_cerr("\"solver\" keyword at line "
3283 << HP.GetLineData()
3284 << " is deprecated; "
3285 "use \"linear solver\" instead"
3286 << std::endl);
3287 case LINEARSOLVER:
3288 ReadLinSol(CurrLinearSolver, HP);
3289 break;
3290
3291 case INTERFACESOLVER:
3292 silent_cerr("\"interface solver\" keyword at line "
3293 << HP.GetLineData()
3294 << " is deprecated; "
3295 "use \"interface linear solver\" "
3296 "instead" << std::endl);
3297 case INTERFACELINEARSOLVER:
3298 ReadLinSol(CurrIntSolver, HP, true);
3299
3300 #ifndef USE_MPI
3301 silent_cerr("Interface solver only allowed "
3302 "when compiled with MPI support" << std::endl);
3303 #endif /* ! USE_MPI */
3304 break;
3305
3306 case NONLINEARSOLVER:
3307 switch (KeyWords(HP.GetWord())) {
3308 case DEFAULT:
3309 NonlinearSolverType = NonlinearSolver::DEFAULT;
3310 break;
3311
3312 case NEWTONRAPHSON:
3313 NonlinearSolverType = NonlinearSolver::NEWTONRAPHSON;
3314 break;
3315
3316 case LINESEARCH:
3317 NonlinearSolverType = NonlinearSolver::LINESEARCH;
3318 break;
3319
3320 case MATRIXFREE:
3321 NonlinearSolverType = NonlinearSolver::MATRIXFREE;
3322 break;
3323
3324 default:
3325 silent_cerr("unknown nonlinear solver "
3326 "at line " << HP.GetLineData()
3327 << std::endl);
3328 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3329 }
3330
3331 switch (NonlinearSolverType) {
3332 case NonlinearSolver::NEWTONRAPHSON:
3333 case NonlinearSolver::LINESEARCH:
3334 bTrueNewtonRaphson = true;
3335 bKeepJac = false;
3336 iIterationsBeforeAssembly = 0;
3337
3338 if (NonlinearSolverType == NonlinearSolver::NEWTONRAPHSON && HP.IsKeyWord("true")) {
3339 break;
3340 }
3341
3342 if (HP.IsKeyWord("modified")) {
3343 bTrueNewtonRaphson = false;
3344 iIterationsBeforeAssembly = HP.GetInt();
3345
3346 if (HP.IsKeyWord("keep" "jacobian")) {
3347 pedantic_cout("Use of deprecated \"keep jacobian\" "
3348 "at line " << HP.GetLineData() << std::endl);
3349 bKeepJac = true;
3350
3351 } else if (HP.IsKeyWord("keep" "jacobian" "matrix")) {
3352 bKeepJac = true;
3353 }
3354
3355 DEBUGLCOUT(MYDEBUG_INPUT, "modified "
3356 "Newton-Raphson "
3357 "will be used; "
3358 "matrix will be "
3359 "assembled at most "
3360 "after "
3361 << iIterationsBeforeAssembly
3362 << " iterations"
3363 << std::endl);
3364 if (HP.IsKeyWord("honor" "element" "requests")) {
3365 bHonorJacRequest = true;
3366 DEBUGLCOUT(MYDEBUG_INPUT,
3367 "honor elements' "
3368 "request to update "
3369 "the preconditioner"
3370 << std::endl);
3371 }
3372 }
3373
3374 if (NonlinearSolver::LINESEARCH == NonlinearSolverType) {
3375 while (HP.IsArg()) {
3376 if (HP.IsKeyWord("tolerance" "x")) {
3377 LineSearch.dTolX = HP.GetReal();
3378 if (LineSearch.dTolX < 0.) {
3379 silent_cerr("tolerance x must be greater than or equal to zero at line " << HP.GetLineData() << std::endl);
3380 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3381 }
3382 } else if (HP.IsKeyWord("tolerance" "min")) {
3383 LineSearch.dTolMin = HP.GetReal();
3384 if (LineSearch.dTolMin < 0.) {
3385 silent_cerr("tolerance min must be greater than or equal to zero at line " << HP.GetLineData() << std::endl);
3386 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3387 }
3388 } else if (HP.IsKeyWord("max" "iterations")) {
3389 LineSearch.iMaxIterations = HP.GetInt();
3390 if (LineSearch.iMaxIterations < 0) {
3391 silent_cerr("max iterations must be greater than or equal to zero at line " << HP.GetLineData() << std::endl);
3392 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3393 }
3394 } else if (HP.IsKeyWord("alpha")) {
3395 LineSearch.dAlpha = HP.GetReal();
3396 if (LineSearch.dAlpha < 0.) {
3397 silent_cerr("alpha must be greater than or equal to zero at line " << HP.GetLineData() << std::endl);
3398 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3399 }
3400 } else if (HP.IsKeyWord("lambda" "min")) {
3401 LineSearch.dLambdaMin = HP.GetReal();
3402 if (LineSearch.dLambdaMin < 0.) {
3403 silent_cerr("lambda min must be greater than or equal to zero at line " << HP.GetLineData() << std::endl);
3404 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3405 }
3406 if (HP.IsKeyWord("relative")) {
3407 if (HP.GetYesNoOrBool()) {
3408 LineSearch.uFlags |= LineSearchParameters::RELATIVE_LAMBDA_MIN;
3409 } else {
3410 LineSearch.uFlags &= ~LineSearchParameters::RELATIVE_LAMBDA_MIN;
3411 }
3412 }
3413 } else if (HP.IsKeyWord("lambda" "factor" "min")) {
3414 LineSearch.dLambdaFactMin = HP.GetReal();
3415 if (LineSearch.dLambdaFactMin <= 0. || LineSearch.dLambdaFactMin >= 1.) {
3416 silent_cerr("lambda factor min must be in between zero and one at line" << HP.GetLineData() << std::endl);
3417 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3418 }
3419 } else if (HP.IsKeyWord("max" "step")) {
3420 LineSearch.dMaxStep = HP.GetReal();
3421 if (LineSearch.dMaxStep <= 0.) {
3422 silent_cerr("max step must be greater than zero at line " << HP.GetLineData() << std::endl);
3423 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3424 }
3425 } else if (HP.IsKeyWord("zero" "gradient")) {
3426 if (HP.IsKeyWord("continue")) {
3427 if (HP.GetYesNoOrBool()) {
3428 LineSearch.uFlags |= LineSearchParameters::ZERO_GRADIENT_CONTINUE;
3429 } else {
3430 LineSearch.uFlags &= ~LineSearchParameters::ZERO_GRADIENT_CONTINUE;
3431 }
3432 }
3433 else {
3434 silent_cerr("Keyword \"continue\" expected at line " << HP.GetLineData() << std::endl);
3435 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3436 }
3437 } else if (HP.IsKeyWord("divergence" "check")) {
3438 if (HP.GetYesNoOrBool()) {
3439 LineSearch.uFlags |= LineSearchParameters::DIVERGENCE_CHECK;
3440 } else {
3441 LineSearch.uFlags &= ~LineSearchParameters::DIVERGENCE_CHECK;
3442 }
3443 if (HP.IsKeyWord("factor")) {
3444 LineSearch.dDivergenceCheck = HP.GetReal();
3445 if (LineSearch.dDivergenceCheck <= 0.) {
3446 silent_cerr("divergence check factor must be greater than zero at line " << HP.GetLineData() << std::endl);
3447 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3448 }
3449 }
3450 } else if (HP.IsKeyWord("algorithm")) {
3451 LineSearch.uFlags &= ~LineSearchParameters::ALGORITHM;
3452 if (HP.IsKeyWord("cubic")) {
3453 LineSearch.uFlags |= LineSearchParameters::ALGORITHM_CUBIC;
3454 } else if (HP.IsKeyWord("factor")) {
3455 LineSearch.uFlags |= LineSearchParameters::ALGORITHM_FACTOR;
3456 } else {
3457 silent_cerr("Keyword \"cubic\" or \"factor\" expected at line " << HP.GetLineData() << std::endl);
3458 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3459 }
3460 } else if (HP.IsKeyWord("scale" "newton" "step")) {
3461 if (HP.GetYesNoOrBool()) {
3462 LineSearch.uFlags |= LineSearchParameters::SCALE_NEWTON_STEP;
3463 } else {
3464 LineSearch.uFlags &= ~LineSearchParameters::SCALE_NEWTON_STEP;
3465 }
3466 if (HP.IsKeyWord("min" "scale")) {
3467 LineSearch.dMinStepScale = HP.GetReal();
3468 if (LineSearch.dMinStepScale < 0. || LineSearch.dMinStepScale > 1.) {
3469 silent_cerr("min scale must be in range [0-1] at line " << HP.GetLineData() << std::endl);
3470 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3471 }
3472 }
3473 } else if (HP.IsKeyWord("print" "convergence" "info")) {
3474 if (HP.GetYesNoOrBool()) {
3475 LineSearch.uFlags |= LineSearchParameters::PRINT_CONVERGENCE_INFO;
3476 } else {
3477 LineSearch.uFlags &= ~LineSearchParameters::PRINT_CONVERGENCE_INFO;
3478 }
3479 } else if (HP.IsKeyWord("verbose")) {
3480 if (HP.GetYesNoOrBool()) {
3481 LineSearch.uFlags |= LineSearchParameters::VERBOSE_MODE;
3482 } else {
3483 LineSearch.uFlags &= ~LineSearchParameters::VERBOSE_MODE;
3484 }
3485 } else if (HP.IsKeyWord("abort" "at" "lambda" "min")) {
3486 if (HP.GetYesNoOrBool()) {
3487 LineSearch.uFlags |= LineSearchParameters::ABORT_AT_LAMBDA_MIN;
3488 } else {
3489 LineSearch.uFlags &= ~LineSearchParameters::ABORT_AT_LAMBDA_MIN;
3490 }
3491 } else if (HP.IsKeyWord("non" "negative" "slope" "continue")) {
3492 if (HP.GetYesNoOrBool()) {
3493 LineSearch.uFlags |= LineSearchParameters::NON_NEGATIVE_SLOPE_CONTINUE;
3494 } else {
3495 LineSearch.uFlags &= ~LineSearchParameters::NON_NEGATIVE_SLOPE_CONTINUE;
3496 }
3497 } else {
3498 silent_cerr("Keyword \"tolerance x\", "
3499 "\"tolerance min\", \"max iterations\", \"alpha\", "
3500 "\"lambda min\" \"lambda factor min\", \"max step\" "
3501 "\"divergence check\", \"algorithm\", \"scale newton step\" "
3502 "\"print convergence info\", \"verbose\", "
3503 "\"abort at lambda min\" "
3504 "or \"zero gradient\" expected at line " << HP.GetLineData() << std::endl);
3505 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3506 }
3507 }
3508 }
3509 break;
3510
3511 case NonlinearSolver::MATRIXFREE:
3512 switch (KeyWords(HP.GetWord())) {
3513 case DEFAULT:
3514 MFSolverType = MatrixFreeSolver::DEFAULT;
3515 break;
3516
3517
3518 case BICGSTAB:
3519 MFSolverType = MatrixFreeSolver::BICGSTAB;
3520 break;
3521
3522 case GMRES:
3523 MFSolverType = MatrixFreeSolver::GMRES;
3524 break;
3525
3526 default:
3527 silent_cerr("unknown iterative "
3528 "solver at line "
3529 << HP.GetLineData()
3530 << std::endl);
3531 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3532 }
3533
3534 if (HP.IsKeyWord("tolerance")) {
3535 dIterTol = HP.GetReal();
3536 DEBUGLCOUT(MYDEBUG_INPUT,"inner "
3537 "iterative solver "
3538 "tolerance: "
3539 << dIterTol
3540 << std::endl);
3541 }
3542
3543 if (HP.IsKeyWord("steps")) {
3544 iIterativeMaxSteps = HP.GetInt();
3545 DEBUGLCOUT(MYDEBUG_INPUT, "maximum "
3546 "number of inner "
3547 "steps for iterative "
3548 "solver: "
3549 << iIterativeMaxSteps
3550 << std::endl);
3551 }
3552
3553 if (HP.IsKeyWord("tau")) {
3554 dIterertiveTau = HP.GetReal();
3555 DEBUGLCOUT(MYDEBUG_INPUT,
3556 "tau scaling "
3557 "coefficient "
3558 "for iterative "
3559 "solver: "
3560 << dIterertiveTau
3561 << std::endl);
3562 }
3563
3564 if (HP.IsKeyWord("eta")) {
3565 dIterertiveEtaMax = HP.GetReal();
3566 DEBUGLCOUT(MYDEBUG_INPUT, "maximum "
3567 "eta coefficient "
3568 "for iterative "
3569 "solver: "
3570 << dIterertiveEtaMax
3571 << std::endl);
3572 }
3573
3574 if (HP.IsKeyWord("preconditioner")) {
3575 KeyWords KPrecond = KeyWords(HP.GetWord());
3576 switch (KPrecond) {
3577 case FULLJACOBIAN:
3578 case FULLJACOBIANMATRIX:
3579 PcType = Preconditioner::FULLJACOBIANMATRIX;
3580 if (HP.IsKeyWord("steps")) {
3581 iPrecondSteps = HP.GetInt();
3582 DEBUGLCOUT(MYDEBUG_INPUT,
3583 "number of steps "
3584 "before recomputing "
3585 "the preconditioner: "
3586 << iPrecondSteps
3587 << std::endl);
3588 }
3589 if (HP.IsKeyWord("honor" "element" "requests")) {
3590 bHonorJacRequest = true;
3591 DEBUGLCOUT(MYDEBUG_INPUT,
3592 "honor elements' "
3593 "request to update "
3594 "the preconditioner"
3595 << std::endl);
3596 }
3597 break;
3598
3599 /* add other preconditioners
3600 * here */
3601
3602 default:
3603 silent_cerr("unknown "
3604 "preconditioner "
3605 "at line "
3606 << HP.GetLineData()
3607 << std::endl);
3608 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3609 }
3610 break;
3611 }
3612 break;
3613
3614 default:
3615 ASSERT(0);
3616 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3617 }
3618 break;
3619
3620 case REALTIME:
3621 pRTSolver = ReadRTSolver(this, HP);
3622 break;
3623
3624 case THREADS:
3625 if (HP.IsKeyWord("auto")) {
3626 #ifdef USE_MULTITHREAD
3627 int n = get_nprocs();
3628 /* sanity checks ... */
3629 if (n <= 0) {
3630 silent_cerr("got " << n << " CPUs "
3631 "at line "
3632 << HP.GetLineData()
3633 << std::endl);
3634 nThreads = 1;
3635 } else {
3636 nThreads = n;
3637 }
3638 #else /* ! USE_MULTITHREAD */
3639 silent_cerr("configure with "
3640 "--enable-multithread "
3641 "for multithreaded assembly"
3642 << std::endl);
3643 #endif /* ! USE_MULTITHREAD */
3644
3645 } else if (HP.IsKeyWord("disable")) {
3646 #ifdef USE_MULTITHREAD
3647 nThreads = 1;
3648 #endif /* USE_MULTITHREAD */
3649
3650 } else {
3651 #ifdef USE_MULTITHREAD
3652 bool bAssembly = false;
3653 bool bSolver = false;
3654 bool bAll = true;
3655 unsigned nt;
3656 #endif // USE_MULTITHREAD
3657
3658 if (HP.IsKeyWord("assembly")) {
3659 #ifdef USE_MULTITHREAD
3660 bAll = false;
3661 bAssembly = true;
3662 #endif // USE_MULTITHREAD
3663
3664 } else if (HP.IsKeyWord("solver")) {
3665 #ifdef USE_MULTITHREAD
3666 bAll = false;
3667 bSolver = true;
3668 #endif // USE_MULTITHREAD
3669 }
3670
3671 #ifdef USE_MULTITHREAD
3672 nt =
3673 #endif // USE_MULTITHREAD
3674 HP.GetInt();
3675
3676 #ifdef USE_MULTITHREAD
3677 if (bAll || bAssembly) {
3678 nThreads = nt;
3679 }
3680
3681 if (bAll || bSolver) {
3682 bSolverThreads = true;
3683 nSolverThreads = nt;
3684 }
3685 #else /* ! USE_MULTITHREAD */
3686 silent_cerr("configure with "
3687 "--enable-multithread "
3688 "for multithreaded assembly"
3689 << std::endl);
3690 #endif /* ! USE_MULTITHREAD */
3691 }
3692 break;
3693
3694
3695 default:
3696 silent_cerr("unknown description at line "
3697 << HP.GetLineData() << "; aborting..."
3698 << std::endl);
3699 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3700 }
3701 }
3702
3703 EndOfCycle: /* esce dal ciclo di lettura */
3704
3705 if (pTSC == 0) {
3706 pedantic_cout("No time step control strategy defined; defaulting to NoChange time step control" );
3707 pTSC = new NoChange(this);
3708 }
3709
3710 if (dFinalTime < dInitialTime) {
3711 eAbortAfter = AFTER_ASSEMBLY;
3712 }
3713
3714 if (dFinalTime == dInitialTime) {
3715 eAbortAfter = AFTER_DERIVATIVES;
3716 }
3717
3718 /* Metodo di integrazione di default */
3719 if (!bMethod) {
3720 ASSERT(RegularType == INT_UNKNOWN);
3721
3722 /* FIXME: maybe we should use a better value
3723 * like 0.6; however, BDF should be conservative */
3724 SAFENEW(pRhoRegular, NullDriveCaller);
3725
3726 /* DriveCaller per Rho asintotico per variabili algebriche */
3727 pRhoAlgebraicRegular = pRhoRegular->pCopy();
3728
3729 RegularType = INT_MS2;
3730 }
3731
3732 /* Metodo di integrazione di default */
3733 if (iDummyStepsNumber && !bDummyStepsMethod) {
3734 ASSERT(DummyType == INT_UNKNOWN);
3735
3736 SAFENEW(pRhoDummy, NullDriveCaller);
3737
3738 /* DriveCaller per Rho asintotico per variabili algebriche */
3739 pRhoAlgebraicDummy = pRhoDummy->pCopy();
3740
3741 DummyType = INT_MS2;
3742 }
3743
3744 /* costruzione dello step solver derivative */
3745 SAFENEWWITHCONSTRUCTOR(pDerivativeSteps,
3746 DerivativeSolver,
3747 DerivativeSolver(dDerivativesTol,
3748 0.,
3749 dInitialTimeStep*dDerivativesCoef,
3750 iDerivativesMaxIterations,
3751 bModResTest,
3752 iDerivativesCoefMaxIter,
3753 dDerivativesCoefFactor));
3754
3755 /* First step prediction must always be Crank-Nicolson for accuracy */
3756 if (iDummyStepsNumber) {
3757 SAFENEWWITHCONSTRUCTOR(pFirstDummyStep,
3758 CrankNicolsonIntegrator,
3759 CrankNicolsonIntegrator(dDummyStepsTolerance,
3760 0.,
3761 iDummyStepsMaxIterations,
3762 bModResTest));
3763
3764 /* costruzione dello step solver dummy */
3765 switch (DummyType) {
3766 case INT_CRANKNICOLSON:
3767 SAFENEWWITHCONSTRUCTOR(pDummySteps,
3768 CrankNicolsonIntegrator,
3769 CrankNicolsonIntegrator(dDummyStepsTolerance,
3770 0.,
3771 iDummyStepsMaxIterations,
3772 bModResTest));
3773 break;
3774
3775 case INT_MS2:
3776 SAFENEWWITHCONSTRUCTOR(pDummySteps,
3777 MultistepSolver,
3778 MultistepSolver(dDummyStepsTolerance,
3779 0.,
3780 iDummyStepsMaxIterations,
3781 pRhoDummy,
3782 pRhoAlgebraicDummy,
3783 bModResTest));
3784 break;
3785
3786 case INT_HOPE:
3787 SAFENEWWITHCONSTRUCTOR(pDummySteps,
3788 HopeSolver,
3789 HopeSolver(dDummyStepsTolerance,
3790 dSolutionTol,
3791 iDummyStepsMaxIterations,
3792 pRhoDummy,
3793 pRhoAlgebraicDummy,
3794 bModResTest));
3795 break;
3796
3797 case INT_THIRDORDER:
3798 if (pRhoDummy == 0) {
3799 SAFENEWWITHCONSTRUCTOR(pDummySteps,
3800 AdHocThirdOrderIntegrator,
3801 AdHocThirdOrderIntegrator(dDummyStepsTolerance,
3802 dSolutionTol,
3803 iDummyStepsMaxIterations,
3804 bModResTest));
3805 } else {
3806 SAFENEWWITHCONSTRUCTOR(pDummySteps,
3807 TunableThirdOrderIntegrator,
3808 TunableThirdOrderIntegrator(dDummyStepsTolerance,
3809 dSolutionTol,
3810 iDummyStepsMaxIterations,
3811 pRhoDummy,
3812 bModResTest));
3813 }
3814 break;
3815
3816 case INT_IMPLICITEULER:
3817 SAFENEWWITHCONSTRUCTOR(pDummySteps,
3818 ImplicitEulerIntegrator,
3819 ImplicitEulerIntegrator(dDummyStepsTolerance,
3820 dSolutionTol, iDummyStepsMaxIterations,
3821 bModResTest));
3822 break;
3823
3824 default:
3825 silent_cerr("unknown dummy steps integration method"
3826 << std::endl);
3827 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3828 break;
3829 }
3830 }
3831
3832 SAFENEWWITHCONSTRUCTOR(pFirstRegularStep,
3833 CrankNicolsonIntegrator,
3834 CrankNicolsonIntegrator(dTol,
3835 dSolutionTol,
3836 iMaxIterations,
3837 bModResTest));
3838
3839 /* costruzione dello step solver per i passi normali */
3840 switch (RegularType) {
3841 case INT_CRANKNICOLSON:
3842 SAFENEWWITHCONSTRUCTOR(pRegularSteps,
3843 CrankNicolsonIntegrator,
3844 CrankNicolsonIntegrator(dTol,
3845 dSolutionTol,
3846 iMaxIterations,
3847 bModResTest));
3848 break;
3849
3850 case INT_MS2:
3851 SAFENEWWITHCONSTRUCTOR(pRegularSteps,
3852 MultistepSolver,
3853 MultistepSolver(dTol,
3854 dSolutionTol,
3855 iMaxIterations,
3856 pRhoRegular,
3857 pRhoAlgebraicRegular,
3858 bModResTest));
3859 break;
3860
3861 case INT_HOPE:
3862 SAFENEWWITHCONSTRUCTOR(pRegularSteps,
3863 HopeSolver,
3864 HopeSolver(dTol,
3865 dSolutionTol,
3866 iMaxIterations,
3867 pRhoRegular,
3868 pRhoAlgebraicRegular,
3869 bModResTest));
3870 break;
3871
3872 case INT_THIRDORDER:
3873 if (pRhoRegular == 0) {
3874 SAFENEWWITHCONSTRUCTOR(pRegularSteps,
3875 AdHocThirdOrderIntegrator,
3876 AdHocThirdOrderIntegrator(dTol,
3877 dSolutionTol,
3878 iMaxIterations,
3879 bModResTest));
3880 } else {
3881 SAFENEWWITHCONSTRUCTOR(pRegularSteps,
3882 TunableThirdOrderIntegrator,
3883 TunableThirdOrderIntegrator(dTol,
3884 dSolutionTol,
3885 iMaxIterations,
3886 pRhoRegular,
3887 bModResTest));
3888 }
3889 break;
3890
3891 case INT_IMPLICITEULER:
3892 SAFENEWWITHCONSTRUCTOR(pRegularSteps,
3893 ImplicitEulerIntegrator,
3894 ImplicitEulerIntegrator(dTol,
3895 dSolutionTol,
3896 iMaxIterations,
3897 bModResTest));
3898 break;
3899
3900 default:
3901 silent_cerr("Unknown integration method" << std::endl);
3902 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3903 break;
3904 }
3905
3906 if (bSetScaleAlgebraic) {
3907 dScaleAlgebraic = 1. / dInitialTimeStep;
3908 }
3909
3910 #ifdef USE_MULTITHREAD
3911 if (bSolverThreads) {
3912 if (CurrLinearSolver.SetNumThreads(nSolverThreads)) {
3913 silent_cerr("linear solver "
3914 << CurrLinearSolver.GetSolverName()
3915 << " does not support "
3916 "threaded solution" << std::endl);
3917 }
3918 }
3919 #endif /* USE_MULTITHREAD */
3920 }
3921
3922 static int
do_eig(const doublereal & b,const doublereal & re,const doublereal & im,const doublereal & h,doublereal & sigma,doublereal & omega,doublereal & csi,doublereal & freq)3923 do_eig(const doublereal& b, const doublereal& re,
3924 const doublereal& im, const doublereal& h,
3925 doublereal& sigma, doublereal& omega,
3926 doublereal& csi, doublereal& freq)
3927 {
3928 // denominator
3929 doublereal d = re + b;
3930 d *= d;
3931 d += im*im;
3932 d *= h/2.;
3933
3934 // real & imag
3935 sigma = (re*re - b*b + im*im)/d;
3936 omega = 2.*b*im/d;
3937
3938 // frequency and damping factor
3939 if (im != 0.) {
3940 d = sigma*sigma + omega*omega;
3941 if (d > std::numeric_limits<doublereal>::epsilon()) {
3942 csi = -100*sigma/sqrt(d);
3943
3944 } else {
3945 csi = 0.;
3946 }
3947
3948 freq = omega/(2*M_PI);
3949
3950 } else {
3951 if (std::abs(sigma) < std::numeric_limits<doublereal>::epsilon()) {
3952 csi = 0.;
3953
3954 } else {
3955 csi = -100.*copysign(1, sigma);
3956 }
3957
3958 freq = 0.;
3959 }
3960
3961 return 0;
3962 }
3963
3964 // Writes eigenvalues to the .out file in human-readable form
3965 static void
output_eigenvalues(const VectorHandler * pBeta,const VectorHandler & R,const VectorHandler & I,const doublereal & dShiftR,DataManager * pDM,const Solver::EigenAnalysis * pEA,integer iLow,integer iHigh,std::vector<bool> & vOut)3966 output_eigenvalues(const VectorHandler *pBeta,
3967 const VectorHandler& R, const VectorHandler& I,
3968 const doublereal& dShiftR,
3969 DataManager* pDM,
3970 const Solver::EigenAnalysis *pEA,
3971 integer iLow, integer iHigh,
3972 std::vector<bool>& vOut)
3973 {
3974 std::ostream& Out = pDM->GetOutFile();
3975
3976 /* Output? */
3977 Out << "Mode n. " " " " Real " " " " Imag " " " " " " Damp % " " Freq Hz" << std::endl;
3978
3979 integer iNVec = R.iGetSize();
3980
3981 for (int iCnt = 1; iCnt <= iNVec; iCnt++) {
3982 doublereal b = pBeta ? (*pBeta)(iCnt) : 1.;
3983 doublereal re = R(iCnt) + dShiftR;
3984 doublereal im = I(iCnt);
3985 doublereal sigma;
3986 doublereal omega;
3987 doublereal csi;
3988 doublereal freq;
3989
3990 if (iCnt < iLow || iCnt > iHigh) {
3991 vOut[iCnt - 1] = false;
3992 continue;
3993 }
3994
3995 const doublereal& h = pEA->dParam;
3996 do_eig(b, re, im, h, sigma, omega, csi, freq);
3997
3998 if (freq < pEA->dLowerFreq || freq > pEA->dUpperFreq) {
3999 vOut[iCnt - 1] = false;
4000 continue;
4001 }
4002
4003 vOut[iCnt - 1] = true;
4004
4005 Out << std::setw(8) << iCnt << ": "
4006 << std::setw(12) << sigma << " + " << std::setw(12) << omega << " j";
4007
4008 if (fabs(csi) > std::numeric_limits<doublereal>::epsilon()) {
4009 Out << " " << std::setw(12) << csi;
4010 } else {
4011 Out << " " << std::setw(12) << 0.;
4012 }
4013
4014 Out << " " << std::setw(12) << freq;
4015
4016 Out << std::endl;
4017 }
4018 }
4019
4020 #ifdef USE_LAPACK
4021 // Computes eigenvalues and eigenvectors using LAPACK's
4022 // generalized non-symmetric eigenanalysis
4023 static void
eig_lapack(const MatrixHandler * pMatA,const MatrixHandler * pMatB,DataManager * pDM,Solver::EigenAnalysis * pEA,bool bNewLine,const unsigned uCurr)4024 eig_lapack(const MatrixHandler* pMatA, const MatrixHandler* pMatB,
4025 DataManager *pDM, Solver::EigenAnalysis *pEA,
4026 bool bNewLine, const unsigned uCurr)
4027 {
4028 const FullMatrixHandler& MatA = dynamic_cast<const FullMatrixHandler &>(*pMatA);
4029 const FullMatrixHandler& MatB = dynamic_cast<const FullMatrixHandler &>(*pMatB);
4030
4031 char sB[2] = "N";
4032 if ((pEA->uFlags & Solver::EigenAnalysis::EIG_BALANCE)
4033 == Solver::EigenAnalysis::EIG_BALANCE)
4034 {
4035 sB[0] = 'B';
4036
4037 } else if (pEA->uFlags & Solver::EigenAnalysis::EIG_PERMUTE) {
4038 sB[0] = 'P';
4039
4040 } else if (pEA->uFlags & Solver::EigenAnalysis::EIG_SCALE) {
4041 sB[0] = 'S';
4042 }
4043
4044 char sL[2] = "V";
4045 char sR[2] = "V";
4046 char sS[2] = "N";
4047
4048 // iNumDof is a member, set after dataman constr.
4049 integer iSize = MatA.iGetNumRows();
4050
4051 // Minimum workspace size. To be improved.
4052 // NOTE: optimal iWorkSize is computed by dggev() and dggevx()
4053 // when called with iWorkSize = -1.
4054 // The computed value is stored in WorkVec[0].
4055 integer iWorkSize = -1;
4056 integer iMinWorkSize = -1;
4057 integer iInfo = 0;
4058
4059 integer iILO = 1, iIHI = iSize;
4060 doublereal dABNRM = -1., dBBNRM = -1.;
4061
4062 doublereal dDmy;
4063 integer iDmy;
4064 logical lDmy;
4065 doublereal dWV;
4066 if (sB[0] == 'N') {
4067 iMinWorkSize = 8*iSize;
4068
4069 __FC_DECL__(dggev)(
4070 sL, // JOBVL
4071 sR, // JOBVR
4072 &iSize, // N
4073 &dDmy, // A
4074 &iSize, // LDA
4075 &dDmy, // B
4076 &iSize, // LDB
4077 &dDmy, // ALPHAR
4078 &dDmy, // ALPHAI
4079 &dDmy, // BETA
4080 &dDmy, // VL
4081 &iSize, // LDVL
4082 &dDmy, // VR
4083 &iSize, // LDVR
4084 &dWV, // WORK
4085 &iWorkSize, // LWORK
4086 &iInfo);
4087
4088 } else {
4089 iMinWorkSize = 6*iSize;
4090
4091 __FC_DECL__(dggevx)(
4092 sB, // BALANC
4093 sL, // JOBVL
4094 sR, // JOBVR
4095 sS, // SENSE
4096 &iSize, // N
4097 &dDmy, // A
4098 &iSize, // LDA
4099 &dDmy, // B
4100 &iSize, // LDB
4101 &dDmy, // ALPHAR
4102 &dDmy, // ALPHAI
4103 &dDmy, // BETA
4104 &dDmy, // VL
4105 &iSize, // LDVL
4106 &dDmy, // VR
4107 &iSize, // LDVR
4108 &iILO, // ILO
4109 &iIHI, // IHI
4110 &dDmy, // LSCALE
4111 &dDmy, // RSCALE
4112 &dABNRM, // ABNRM
4113 &dBBNRM, // BBNRM
4114 &dDmy, // RCONDE
4115 &dDmy, // RCONDV
4116 &dWV, // WORK
4117 &iWorkSize, // LWORK
4118 &iDmy, // IWORK
4119 &lDmy, // BWORK
4120 &iInfo);
4121 }
4122
4123 if (iInfo != 0) {
4124 if (bNewLine && silent_err) {
4125 silent_cerr(std::endl);
4126 bNewLine = false;
4127 }
4128 silent_cerr("dggev[x]() query for worksize failed "
4129 "INFO=" << iInfo << std::endl);
4130 iInfo = 0;
4131 }
4132
4133 iWorkSize = (integer)dWV;
4134 if (iWorkSize < iMinWorkSize) {
4135 if (bNewLine && silent_err) {
4136 silent_cerr(std::endl);
4137 bNewLine = false;
4138 }
4139 silent_cerr("dggev[x]() asked for a worksize " << iWorkSize
4140 << " less than the minimum, " << iMinWorkSize
4141 << "; using the minimum" << std::endl);
4142 iWorkSize = iMinWorkSize;
4143 }
4144
4145 // Workspaces
4146 // 2 matrices iSize x iSize
4147 // 5 vectors iSize x 1
4148 // 1 vector iWorkSize x 1
4149 doublereal* pd = 0;
4150 int iTmpSize = 2*(iSize*iSize) + 3*iSize + iWorkSize;
4151 if (sB[0] != 'N') {
4152 iTmpSize += 2*iSize;
4153 }
4154 SAFENEWARR(pd, doublereal, iTmpSize);
4155 #if defined HAVE_MEMSET
4156 memset(pd, 0, iTmpSize*sizeof(doublereal));
4157 #else // !HAVE_MEMSET
4158 for (int iCnt = iTmpSize; iCnt-- > 0; ) {
4159 pd[iCnt] = 0.;
4160 }
4161 #endif // !HAVE_MEMSET
4162
4163 // 2 pointer arrays iSize x 1 for the matrices
4164 doublereal** ppd = 0;
4165 SAFENEWARR(ppd, doublereal*, 2*iSize);
4166
4167 // Data Handlers
4168 doublereal* pdTmp = pd;
4169 doublereal** ppdTmp = ppd;
4170
4171 FullMatrixHandler MatVL(pdTmp, ppdTmp, iSize*iSize, iSize, iSize);
4172 pdTmp += iSize*iSize;
4173 ppdTmp += iSize;
4174
4175 FullMatrixHandler MatVR(pdTmp, ppdTmp, iSize*iSize, iSize, iSize);
4176 pdTmp += iSize*iSize;
4177
4178 MyVectorHandler AlphaR(iSize, pdTmp);
4179 pdTmp += iSize;
4180
4181 MyVectorHandler AlphaI(iSize, pdTmp);
4182 pdTmp += iSize;
4183
4184 MyVectorHandler Beta(iSize, pdTmp);
4185 pdTmp += iSize;
4186
4187 MyVectorHandler LScale;
4188 MyVectorHandler RScale;
4189 if (sB[0] != 'N') {
4190 LScale.Attach(iSize, pdTmp);
4191 pdTmp += iSize;
4192
4193 RScale.Attach(iSize, pdTmp);
4194 pdTmp += iSize;
4195 }
4196
4197 MyVectorHandler WorkVec(iWorkSize, pdTmp);
4198
4199 // Eigenanalysis
4200 // NOTE: according to lapack's documentation, dgegv() is deprecated
4201 // in favour of dggev()... I find dgegv() a little bit faster (10%?)
4202 // for typical problems (N ~ 1000).
4203 if (sB[0] == 'N') {
4204 __FC_DECL__(dggev)(
4205 sL, // JOBVL
4206 sR, // JOBVR
4207 &iSize, // N
4208 const_cast<doublereal *>(MatA.pdGetMat()), // A; remove const'ness (hack)
4209 &iSize, // LDA
4210 const_cast<doublereal *>(MatB.pdGetMat()), // B; remove const'ness (hack)
4211 &iSize, // LDB
4212 AlphaR.pdGetVec(), // ALPHAR
4213 AlphaI.pdGetVec(), // ALPHAI
4214 Beta.pdGetVec(), // BETA
4215 MatVL.pdGetMat(), // VL
4216 &iSize, // LDVL
4217 MatVR.pdGetMat(), // VR
4218 &iSize, // LDVR
4219 WorkVec.pdGetVec(), // WORK
4220 &iWorkSize, // LWORK
4221 &iInfo); // INFO
4222
4223 } else {
4224 std::vector<integer> iIWORK(iSize + 6);
4225
4226 __FC_DECL__(dggevx)(
4227 sB, // BALANCE
4228 sL, // JOBVL
4229 sR, // JOBVR
4230 sS, // SENSE
4231 &iSize, // N
4232 const_cast<doublereal *>(MatA.pdGetMat()), // A; remove const'ness (hack)
4233 &iSize, // LDA
4234 const_cast<doublereal *>(MatB.pdGetMat()), // B; remove const'ness (hack)
4235 &iSize, // LDB
4236 AlphaR.pdGetVec(), // ALPHAR
4237 AlphaI.pdGetVec(), // ALPHAI
4238 Beta.pdGetVec(), // BETA
4239 MatVL.pdGetMat(), // VL
4240 &iSize, // LDVL
4241 MatVR.pdGetMat(), // VR
4242 &iSize, // LDVR
4243 &iILO, // ILO
4244 &iIHI, // IHI
4245 LScale.pdGetVec(), // LSCALE
4246 RScale.pdGetVec(), // RSCALE
4247 &dABNRM, // ABNRM
4248 &dBBNRM, // BBNRM
4249 &dDmy, // RCONDE
4250 &dDmy, // RCONDV
4251 WorkVec.pdGetVec(), // WORK
4252 &iWorkSize, // LWORK
4253 &iIWORK[0], // IWORK
4254 &lDmy, // BWORK
4255 &iInfo); // INFO
4256 }
4257
4258 std::ostream& Out = pDM->GetOutFile();
4259 Out << "Info: " << iInfo << ", ";
4260
4261 if (iInfo == 0) {
4262 // = 0: successful exit
4263 Out << "success"
4264 << " BALANC=\"" << sB << "\""
4265 << " ILO=" << iILO
4266 << " IHI=" << iIHI
4267 << " ABNRM=" << dABNRM
4268 << " BBNRM=" << dBBNRM
4269 << std::endl;
4270
4271 } else if (iInfo < 0) {
4272 const char *arg[] = {
4273 "JOBVL",
4274 "JOBVR",
4275 "N",
4276 "A",
4277 "LDA",
4278 "B",
4279 "LDB",
4280 "ALPHAR",
4281 "ALPHAI",
4282 "BETA",
4283 "VL",
4284 "LDVL",
4285 "VR",
4286 "LDVR",
4287 "WORK",
4288 "LWORK",
4289 "INFO",
4290 NULL
4291 };
4292
4293 const char *argx[] = {
4294 "BALANCE",
4295 "JOBVL",
4296 "JOBVR",
4297 "SENSE",
4298 "N",
4299 "A",
4300 "LDA",
4301 "B",
4302 "LDB",
4303 "ALPHAR",
4304 "ALPHAI",
4305 "BETA",
4306 "VL",
4307 "LDVL",
4308 "VR",
4309 "LDVR",
4310 "ILO",
4311 "IHI",
4312 "LSCALE",
4313 "RSCALE",
4314 "ABNRM",
4315 "BBNRM",
4316 "RCONDE",
4317 "RCONDV",
4318 "WORK",
4319 "LWORK",
4320 "IWORK",
4321 "BWORK",
4322 "INFO",
4323 NULL
4324 };
4325
4326 const char **argv = (sB[0] == 'N' ? arg : argx );
4327
4328 // < 0: if INFO = -i, the i-th argument had an illegal value.
4329 Out << "argument #" << -iInfo
4330 << " (" << argv[-iInfo - 1] << ") "
4331 << "was passed an illegal value" << std::endl;
4332
4333 } else if (iInfo > 0 && iInfo <= iSize) {
4334 /* = 1,...,N:
4335 * The QZ iteration failed. No eigenvectors have been
4336 * calculated, but ALPHAR(j), ALPHAI(j), and BETA(j)
4337 * should be correct for j=INFO+1,...,N. */
4338 Out << "the QZ iteration failed, but eigenvalues "
4339 << iInfo + 1 << "->" << iSize << " should be correct"
4340 << std::endl;
4341
4342 } else if (iInfo > iSize) {
4343 const char* const sErrs[] = {
4344 "DHGEQZ (other than QZ iteration failed in DHGEQZ)",
4345 "DTGEVC",
4346 NULL
4347 };
4348
4349 Out << "error return from " << sErrs[iInfo - iSize - 1]
4350 << std::endl;
4351 }
4352
4353 std::vector<bool> vOut(iSize);
4354 output_eigenvalues(&Beta, AlphaR, AlphaI, 0., pDM, pEA, iILO, iIHI, vOut);
4355
4356 if (pEA->uFlags & Solver::EigenAnalysis::EIG_OUTPUT_GEOMETRY) {
4357 pDM->OutputEigGeometry(uCurr, pEA->iResultsPrecision);
4358 }
4359
4360 if (pEA->uFlags & Solver::EigenAnalysis::EIG_OUTPUT_EIGENVECTORS) {
4361 pDM->OutputEigenvectors(&Beta, AlphaR, AlphaI, 0.,
4362 &MatVL, MatVR, vOut, uCurr, pEA->iResultsPrecision);
4363 }
4364
4365 SAFEDELETEARR(pd);
4366 SAFEDELETEARR(ppd);
4367 }
4368 #endif // USE_LAPACK
4369
4370 #ifdef USE_ARPACK
4371 // Computes eigenvalues and eigenvectors using ARPACK's
4372 // canonical non-symmetric eigenanalysis
4373 static void
eig_arpack(const MatrixHandler * pMatA,SolutionManager * pSM,DataManager * pDM,Solver::EigenAnalysis * pEA,bool bNewLine,const unsigned uCurr)4374 eig_arpack(const MatrixHandler* pMatA, SolutionManager* pSM,
4375 DataManager *pDM, Solver::EigenAnalysis *pEA,
4376 bool bNewLine, const unsigned uCurr)
4377 {
4378 NaiveSparsePermSolutionManager<Colamd_ordering>& sm
4379 = dynamic_cast<NaiveSparsePermSolutionManager<Colamd_ordering> &>(*pSM);
4380 const NaiveMatrixHandler& MatA = dynamic_cast<const NaiveMatrixHandler &>(*pMatA);
4381
4382 // shift
4383 doublereal SIGMAR = 0.;
4384 doublereal SIGMAI = 0.;
4385
4386 // arpack-related vars
4387 integer IDO; // 0 at first iteration; then set by dnaupd
4388 const char *BMAT; // 'I' for standard problem
4389 integer N; // size of problem
4390 const char *WHICH; // "SM" to request smallest eigenvalues
4391 integer NEV; // number of eigenvalues
4392 doublereal TOL; // -1 to use machine precision
4393 std::vector<doublereal> RESID; // residual vector (ignored if IDO==0)
4394 integer NCV; // number of vectors in subspace
4395 std::vector<doublereal> V; // Schur basis
4396 integer LDV; // leading dimension of V (==N!)
4397 integer IPARAM[11] = { 0 };
4398 integer IPNTR[14] = { 0 };
4399 std::vector<doublereal> WORKD;
4400 std::vector<doublereal> WORKL;
4401 integer LWORKL;
4402 integer INFO;
4403
4404 IDO = 0;
4405 BMAT = "I";
4406 N = MatA.iGetNumRows();
4407 WHICH = "SM";
4408 NEV = pEA->arpack.iNEV;
4409 if (NEV > N) {
4410 silent_cerr("eig_arpack: invalid NEV=" << NEV << " > size of problem (=" << N << ")" << std::endl);
4411 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
4412 }
4413
4414 TOL = pEA->arpack.dTOL;
4415 RESID.resize(N, 0.);
4416 NCV = pEA->arpack.iNCV;
4417 if (NCV > N) {
4418 silent_cerr("eig_arpack: invalid NCV=" << NCV << " > size of problem (=" << N << ")" << std::endl);
4419 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
4420 }
4421 // NOTE: we accommodate NCV + 1 vectors to handle the case
4422 // of a complex last eigenvalue which requires a vector
4423 // for the real part and one for the imaginary part
4424 V.resize(N*(NCV + 1), 0.);
4425 LDV = N;
4426 IPARAM[0] = 1;
4427 IPARAM[2] = 300; // configurable?
4428 IPARAM[3] = 1;
4429 IPARAM[6] = 1; // mode 1: canonical problem
4430 WORKD.resize(3*N, 0.);
4431 LWORKL = 3*NCV*NCV + 6*NCV;
4432 WORKL.resize(LWORKL, 0.);
4433 INFO = 0;
4434
4435 int cnt = 0;
4436
4437 const bool bOutputStatus = isatty(fileno(stderr));
4438
4439 do {
4440 __FC_DECL__(dnaupd)(&IDO, &BMAT[0], &N, &WHICH[0], &NEV,
4441 &TOL, &RESID[0], &NCV, &V[0], &LDV, &IPARAM[0], &IPNTR[0],
4442 &WORKD[0], &WORKL[0], &LWORKL, &INFO);
4443
4444 #if 0
4445 std::cout << "cnt=" << cnt << ": IDO=" << IDO << ", INFO=" << INFO << std::endl;
4446 #endif
4447
4448 // compute Y = OP*X
4449 ASSERT(IPNTR[0] - 1 <= 2*N);
4450 ASSERT(IPNTR[1] - 1 <= 2*N);
4451 MyVectorHandler X(N, &WORKD[IPNTR[0] - 1]);
4452 MyVectorHandler Y(N, &WORKD[IPNTR[1] - 1]);
4453
4454 /*
4455 * NOTE: we are solving the problem
4456
4457 MatB * X * Lambda = MatA * X
4458
4459 * and we want to focus on Ritz parameters Lambda
4460 * as close as possible to (1., 0.), which maps
4461 * to (0., 0.) in continuous time.
4462 *
4463 * We are casting the problem in the form
4464
4465 X * Alpha = A * X
4466
4467 * by putting the problem in canonical form
4468
4469 X * Lambda = MatB \ MatA * X
4470
4471 * and then subtracting a shift Sigma = (1., 0) after :
4472
4473 X * Lambda - X = MatB \ MatA * X - X
4474
4475 X * (Lambda - 1.) = (MatB \ MatA - I) * X
4476
4477 * so
4478
4479 Alpha = Lambda - 1.
4480
4481 A = MatB \ MatA - I
4482
4483 * and the sequence of operations for Y = A * X is
4484
4485 X' = MatA * X
4486 X'' = MatB \ X'
4487 Y = X'' - X
4488
4489 * the eigenvalues need to be modified by adding 1.
4490 */
4491
4492
4493 MatA.MatVecMul(*sm.pResHdl(), X);
4494 sm.Solve();
4495 *sm.pSolHdl() -= X;
4496
4497 Y = *sm.pSolHdl();
4498
4499 static const int CNT = 100;
4500 cnt++;
4501 if (bOutputStatus && !(cnt % CNT)) {
4502 if (bNewLine && silent_err) {
4503 silent_cerr(std::endl);
4504 bNewLine = false;
4505 }
4506 silent_cerr("\r" "cnt=" << cnt);
4507 }
4508
4509 if (mbdyn_stop_at_end_of_iteration()) {
4510 if (bNewLine && silent_err) {
4511 silent_cerr(std::endl);
4512 bNewLine = false;
4513 }
4514 silent_cerr((cnt >= CNT ? "\n" : "")
4515 << "ARPACK: interrupted" << std::endl);
4516 return;
4517 }
4518 } while (IDO == 1 || IDO == -1);
4519
4520 if (bOutputStatus) {
4521 if (bNewLine && silent_err) {
4522 silent_cerr(std::endl);
4523 bNewLine = false;
4524 }
4525 silent_cerr("\r" "cnt=" << cnt << std::endl);
4526 }
4527
4528 // NOTE: improve diagnostics
4529 if (INFO < 0) {
4530 if (bNewLine && silent_err) {
4531 silent_cerr(std::endl);
4532 bNewLine = false;
4533 }
4534 silent_cerr("ARPACK error after " << cnt << " iterations; "
4535 "IDO=" << IDO << ", INFO=" << INFO << std::endl);
4536 return;
4537 }
4538
4539 switch (INFO) {
4540 case 0:
4541 break;
4542
4543 case 1:
4544 if (bNewLine && silent_err) {
4545 silent_cerr(std::endl);
4546 bNewLine = false;
4547 }
4548 silent_cerr("INFO=1: Maximum number of iterations taken. "
4549 "All possible eigenvalues of OP have been found. IPARAM(5) "
4550 "returns the number of wanted converged Ritz values "
4551 "(currently = " << IPARAM[4] << "; requested NEV = " << NEV << ")."
4552 << std::endl);
4553 break;
4554
4555 case 2:
4556 if (bNewLine && silent_err) {
4557 silent_cerr(std::endl);
4558 bNewLine = false;
4559 }
4560 silent_cerr("INFO=2: No longer an informational error. Deprecated starting "
4561 "with release 2 of ARPACK."
4562 << std::endl);
4563 break;
4564
4565 case 3:
4566 if (bNewLine && silent_err) {
4567 silent_cerr(std::endl);
4568 bNewLine = false;
4569 }
4570 silent_cerr("INFO=3: No shifts could be applied during a cycle of the "
4571 "implicitly restarted Arnoldi iteration. One possibility "
4572 "is to increase the size of NCV (currently = " << NCV << ") "
4573 "relative to NEV (currently = " << NEV << "). "
4574 "See remark 4 in dnaupd(3)."
4575 << std::endl);
4576 break;
4577
4578 default:
4579 if (bNewLine && silent_err) {
4580 silent_cerr(std::endl);
4581 bNewLine = false;
4582 }
4583 silent_cerr("INFO=" << INFO << ": undocumented value." << std::endl);
4584 break;
4585 }
4586
4587 std::ostream& Out = pDM->GetOutFile();
4588 Out << "INFO: " << INFO << std::endl;
4589
4590 #if 0
4591 for (int ii = 0; ii < 11; ii++) {
4592 std::cerr << "IPARAM(" << ii + 1 << ")=" << IPARAM[ii] << std::endl;
4593 }
4594 #endif
4595
4596 logical RVEC = true;
4597 const char *HOWMNY = "A";
4598 std::vector<logical> SELECT(NCV);
4599 std::vector<doublereal> D(2*NCV);
4600 doublereal *DR = &D[0], *DI = &D[NCV];
4601 std::vector<doublereal> Z(N*(NCV + 1));
4602 integer LDZ = N;
4603 std::vector<doublereal> WORKEV(3*NCV);
4604
4605 #if 0
4606 std::cerr << "dneupd:" << std::endl
4607 << "RVEC = " << RVEC << std::endl
4608 << "HOWMNY = " << HOWMNY << std::endl
4609 << "SELECT(" << NCV << ")" << std::endl
4610 << "DR(" << NEV + 1 << ")" << std::endl
4611 << "DI(" << NEV + 1 << ")" << std::endl
4612 << "Z(" << N << ", " << NEV + 1 << ")" << std::endl
4613 << "LDZ = " << LDZ << std::endl
4614 << "SIGMAR = " << SIGMAR << std::endl
4615 << "SIGMAI = " << SIGMAI << std::endl
4616 << "WORKEV(" << 3*NCV << ")" << std::endl
4617 << "BMAT = " << BMAT << std::endl
4618 << "N = " << N << std::endl
4619 << "WHICH = " << WHICH << std::endl
4620 << "NEV = " << NEV << std::endl
4621 << "TOL = " << TOL << std::endl
4622 << "RESID(" << N << ")" << std::endl
4623 << "NCV = " << NCV << std::endl
4624 << "V(" << N << ", " << NCV << ")" << std::endl
4625 << "LDV = " << LDV << std::endl
4626 << "IPARAM(" << sizeof(IPARAM)/sizeof(IPARAM[0]) << ")" << std::endl
4627 << "IPNTR(" << sizeof(IPNTR)/sizeof(IPNTR[0]) << ")" << std::endl
4628 << "WORKD(" << 3*N << ")" << std::endl
4629 << "WORKL(" << LWORKL << ")" << std::endl
4630 << "LWORKL = " << LWORKL << std::endl
4631 << "INFO = " << INFO << std::endl;
4632 #endif
4633
4634 __FC_DECL__(dneupd)(&RVEC, &HOWMNY[0], &SELECT[0], &DR[0], &DI[0],
4635 &Z[0], &LDZ, &SIGMAR, &SIGMAI, &WORKEV[0],
4636 &BMAT[0], &N, &WHICH[0], &NEV,
4637 &TOL, &RESID[0], &NCV, &V[0], &LDV, &IPARAM[0], &IPNTR[0],
4638 &WORKD[0], &WORKL[0], &LWORKL, &INFO);
4639
4640 int nconv = IPARAM[4];
4641 if (nconv > 0) {
4642 ASSERT(nconv <= NCV);
4643 MyVectorHandler AlphaR(nconv, DR);
4644 MyVectorHandler AlphaI(nconv, DI);
4645 std::vector<bool> vOut(nconv);
4646 output_eigenvalues(0, AlphaR, AlphaI, 1., pDM, pEA, 1, nconv, vOut);
4647
4648 if (pEA->uFlags & Solver::EigenAnalysis::EIG_OUTPUT_GEOMETRY) {
4649 pDM->OutputEigGeometry(uCurr, pEA->iResultsPrecision);
4650 }
4651
4652 if (pEA->uFlags & Solver::EigenAnalysis::EIG_OUTPUT_EIGENVECTORS) {
4653 std::vector<doublereal *> ZC(nconv);
4654 FullMatrixHandler VR(&Z[0], &ZC[0], N*nconv, N, nconv);
4655 pDM->OutputEigenvectors(0, AlphaR, AlphaI, 1.,
4656 0, VR, vOut, uCurr, pEA->iResultsPrecision);
4657 }
4658
4659 } else {
4660 if (bNewLine && silent_err) {
4661 silent_cerr(std::endl);
4662 bNewLine = false;
4663 }
4664 silent_cerr("no converged Ritz coefficients" << std::endl);
4665 }
4666 }
4667 #endif // USE_ARPACK
4668
4669 #ifdef USE_JDQZ
4670 // Computes eigenvalues and eigenvectors using ARPACK's
4671 // canonical non-symmetric eigenanalysis
4672 static void
eig_jdqz(const MatrixHandler * pMatA,const MatrixHandler * pMatB,DataManager * pDM,Solver::EigenAnalysis * pEA,bool bNewLine,const unsigned uCurr)4673 eig_jdqz(const MatrixHandler *pMatA, const MatrixHandler *pMatB,
4674 DataManager *pDM, Solver::EigenAnalysis *pEA,
4675 bool bNewLine, const unsigned uCurr)
4676 {
4677 const NaiveMatrixHandler& MatA = dynamic_cast<const NaiveMatrixHandler &>(*pMatA);
4678 const NaiveMatrixHandler& MatB = dynamic_cast<const NaiveMatrixHandler &>(*pMatB);
4679
4680 MBJDQZ mbjdqz(MatA, MatB);
4681 mbjdqzp = &mbjdqz;
4682
4683 /*
4684
4685 alpha, beta Obvious from equation (1)
4686 wanted Compute the converged eigenvectors (if wanted =
4687 .true.)
4688 eivec Converged eigenvectors if wanted = .true., else con-
4689 verged Schur vectors
4690 n The size of the problem
4691 target The value near which the eigenvalues are sought
4692 eps Tolerance of the eigensolutions, Ax-Bx /|/| <
4693 kmax Number of wanted eigensolutions, on output: number of
4694 converged eigenpairs
4695 jmax Maximum size of the search space
4696 jmin Minimum size of the search space
4697 method Linear equation solver:
4698 1: GMRESm , [2]
4699 2: BiCGstab(), [3]
4700 m Maximum dimension of searchspace of GMRESm
4701 l Degree of GMRES-polynomial in Bi-CGstab()
4702 mxmv Maximum number of matrix-vector multiplications in
4703 GMRESm or BiCGstab()
4704 maxstep Maximum number of Jacobi-Davidson iterations
4705 lock Tracking parameter (section 2.5.1)
4706 order Selection criterion for Ritz values:
4707 0: nearest to target
4708 -1: smallest real part
4709 1: largest real part
4710 -2: smallest imaginary part
4711 2: largest imaginary part
4712 testspace Determines how to expand the testspace W
4713 1: w = "Standard Petrov" ×v (Section 3.1.1)
4714 2: w = "Standard 'variable' Petrov" ×v (Section 3.1.2)
4715 3: w = "Harmonic Petrov" ×v (Section 3.5.1)
4716 zwork Workspace
4717 lwork Size of workspace, >= 4+m+5jmax+3kmax if GMRESm
4718 is used, >= 10 + 6 + 5jmax + 3kmax if Bi-CGstab() is
4719 used.
4720
4721 */
4722
4723 std::vector<doublecomplex> alpha;
4724 std::vector<doublecomplex> beta;
4725 std::vector<doublecomplex> eivec;
4726 logical wanted = 1;
4727 integer n = pMatA->iGetNumRows();
4728 doublecomplex target = { 1., 0. };
4729 doublereal eps = pEA->jdqz.eps;
4730 integer kmax = pEA->jdqz.kmax;
4731 integer jmax = pEA->jdqz.jmax;
4732 integer jmin = pEA->jdqz.jmin;
4733 integer method = pEA->jdqz.method;
4734 integer m = pEA->jdqz.m;
4735 integer l = pEA->jdqz.l;
4736 integer mxmv = pEA->jdqz.mxmv;
4737 integer maxstep = pEA->jdqz.maxstep;
4738 doublereal lock = pEA->jdqz.lock;
4739 integer order = pEA->jdqz.order;
4740 integer testspace = pEA->jdqz.testspace;
4741 std::vector<doublecomplex> zwork;
4742 integer lwork;
4743
4744 switch (method) {
4745 case Solver::EigenAnalysis::JDQZ::GMRES:
4746 lwork = 4 + m + 5*jmax + 3*kmax;
4747 break;
4748
4749 case Solver::EigenAnalysis::JDQZ::BICGSTAB:
4750 lwork = 10 + 6*l + 5*jmax + 3*kmax;
4751 break;
4752
4753 default:
4754 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
4755 }
4756
4757 alpha.resize(jmax);
4758 beta.resize(jmax);
4759 eivec.resize(n*kmax);
4760 zwork.resize(n*lwork);
4761
4762 __FC_DECL__(jdqz)(
4763 &alpha[0],
4764 &beta[0],
4765 &eivec[0],
4766 &wanted,
4767 &n,
4768 &target,
4769 &eps,
4770 &kmax,
4771 &jmax,
4772 &jmin,
4773 &method,
4774 &m,
4775 &l,
4776 &mxmv,
4777 &maxstep,
4778 &lock,
4779 &order,
4780 &testspace,
4781 &zwork[0],
4782 &lwork);
4783
4784 if (bNewLine && silent_err) {
4785 silent_cerr(std::endl);
4786 bNewLine = false;
4787 }
4788 silent_cerr("\r" "cnt=" << mbjdqz.Cnt() << std::endl);
4789
4790 if (kmax > 0) {
4791 int nconv = kmax;
4792 MyVectorHandler AlphaR(nconv);
4793 MyVectorHandler AlphaI(nconv);
4794 MyVectorHandler Beta(nconv);
4795 for (integer c = 0; c < nconv; c++) {
4796 if (beta[c].i != 0.) {
4797 doublereal d = std::sqrt(beta[c].r*beta[c].r + beta[c].i*beta[c].i);
4798 Beta(c + 1) = d;
4799 AlphaR(c + 1) = (alpha[c].r*beta[c].r + alpha[c].i*beta[c].i)/d;
4800 AlphaI(c + 1) = (alpha[c].i*beta[c].r - alpha[c].r*beta[c].i)/d;
4801
4802 } else {
4803 Beta(c + 1) = beta[c].r;
4804 AlphaR(c + 1) = alpha[c].r;
4805 AlphaI(c + 1) = alpha[c].i;
4806 }
4807 }
4808 std::vector<bool> vOut(nconv);
4809 output_eigenvalues(0, AlphaR, AlphaI, 1., pDM, pEA, 1, nconv, vOut);
4810
4811 if (pEA->uFlags & Solver::EigenAnalysis::EIG_OUTPUT_GEOMETRY) {
4812 pDM->OutputGeometry(pEA->iResultsPrecision);
4813 }
4814
4815 if (pEA->uFlags & Solver::EigenAnalysis::EIG_OUTPUT_EIGENVECTORS) {
4816 FullMatrixHandler VR(n, nconv);
4817 doublecomplex *p = &eivec[0] - 1;
4818 for (integer c = 1; c <= nconv; c++) {
4819
4820 if (AlphaI(c) == 0.) {
4821 for (integer r = 1; r <= n; r++) {
4822 VR(r, c) = p[r].r;
4823 }
4824
4825 } else {
4826 for (integer r = 1; r <= n; r++) {
4827 VR(r, c) = p[r].r;
4828 VR(r, c + 1) = p[n + r].i;
4829 }
4830
4831 p += n;
4832 c++;
4833 }
4834
4835 p += n;
4836 }
4837
4838 pDM->OutputEigenvectors(&Beta, AlphaR, AlphaI, 1.,
4839 0, VR, vOut, uCurr, pEA->iResultsPrecision);
4840 }
4841
4842 } else {
4843 if (bNewLine && silent_err) {
4844 silent_cerr(std::endl);
4845 bNewLine = false;
4846 }
4847 silent_cerr("no converged eigenpairs" << std::endl);
4848 }
4849 }
4850 #endif // USE_JDQZ
4851
4852 // Driver for eigenanalysis
4853 void
Eig(bool bNewLine)4854 Solver::Eig(bool bNewLine)
4855 {
4856 DEBUGCOUTFNAME("Solver::Eig");
4857
4858 /*
4859 * MatA, MatB: MatrixHandlers to eigenanalysis matrices
4860 * MatVL, MatVR: MatrixHandlers to eigenvectors, if required
4861 * AlphaR, AlphaI Beta: eigenvalues
4862 * WorkVec: Workspace
4863 * iWorkSize: Size of the workspace
4864 */
4865
4866 DEBUGCOUT("Solver::Eig(): performing eigenanalysis" << std::endl);
4867
4868 integer iSize = iNumDofs;
4869
4870 SolutionManager *pSM = 0;
4871 MatrixHandler *pMatA = 0;
4872 MatrixHandler *pMatB = 0;
4873
4874 if (EigAn.uFlags & EigenAnalysis::EIG_USE_LAPACK) {
4875 SAFENEWWITHCONSTRUCTOR(pMatA, FullMatrixHandler,
4876 FullMatrixHandler(iSize));
4877
4878 SAFENEWWITHCONSTRUCTOR(pMatB, FullMatrixHandler,
4879 FullMatrixHandler(iSize));
4880
4881 } else if (EigAn.uFlags & EigenAnalysis::EIG_USE_ARPACK) {
4882 SAFENEWWITHCONSTRUCTOR(pSM, NaiveSparsePermSolutionManager<Colamd_ordering>,
4883 NaiveSparsePermSolutionManager<Colamd_ordering>(iSize));
4884 SAFENEWWITHCONSTRUCTOR(pMatA, NaiveMatrixHandler,
4885 NaiveMatrixHandler(iSize));
4886 pMatB = pSM->pMatHdl();
4887
4888 } else if (EigAn.uFlags & EigenAnalysis::EIG_USE_JDQZ) {
4889 SAFENEWWITHCONSTRUCTOR(pMatA, NaiveMatrixHandler,
4890 NaiveMatrixHandler(iSize));
4891 SAFENEWWITHCONSTRUCTOR(pMatB, NaiveMatrixHandler,
4892 NaiveMatrixHandler(iSize));
4893
4894 } else if (EigAn.uFlags & EigenAnalysis::EIG_OUTPUT_SPARSE_MATRICES) {
4895 SAFENEWWITHCONSTRUCTOR(pMatA, SpMapMatrixHandler,
4896 SpMapMatrixHandler(iSize));
4897 SAFENEWWITHCONSTRUCTOR(pMatB, SpMapMatrixHandler,
4898 SpMapMatrixHandler(iSize));
4899 }
4900
4901 pMatA->Reset();
4902 pMatB->Reset();
4903
4904 // Matrices assembly (see eig.ps)
4905 doublereal h = EigAn.dParam;
4906 pDM->AssJac(*pMatA, -h/2.);
4907 pDM->AssJac(*pMatB, h/2.);
4908
4909 #ifdef DEBUG
4910 DEBUGCOUT(std::endl
4911 << "Matrix A:" << std::endl << *pMatA << std::endl
4912 << "Matrix B:" << std::endl << *pMatB << std::endl);
4913 #endif /* DEBUG */
4914
4915 unsigned uCurr = EigAn.currAnalysis - EigAn.Analyses.begin();
4916 if (EigAn.uFlags & EigenAnalysis::EIG_OUTPUT) {
4917 unsigned uSize = EigAn.Analyses.size();
4918 if (uSize >= 1) {
4919 std::stringstream postfix_ss;
4920
4921 postfix_ss << '_';
4922
4923 if (EigAn.iFNameWidth > 0) {
4924 postfix_ss << std::setw(EigAn.iFNameWidth) << std::setfill('0') << uCurr;
4925
4926 } else if (!EigAn.iFNameFormat.empty()) {
4927 char buf[BUFSIZ];
4928
4929 snprintf(buf, sizeof(buf), EigAn.iFNameFormat.c_str(), (int)uCurr);
4930 postfix_ss << buf;
4931
4932 } else {
4933 postfix_ss << uCurr;
4934 }
4935
4936 pDM->OutputEigOpen(postfix_ss.str());
4937 pDM->OutputEigParams(dTime, h/2., uCurr, EigAn.iResultsPrecision);
4938 }
4939 }
4940 if (EigAn.uFlags & EigenAnalysis::EIG_OUTPUT_FULL_MATRICES) {
4941 pDM->OutputEigFullMatrices(pMatA, pMatB, uCurr, EigAn.iMatrixPrecision);
4942 }
4943 else if (EigAn.uFlags & EigenAnalysis::EIG_OUTPUT_SPARSE_MATRICES) {
4944 if (dynamic_cast<const NaiveMatrixHandler *>(pMatB)) {
4945 pDM->OutputEigNaiveMatrices(pMatA, pMatB, uCurr, EigAn.iMatrixPrecision);
4946 } else {
4947 pDM->OutputEigSparseMatrices(pMatA, pMatB, uCurr, EigAn.iMatrixPrecision);
4948 }
4949 }
4950
4951 switch (EigAn.uFlags & EigenAnalysis::EIG_USE_MASK) {
4952 #ifdef USE_LAPACK
4953 case EigenAnalysis::EIG_USE_LAPACK:
4954 eig_lapack(pMatA, pMatB, pDM, &EigAn, bNewLine, uCurr);
4955 break;
4956 #endif // USE_LAPACK
4957
4958 #ifdef USE_ARPACK
4959 case EigenAnalysis::EIG_USE_ARPACK:
4960 eig_arpack(pMatA, pSM, pDM, &EigAn, bNewLine, uCurr);
4961 break;
4962 #endif // USE_ARPACK
4963
4964 #ifdef USE_JDQZ
4965 case EigenAnalysis::EIG_USE_JDQZ:
4966 eig_jdqz(pMatA, pMatB, pDM, &EigAn, bNewLine, uCurr);
4967 break;
4968 #endif // USE_JDQZ
4969
4970 default:
4971 // only output matrices, use external eigenanalysis
4972 break;
4973 }
4974
4975 pDM->OutputEigClose();
4976
4977 if (pSM) {
4978 pMatB = 0;
4979 SAFEDELETE(pSM);
4980 }
4981
4982 if (pMatA) {
4983 SAFEDELETE(pMatA);
4984 }
4985
4986 if (pMatB) {
4987 SAFEDELETE(pMatB);
4988 }
4989 }
4990
dGetInitialMaxTimeStep() const4991 doublereal Solver::dGetInitialMaxTimeStep() const
4992 {
4993 if (typeid(*MaxTimeStep.pGetDriveCaller()) == typeid(PostponedDriveCaller)) {
4994 return ::dDefaultMaxTimeStep;
4995 }
4996
4997 // The same behavior like in previous releases
4998 return MaxTimeStep.dGet();
4999 }
5000
5001 SolutionManager *const
AllocateSolman(integer iNLD,integer iLWS)5002 Solver::AllocateSolman(integer iNLD, integer iLWS)
5003 {
5004 SolutionManager *pCurrSM = CurrLinearSolver.GetSolutionManager(iNLD, iLWS);
5005
5006 /* special extra parameters if required */
5007 switch (CurrLinearSolver.GetSolver()) {
5008 case LinSol::UMFPACK_SOLVER:
5009 #ifdef HAVE_UMFPACK_TIC_DISABLE
5010 if (pRTSolver) {
5011 /* disable profiling, to avoid times() system call
5012 *
5013 * This function has been introduced in Umfpack 4.1
5014 * by our patch at
5015 *
5016 * http://mbdyn.aero.polimi.it/~masarati/Download/\
5017 * mbdyn/umfpack-4.1-nosyscalls.patch
5018 *
5019 * but since Umfpack 4.3 is no longer required,
5020 * provided the library is compiled with -DNO_TIMER
5021 * to disable run-time syscalls to timing routines.
5022 */
5023 umfpack_tic_disable();
5024 }
5025 #endif // HAVE_UMFPACK_TIC_DISABLE
5026 break;
5027
5028 default:
5029 break;
5030 }
5031
5032 return pCurrSM;
5033 };
5034
5035
5036 SolutionManager *const
AllocateSchurSolman(integer iStates)5037 Solver::AllocateSchurSolman(integer iStates)
5038 {
5039 SolutionManager *pSSM(0);
5040
5041 #ifdef USE_MPI
5042 switch (CurrIntSolver.GetSolver()) {
5043 case LinSol::LAPACK_SOLVER:
5044 case LinSol::MESCHACH_SOLVER:
5045 case LinSol::NAIVE_SOLVER:
5046 case LinSol::UMFPACK_SOLVER:
5047 case LinSol::Y12_SOLVER:
5048 break;
5049
5050 default:
5051 silent_cerr("apparently solver "
5052 << CurrIntSolver.GetSolverName()
5053 << " is not allowed as interface solver "
5054 "for SchurSolutionManager" << std::endl);
5055 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
5056 }
5057
5058 SAFENEWWITHCONSTRUCTOR(pSSM,
5059 SchurSolutionManager,
5060 SchurSolutionManager(iNumDofs, iStates, pLocDofs,
5061 iNumLocDofs,
5062 pIntDofs, iNumIntDofs,
5063 pLocalSM, CurrIntSolver));
5064
5065 #else /* !USE_MPI */
5066 silent_cerr("Configure --with-mpi to enable Schur solver" << std::endl);
5067 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
5068 #endif /* !USE_MPI */
5069
5070 return pSSM;
5071 };
5072
5073 NonlinearSolver *const
AllocateNonlinearSolver()5074 Solver::AllocateNonlinearSolver()
5075 {
5076 NonlinearSolver *pNLS = 0;
5077
5078 switch (NonlinearSolverType) {
5079 case NonlinearSolver::MATRIXFREE:
5080 switch (MFSolverType) {
5081 case MatrixFreeSolver::BICGSTAB:
5082 SAFENEWWITHCONSTRUCTOR(pNLS,
5083 BiCGStab,
5084 BiCGStab(PcType,
5085 iPrecondSteps,
5086 dIterTol,
5087 iIterativeMaxSteps,
5088 dIterertiveEtaMax,
5089 dIterertiveTau,
5090 *this));
5091 break;
5092
5093 default:
5094 pedantic_cout("unknown matrix free solver type; "
5095 "using default" << std::endl);
5096 /* warning: should be unreachable */
5097
5098 case MatrixFreeSolver::GMRES:
5099 SAFENEWWITHCONSTRUCTOR(pNLS,
5100 Gmres,
5101 Gmres(PcType,
5102 iPrecondSteps,
5103 dIterTol,
5104 iIterativeMaxSteps,
5105 dIterertiveEtaMax,
5106 dIterertiveTau,
5107 *this));
5108 break;
5109 }
5110 break;
5111
5112 default:
5113 pedantic_cout("unknown nonlinear solver type; using default"
5114 << std::endl);
5115
5116 case NonlinearSolver::NEWTONRAPHSON:
5117 SAFENEWWITHCONSTRUCTOR(pNLS,
5118 NewtonRaphsonSolver,
5119 NewtonRaphsonSolver(bTrueNewtonRaphson,
5120 bKeepJac,
5121 iIterationsBeforeAssembly,
5122 *this));
5123 break;
5124 case NonlinearSolver::LINESEARCH:
5125 SAFENEWWITHCONSTRUCTOR(pNLS,
5126 LineSearchSolver,
5127 LineSearchSolver(pDM,
5128 *this,
5129 LineSearch));
5130 break;
5131 }
5132 return pNLS;
5133 }
5134
5135 void
SetupSolmans(integer iStates,bool bCanBeParallel)5136 Solver::SetupSolmans(integer iStates, bool bCanBeParallel)
5137 {
5138 DEBUGLCOUT(MYDEBUG_MEM, "creating SolutionManager\n\tsize = "
5139 << iNumDofs*iUnkStates <<
5140 "\n\tnumdofs = " << iNumDofs
5141 << "\n\tnumstates = " << iStates << std::endl);
5142
5143 /* delete previous solmans */
5144 if (pSM != 0) {
5145 SAFEDELETE(pSM);
5146 pSM = 0;
5147 }
5148 if (pLocalSM != 0) {
5149 SAFEDELETE(pLocalSM);
5150 pLocalSM = 0;
5151 }
5152
5153 integer iWorkSpaceSize = CurrLinearSolver.iGetWorkSpaceSize();
5154 integer iLWS = iWorkSpaceSize;
5155 integer iNLD = iNumDofs*iStates;
5156 if (bCanBeParallel && bParallel) {
5157 /* FIXME BEPPE! */
5158 iLWS = iWorkSpaceSize*iNumLocDofs/(iNumDofs*iNumDofs);
5159 /* FIXME: GIUSTO QUESTO? */
5160 iNLD = iNumLocDofs*iStates;
5161 }
5162
5163 SolutionManager *pCurrSM = AllocateSolman(iNLD, iLWS);
5164
5165 /*
5166 * This is the LOCAL solver if instantiating a parallel
5167 * integrator; otherwise it is the MAIN solver
5168 */
5169 if (bCanBeParallel && bParallel) {
5170 pLocalSM = pCurrSM;
5171
5172 /* Crea il solutore di Schur globale */
5173 pSM = AllocateSchurSolman(iStates);
5174
5175 } else {
5176 pSM = pCurrSM;
5177 }
5178 /*
5179 * FIXME: at present there MUST be a pSM
5180 * (even for matrix-free nonlinear solvers)
5181 */
5182 if (pSM == 0) {
5183 silent_cerr("No linear solver defined" << std::endl);
5184 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
5185 }
5186 }
5187
5188 clock_t
GetCPUTime(void) const5189 Solver::GetCPUTime(void) const
5190 {
5191 return pDM->GetCPUTime();
5192 }
5193
5194 void
PrintResidual(const VectorHandler & Res,integer iIterCnt) const5195 Solver::PrintResidual(const VectorHandler& Res, integer iIterCnt) const
5196 {
5197 pDM->PrintResidual(Res, iIterCnt);
5198 }
5199
5200 void
PrintSolution(const VectorHandler & Sol,integer iIterCnt) const5201 Solver::PrintSolution(const VectorHandler& Sol, integer iIterCnt) const
5202 {
5203 pDM->PrintSolution(Sol, iIterCnt);
5204 }
5205
CheckTimeStepLimit(doublereal dErr,doublereal dErrDiff) const5206 void Solver::CheckTimeStepLimit(doublereal dErr, doublereal dErrDiff) const throw(NonlinearSolver::MaxResidualExceeded, NonlinearSolver::TimeStepLimitExceeded)
5207 {
5208 if (pDerivativeSteps) {
5209 // Time step cannot be reduced
5210 return;
5211 }
5212
5213 if (dErr > dMaxResidual) {
5214 if (dCurrTimeStep > 2 * dMinTimeStep) { // FIXME
5215 if (outputIters()) {
5216 #ifdef USE_MPI
5217 if (!bParallel || MBDynComm.Get_rank() == 0)
5218 #endif /* USE_MPI */
5219 {
5220 silent_cerr("warning: current residual = " << dErr
5221 << " > maximum residual = " << dMaxResidual
5222 << std::endl);
5223 }
5224 }
5225
5226 throw NonlinearSolver::MaxResidualExceeded(MBDYN_EXCEPT_ARGS);
5227 }
5228 }
5229
5230 if (dErrDiff > dMaxResidualDiff) {
5231 if (dCurrTimeStep > 2 * dMinTimeStep) { // FIXME
5232 if (outputIters()) {
5233 #ifdef USE_MPI
5234 if (!bParallel || MBDynComm.Get_rank() == 0)
5235 #endif /* USE_MPI */
5236 {
5237 silent_cerr("warning: current residual for differential equations = " << dErrDiff
5238 << " > maximum residual = " << dMaxResidualDiff
5239 << std::endl);
5240 }
5241 }
5242
5243 throw NonlinearSolver::MaxResidualExceeded(MBDYN_EXCEPT_ARGS);
5244 }
5245 }
5246
5247 switch (eTimeStepLimit) {
5248 case TS_SOFT_LIMIT:
5249 break;
5250
5251 case TS_HARD_LIMIT: {
5252 const doublereal dMaxTS = MaxTimeStep.dGet();
5253
5254 if (dCurrTimeStep > dMaxTS && dCurrTimeStep > dMinTimeStep) {
5255 if (outputIters()) {
5256 #ifdef USE_MPI
5257 if (!bParallel || MBDynComm.Get_rank() == 0)
5258 #endif /* USE_MPI */
5259 {
5260 silent_cerr("warning: current time step = "
5261 << dCurrTimeStep
5262 << " > hard limit of the maximum time step = "
5263 << dMaxTS << std::endl);
5264 }
5265 }
5266
5267 throw NonlinearSolver::TimeStepLimitExceeded(MBDYN_EXCEPT_ARGS);
5268 }
5269 } break;
5270
5271 default:
5272 ASSERT(0);
5273 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
5274 }
5275 }
5276