MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
invsolver.cc
Go to the documentation of this file.
1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/base/invsolver.cc,v 1.58 2017/09/09 09:20:12 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) 2008
43  * Alessandro Fumagalli <alessandro.fumagalli@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 <unistd.h>
58 #include <cerrno>
59 
60 #include <cfloat>
61 #include <cmath>
62 #include "ac/sys_sysinfo.h"
63 
64 #include "invsolver.h"
65 #include "dataman.h"
66 #include "mtdataman.h"
67 #include "thirdorderstepsol.h"
68 #include "nr.h"
69 #include "bicg.h"
70 #include "gmres.h"
71 #include "solman.h"
72 #include "stepsol.h"
73 #include <vector>
74 #include "readlinsol.h"
75 #include "drive_.h"
76 
77 #include "solver_impl.h"
78 
80  const std::string& sInFName,
81  const std::string& sOutFName,
82  unsigned int nThreads,
83  bool bPar)
84 : Solver(HPar, sInFName, sOutFName, nThreads, bPar),
85 ProblemType(InverseDynamics::FULLY_ACTUATED_COLLOCATED),
86 pXPrimePrime(0), pLambda(0),
87 bFullResTest(false)
88 {
89  DEBUGCOUTFNAME("InverseSolver::InverseSolver");
90 }
91 
92 bool
94 {
95  DEBUGCOUTFNAME("InverseSolver::Prepare");
96 
97  // consistency check
99  silent_cerr("Prepare() must be called first" << std::endl);
101  }
102 
104 
105  /* Legge i dati relativi al metodo di integrazione */
106  ReadData(HP);
107 
108 /*FIXME:*/
109 // bParallel = false;
110 
111 #ifdef USE_MULTITHREAD
112  ThreadPrepare();
113 #endif /* USE_MULTITHREAD */
114 
115  if (pRTSolver) {
116  pRTSolver->Setup();
117  }
118 
119  int mpi_finalize = 0;
120 
121 #ifdef USE_SCHUR
122  if (bParallel) {
123  DEBUGLCOUT(MYDEBUG_MEM, "creating parallel SchurDataManager"
124  << std::endl);
125 
129  OutputFlags,
130  this,
131  dInitialTime,
132  sOutputFileName.c_str(),
133  sInputFileName.c_str(),
135 
136  /* FIXME: who frees sNewOutname? */
137 
138  pDM = pSDM;
139 
140  } else
141 #endif // USE_SCHUR
142  {
143  /* chiama il gestore dei dati generali della simulazione */
144 #ifdef USE_MULTITHREAD
145  if (nThreads > 1) {
147  /* conservative: dir may use too much memory */
149  bool b;
150 
151 #if defined(USE_UMFPACK)
154 #elif defined(USE_Y12)
157 #else
158  b = false;
159 #endif
160  if (!b) {
161  silent_cerr("unable to select a CC-capable solver"
162  << std::endl);
164  }
165  }
166  }
167 
168  silent_cout("Creating multithread solver "
169  "with " << nThreads << " threads "
170  "and "
172  << " linear solver"
173  << std::endl);
174 
176  MultiThreadDataManager,
177  MultiThreadDataManager(HP,
178  OutputFlags,
179  this,
180  dInitialTime,
181  sOutputFileName.c_str(),
182  sInputFileName.c_str(),
184  nThreads));
185 
186  } else
187 #endif /* USE_MULTITHREAD */
188  {
189  DEBUGLCOUT(MYDEBUG_MEM, "creating DataManager"
190  << std::endl);
191 
192  silent_cout("Creating scalar solver "
193  "with "
195  << " linear solver"
196  << std::endl);
197 
199  DataManager,
200  DataManager(HP,
201  OutputFlags,
202  this,
203  dInitialTime,
204  sOutputFileName.c_str(),
205  sInputFileName.c_str(),
207  }
208  }
209 
210  { // log of symbol table
211  std::ostream& out = pDM->GetLogFile();
212  out << HP.GetMathParser().GetSymbolTable();
213  }
214  HP.Close();
215 
216  /* Si fa dare l'std::ostream al file di output per il log */
217  std::ostream& Out = pDM->GetOutFile();
218 
219  /* Qui crea le partizioni: principale fra i processi, se parallelo */
220 #ifdef USE_SCHUR
221  if (bParallel) {
223  }
224 #endif // USE_SCHUR
225 
226  const DriveHandler* pDH = pDM->pGetDrvHdl();
228  /* Costruisce i vettori della soluzione ai vari passi */
229  DEBUGLCOUT(MYDEBUG_MEM, "creating solution vectors" << std::endl);
230 
231 #ifdef USE_SCHUR
232  if (bParallel) {
234  pDofs = pSDM->pGetDofsList();
235 
238 
241 
242  } else
243 #endif // USE_SCHUR
244  {
245  iNumDofs = pDM->iGetNumDofs();
246  }
247 
248  /* relink those known drive callers that might need
249  * the data manager, but were verated ahead of it */
250 
251  ASSERT(iNumDofs > 0);
252 
253 #if 0
254  /* sono i passi precedenti usati dall'integratore */
257 #endif // 0
258 
259  /* FIXME: pdWorkspace?*/
260 
261 #if 1
262  /* allocate workspace for previous time steps */
265  /* allocate MyVectorHandlers for previous time steps: use workspace */
266  for (int ivec = 0; ivec < iNumPreviousVectors; ivec++) {
269  MyVectorHandler(iNumDofs, pdWorkSpace+ivec*iNumDofs));
270  qX.push_back(pX);
273  MyVectorHandler(iNumDofs,
274  pdWorkSpace+((iNumPreviousVectors)+ivec)*iNumDofs));
275  qXPrime.push_back(pXPrime);
276  pX = NULL;
277  pXPrime = NULL;
278  }
279  /* allocate MyVectorHandlers for unknown time step(s): own memory */
280 #endif
281 
284  MyVectorHandler(iUnkStates*iNumDofs));
287  MyVectorHandler(iUnkStates*iNumDofs));
290  MyVectorHandler(iUnkStates*iNumDofs));
293  MyVectorHandler(iUnkStates*iNumDofs));
294 
295  pX->Reset();
296  pXPrime->Reset();
297  pXPrimePrime->Reset();
298  pLambda->Reset();
299 
300  /* relink those known drive callers that might need
301  * the data manager, but were verated ahead of it */
303 
304  /*
305  * Immediately link DataManager to current solution
306  *
307  * this should work as long as the last unknown time step is put
308  * at the beginning of pX, pXPrime
309  */
310 
312 
313  /* a questo punto si costruisce il nonlinear solver */
315 
316  /* FIXME: Serve?*/
317  MyVectorHandler Scale(iNumDofs);
318  if (bScale) {
319  /* collects scale factors from data manager */
320  pDM->SetScale(Scale);
321  }
322 
323 
324  /*
325  * prepare tests for nonlinear solver;
326  *
327  * test on residual may allow pre-scaling;
328  * test on solution (difference between two iterations) does not
329  */
330  NonlinearSolverTest *pResTest = NULL;
331  if (bScale) {
332  NonlinearSolverTestScale *pResTestScale = NULL;
333 
334  switch (ResTest) {
336  SAFENEW(pResTestScale, NonlinearSolverTestScaleNorm);
337  break;
338 
340  SAFENEW(pResTestScale, NonlinearSolverTestScaleMinMax);
341  break;
342 
343  default:
344  ASSERT(0);
346  }
347 
348  /* registers scale factors at nonlinear solver */
349  pResTestScale->SetScale(&Scale);
350 
351  pResTest = pResTestScale;
352 
353 
354  } else {
355  switch (ResTest) {
357  SAFENEW(pResTest, NonlinearSolverTestNone);
358  break;
359 
361  SAFENEW(pResTest, NonlinearSolverTestNorm);
362  break;
363 
366  break;
367 
368  default:
369  ASSERT(0);
371  }
372  }
373 
374  NonlinearSolverTest *pSolTest = NULL;
375  switch (SolTest) {
377  SAFENEW(pSolTest, NonlinearSolverTestNone);
378  break;
379 
381  SAFENEW(pSolTest, NonlinearSolverTestNorm);
382  break;
383 
386  break;
387 
388  default:
389  ASSERT(0);
391  }
392 
393  // pippero
394  switch (GetProblemType()) {
397  {
400  pDM->IDSetTest(pRT, pST, bFullResTest);
401  pResTest = pRT;
402  pSolTest = pST;
403  } break;
404 
405  default:
406  break;
407  }
408 
409  /* registers tests in nonlinear solver */
410  pNLS->SetTest(pResTest, pSolTest);
411 
412  /*
413  * Dell'assemblaggio iniziale dei vincoli se ne occupa il DataManager
414  * in quanto e' lui il responsabile dei dati della simulazione,
415  * e quindi anche della loro coerenza. Inoltre e' lui a sapere
416  * quali equazioni sono di vincolo o meno.
417  */
418 
419  pDM->SetValue(*pX, *pXPrime);
420 
421 
422  /*
423  * Prepare output
424  */
425  pDM->OutputPrepare();
426 
427  /*
428  * Dialoga con il DataManager per dargli il tempo iniziale
429  * e per farsi inizializzare i vettori di soluzione e derivata
430  */
431  /* FIXME: the time is already set by DataManager, but FileDrivers
432  * have not been ServePending'd
433  */
435  pDM->SetTime(dTime + dInitialTimeStep, dInitialTimeStep, 0);
436 
437  // dTest = 0.;
438  // dSolTest = 0.;
439 
441 
442  /* settaggio degli output Types */
443  unsigned OF = OutputFlags;
445  OF |= OUTPUT_RES;
446  }
448  OF |= OUTPUT_JAC;
449  }
451  OF |= OUTPUT_SOL;
452  }
453  pNLS->SetOutputFlags(OF);
454  if (pOutputMeter) {
457  }
458 
461 
462 #ifdef USE_EXTERNAL
463  pNLS->SetExternal(External::EMPTY);
464 #endif /* USE_EXTERNAL */
465 
466  dCurrTimeStep = 0.;
467 
468 #ifdef USE_EXTERNAL
469  /* comunica che gli ultimi dati inviati sono la condizione iniziale */
470  External::SendInitial();
471 #endif /* USE_EXTERNAL */
472 
473 #ifdef USE_EXTERNAL
474  /* il prossimo passo e' un regular */
475  pNLS->SetExternal(External::REGULAR);
476 #endif /* USE_EXTERNAL */
477 
478  lStep = -1;
479  DEBUGCOUT("Current time step: " << dCurrTimeStep << std::endl);
480 
482  /* Fa l'output della soluzione al primo passo ed esce */
483  Out << "Interrupted during first step." << std::endl;
485  }
486 
487  bSolConv = false;
490 
491  if (pRTSolver) {
492  pRTSolver->Init();
493  }
494 
495  bOutputCounter = outputCounter() && isatty(fileno(stderr));
496  outputCounterPrefix = bOutputCounter ? "\n" : "";
497  outputCounterPostfix = outputStep() ? "\n" : "\r";
498 
499  pTSC->Init(iMaxIterations, dMinTimeStep, MaxTimeStep, dInitialTimeStep);
500 
501  /* Setup SolutionManager(s) */
502  ASSERT(pRegularSteps != 0);
504 
505  /* this is here to force AfterConvergence()
506  * since we need to explicitly SetTime() to the next time step */
507  pDM->SetValue(*pX, *pXPrime);
508 
510 
512 
513  return true;
514 }
515 
516 bool
518 {
519  DEBUGCOUTFNAME("InverseSolver::Start");
520 
521  // consistency check
523  silent_cerr("Start() must be called after Prepare()" << std::endl);
525  }
526 
527  // temporarily push status forward
529  bool b = Advance();
530  if (!b) {
531  // roll status back in case of failure
533  }
534  return b;
535 }
536 
537 bool
539 {
540  DEBUGCOUTFNAME("InverseSolver::Advance");
541 
542  // consistency check
544  silent_cerr("Started() must be called first" << std::endl);
546  }
547 
549 
550  if (dTime >= dFinalTime) {
551  if (pRTSolver) {
553  }
554 
555  silent_cout(outputCounterPrefix
556  << "End of simulation at time "
557  << dTime << " after "
558  << lStep << " steps;" << std::endl
559  << "output in file \"" << sOutputFileName << "\"" << std::endl
560  << "total iterations: " << iTotIter << std::endl
561  << "total Jacobian matrices: " << pNLS->TotalAssembledJacobian() << std::endl
562  << "total error: " << dTotErr << std::endl);
563 
564  if (pRTSolver) {
565  pRTSolver->Log();
566  }
567 
568  return false;
569 
570  } else if (pRTSolver && pRTSolver->IsStopCommanded()) {
571  silent_cout(outputCounterPrefix
572  << "Simulation is stopped by RTAI task" << std::endl
573  << "Simulation ended at time "
574  << dTime << " after "
575  << lStep << " steps;" << std::endl
576  << "total iterations: " << iTotIter << std::endl
577  << "total Jacobian matrices: " << pNLS->TotalAssembledJacobian() << std::endl
578  << "total error: " << dTotErr << std::endl);
579  pRTSolver->Log();
580  return false;
581 
582  } else if (mbdyn_stop_at_end_of_time_step()
583 #ifdef USE_MPI
584  || (MPI_Finalized(&mpi_finalize), mpi_finalize)
585 #endif /* USE_MPI */
586  )
587  {
588  if (pRTSolver) {
590  }
591 
592  silent_cout(outputCounterPrefix
593  << "Interrupted!" << std::endl
594  << "Simulation ended at time "
595  << dTime << " after "
596  << lStep << " steps;" << std::endl
597  << "total iterations: " << iTotIter << std::endl
598  << "total Jacobian matrices: " << pNLS->TotalAssembledJacobian() << std::endl
599  << "total error: " << dTotErr << std::endl);
600 
601  if (pRTSolver) {
602  pRTSolver->Log();
603  }
604 
606  }
607 
608  lStep++;
609 
610  if (pRTSolver) {
611  pRTSolver->Wait();
612  }
613 
614  int retries = -1;
615 IfStepIsToBeRepeated:
616  try {
617  retries++;
619  if (outputStep()) {
620  silent_cout("Step(" << lStep << ':' << retries << ") t=" << dTime + dCurrTimeStep << " dt=" << dCurrTimeStep << std::endl);
621  }
623  CurrStep, pX, pXPrime, pXPrimePrime, pLambda,
625  }
626 
628  if (dCurrTimeStep > dMinTimeStep) {
629  /* Riduce il passo */
630  CurrStep = StepIntegrator::REPEATSTEP;
631  doublereal dOldCurrTimeStep = dCurrTimeStep;
632  dCurrTimeStep = pTSC->dGetNewStepTime(CurrStep, iStIter);
633  if (dCurrTimeStep < dOldCurrTimeStep) {
634  DEBUGCOUT("Changing time step"
635  " during step "
636  << lStep << " after "
637  << iStIter << " iterations"
638  << std::endl);
639  goto IfStepIsToBeRepeated;
640  }
641  }
642 
643  silent_cerr(outputCounterPrefix
644  << "Max iterations number "
645  << std::abs(pRegularSteps->GetIntegratorMaxIters())
646  << " has been reached during"
647  " step " << lStep << ';'
648  << std::endl
649  << "time step dt="
650  << dCurrTimeStep
651  << " cannot be reduced"
652  " further;" << std::endl
653  << "aborting..." << std::endl);
655  }
656 
658  /*
659  * Mettere qui eventuali azioni speciali
660  * da intraprendere in caso di errore ...
661  */
662  silent_cerr(outputCounterPrefix
663  << "Simulation diverged after "
664  << iStIter << " iterations, before "
665  "reaching max iteration number "
667  << " during step " << lStep << ';'
668  << std::endl
669  << "time step dt="
670  << dCurrTimeStep
671  << " cannot be reduced"
672  " further;" << std::endl
673  << "aborting..." << std::endl);
675  }
677  bSolConv = true;
678  }
679 
680  catch (EndOfSimulation& eos) {
681  silent_cerr("Simulation ended during a regular step:\n"
682  << eos.what() << "\n");
684 #ifdef USE_MPI
685  MBDynComm.Abort(0);
686 #endif /* USE_MPI */
687 
688  if (pRTSolver) {
690  }
691 
692  silent_cout(outputCounterPrefix
693  << "Simulation ended at time "
694  << dTime << " after "
695  << lStep << " steps;" << std::endl
696  << "total iterations: " << iTotIter << std::endl
697  << "total Jacobian matrices: " << pNLS->TotalAssembledJacobian() << std::endl
698  << "total error: " << dTotErr << std::endl);
699 
700  if (pRTSolver) {
701  pRTSolver->Log();
702  }
703 
704  return false;
705  }
706  catch (...) {
707  throw;
708  }
709 
710  dTotErr += dTest;
711  iTotIter += iStIter;
712 
714 
715  /* Si fa dare l'std::ostream al file di output per il log */
716  std::ostream& Out = pDM->GetOutFile();
717 
718  if (outputMsg()) {
719  Out << "Step " << lStep
720  << " " << dTime+dCurrTimeStep
721  << " " << dCurrTimeStep
722  << " " << iStIter
723  << " " << dTest
724  << " " << dSolTest
725  << " " << bSolConv
726  << " " << bOut
727  << std::endl;
728  }
729 
730  if (bOutputCounter) {
731  silent_cout("Step " << std::setw(5) << lStep
732  << " " << std::setw(13) << dTime + dCurrTimeStep
733  << " " << std::setw(13) << dCurrTimeStep
734  << " " << std::setw(4) << iStIter
735  << " " << std::setw(13) << dTest
736  << " " << std::setw(13) << dSolTest
737  << " " << bSolConv
738  << " " << bOut
740  }
741 
742  DEBUGCOUT("Step " << lStep
743  << " has been successfully completed "
744  "in " << iStIter << " iterations" << std::endl);
745 
747  dTime += dRefTimeStep;
748 
749  bSolConv = false;
750 
751  /* Calcola il nuovo timestep */
753  DEBUGCOUT("Current time step: " << dCurrTimeStep << std::endl);
754 
755  return true;
756 }
757 
758 /* Distruttore */
760 {
761  DEBUGCOUTFNAME("Solver::~Solver");
762 
763  SAFEDELETE(pX);
767 
768  if (pdWorkSpace != NULL) {
770  }
771 
772  if (pDM != NULL) {
773  SAFEDELETE(pDM);
774  }
775 
776  if (pRegularSteps) {
778  }
779 
780  if (pSM) {
781  SAFEDELETE(pSM);
782  }
783 
784  if (pNLS) {
785  SAFEDELETE(pNLS);
786  }
787 }
788 
789 /* scrive il contributo al file di restart */
790 std::ostream &
791 InverseSolver::Restart(std::ostream& out,DataManager::eRestart type) const
792 {
793 #if 0
794  out << "begin: inverse dynamics;" << std::endl;
795  switch(type) {
796  case DataManager::ATEND:
797  out << " # initial time: " << pDM->dGetTime() << ";"
798  << std::endl
799  << " # final time: " << dFinalTime << ";"
800  << std::endl
801  << " # time step: " << dInitialTimeStep << ";"
802  << std::endl;
803  break;
805  case DataManager::TIME:
806  case DataManager::TIMES:
807  out << " initial time: " << pDM->dGetTime()<< ";" << std::endl
808  << " final time: " << dFinalTime << ";" << std::endl
809  << " time step: " << dInitialTimeStep << ";"
810  << std::endl;
811  break;
812  default:
813  ASSERT(0);
814  }
815 
816  out << " max iterations: " << pRegularSteps->GetIntegratorMaxIters()
817  << ";" << std::endl
818  << " tolerance: " << pRegularSteps->GetIntegratorDTol();
819  switch(ResTest) {
821  out << ", test, norm" ;
822  break;
824  out << ", test, minmax" ;
825  break;
827  NO_OP;
828  default:
829  silent_cerr("unhandled nonlinear solver test type" << std::endl);
831  }
832 
833  if (bScale) {
834  out << ", scale"
835  << ", " << pRegularSteps->GetIntegratorDSolTol();
836  }
837 
838  switch (SolTest) {
840  out << ", test, norm" ;
841  break;
843  out << ", test, minmax" ;
844  break;
846  NO_OP;
847  default:
848  ASSERT(0);
849  }
850  out
851  << ";" << std::endl
852  << " derivatives max iterations: " << pDerivativeSteps->GetIntegratorMaxIters() << ";" << std::endl
853  << " derivatives tolerance: " << pDerivativeSteps->GetIntegratorDTol() << ";" << std::endl
854  << " derivatives coefficient: " << dDerivativesCoef << ";" << std::endl;
855  if (iDummyStepsNumber) {
856  out
857  << " dummy steps max iterations: " << pDummySteps->GetIntegratorMaxIters() << ";" << std::endl
858  << " dummy steps tolerance: " << pDummySteps->GetIntegratorDTol() << ";" << std::endl;
859  }
860  out
861  << " dummy steps number: " << iDummyStepsNumber << ";" << std::endl
862  << " dummy steps ratio: " << dDummyStepsRatio << ";" << std::endl;
863  switch (NonlinearSolverType) {
865  out << " # nonlinear solver: matrix free;" << std::endl;
866  break;
868  default :
869  out << " nonlinear solver: newton raphson";
870  if (!bTrueNewtonRaphson) {
871  out << ", modified, " << iIterationsBeforeAssembly;
872  if (bKeepJac) {
873  out << ", keep jacobian matrix";
874  }
875  if (bHonorJacRequest) {
876  out << ", honor element requests";
877  }
878  }
879  out << ";" << std::endl;
880  }
881  out << " solver: ";
883  out << "end: multistep;" << std::endl << std::endl;
884 #endif
885  return out;
886 }
887 
890 {
891  return ProblemType;
892 }
893 
894 void
896 {
897  switch (iOrder) {
899  dw1 = this->dw1[0];
900  dw2 = this->dw2[0];
901  break;
902 
904  dw1 = this->dw1[1];
905  dw2 = this->dw2[1];
906  break;
907 
909 #if 0
910  dw1 = 0.;
911  dw2 = this->dw1[2] + this->dw2[2];
912 #endif
913  dw1 = this->dw1[2];
914  dw2 = this->dw2[2];
915  break;
916 
917  default:
919  }
920 }
921 
922 /* Dati dell'integratore */
923 void
925 {
926  DEBUGCOUTFNAME("InverseDynamics::ReadData");
927 
928  /* parole chiave */
929  const char* sKeyWords[] = {
930  "begin",
931  "inverse" "dynamics",
932  "end",
933 
934  "initial" "time",
935  "final" "time",
936  "time" "step",
937  "min" "time" "step",
938  "max" "time" "step",
939  "tolerance",
940  "max" "iterations",
941  "modify" "residual" "test",
942  "full" "residual" "test",
943 
944  "abort" "after",
945  "input",
946  "assembly",
947  "derivatives",
948 
949  /* DEPRECATED */ "fictitious" "steps" /* END OF DEPRECATED */ ,
950  "dummy" "steps",
951 
952  "output",
953  "none",
954  "iterations",
955  "residual",
956  "solution",
957  /* DEPRECATED */ "jacobian" /* END OF DEPRECATED */ ,
958  "jacobian" "matrix",
959  "bailout",
960  "messages",
961  "output" "meter",
962 
963  "method",
964  "inverse" "default",
965 
966 
967  /* DEPRECATED */
968  "true",
969  "modified",
970  /* END OF DEPRECATED */
971 
972  "strategy",
973 
974  /* DEPRECATED */
975  "solver",
976  "interface" "solver",
977  /* END OF DEPRECATED */
978  "linear" "solver",
979  "interface" "linear" "solver",
980 
981  /* DEPRECATED */
982  "preconditioner",
983  /* END OF DEPRECATED */
984 
985  "nonlinear" "solver",
986  "default",
987  "newton" "raphson",
988  "matrix" "free",
989  "bicgstab",
990  "gmres",
991  /* DEPRECATED */ "full" "jacobian" /* END OF DEPRECATED */ ,
992  "full" "jacobian" "matrix",
993 
994  /* RTAI stuff */
995  "real" "time",
996 
997  /* multithread stuff */
998  "threads",
999 
1000  "problem" "type",
1001  "fully" "actuated" "collocated",
1002  "fully" "actuated" "non" "collocated",
1003  "underdetermined" "underactuated" "collocated",
1004  "underdetermined" "fully" "actuated",
1005 
1006  NULL
1007  };
1008 
1009  /* enum delle parole chiave */
1010  enum KeyWords {
1011  UNKNOWN = -1,
1012  BEGIN = 0,
1013  INVERSEDYNAMICS,
1014  END,
1015 
1016  INITIALTIME,
1017  FINALTIME,
1018  TIMESTEP,
1019  MINTIMESTEP,
1020  MAXTIMESTEP,
1021  TOLERANCE,
1022  MAXITERATIONS,
1023  MODIFY_RES_TEST,
1024  FULL_RES_TEST,
1025 
1026  ABORTAFTER,
1027  INPUT,
1028  ASSEMBLY,
1029  DERIVATIVES,
1030  FICTITIOUSSTEPS,
1031  DUMMYSTEPS,
1032 
1033  OUTPUT,
1034  NONE,
1035  ITERATIONS,
1036  RESIDUAL,
1037  SOLUTION,
1038  JACOBIAN,
1039  JACOBIANMATRIX,
1040  BAILOUT,
1041  MESSAGES,
1042  OUTPUTMETER,
1043 
1044  METHOD,
1045  INVERSEDEFAULT,
1046 
1047  /* DEPRECATED */
1048  NR_TRUE,
1049  MODIFIED,
1050  /* END OF DEPRECATED */
1051 
1052  STRATEGY,
1053 
1054  /* DEPRECATED */
1055  SOLVER,
1056  INTERFACESOLVER,
1057  /* END OF DEPRECATED */
1058  LINEARSOLVER,
1059  INTERFACELINEARSOLVER,
1060 
1061  /* DEPRECATED */
1062  PRECONDITIONER,
1063  /* END OF DEPRECATED */
1064 
1065  NONLINEARSOLVER,
1066  DEFAULT,
1067  NEWTONRAPHSON,
1068  MATRIXFREE,
1069  BICGSTAB,
1070  GMRES,
1071  FULLJACOBIAN,
1072  FULLJACOBIANMATRIX,
1073 
1074  /* RTAI stuff */
1075  REALTIME,
1076 
1077  THREADS,
1078 
1079  PROBLEMTYPE,
1080  FAC,
1081  FANC,
1082  UDUAC,
1083  UDFA,
1084 
1085  LASTKEYWORD
1086  };
1087 
1088  /* tabella delle parole chiave */
1089  KeyTable K(HP, sKeyWords);
1090 
1091  /* legge i dati della simulazione */
1092  if (KeyWords(HP.GetDescription()) != BEGIN) {
1093  silent_cerr("Error: <begin> expected at line "
1094  << HP.GetLineData() << "; aborting..." << std::endl);
1096  }
1097 
1098  if (KeyWords(HP.GetWord()) != INVERSEDYNAMICS) {
1099  silent_cerr("Error: <begin: inverse dynamics;> expected at line "
1100  << HP.GetLineData() << "; aborting..." << std::endl);
1102  }
1103 
1104  /* dati letti qui ma da passare alle classi
1105  * StepIntegration e NonlinearSolver
1106  */
1107 
1108  doublereal dTol = ::dDefaultTol;
1109  doublereal dSolutionTol = 0.;
1111  bool bModResTest(false);
1112 
1113 #ifdef USE_MULTITHREAD
1114  bool bSolverThreads(false);
1115  unsigned nSolverThreads(0);
1116 #endif /* USE_MULTITHREAD */
1117 
1118 
1119  /* Ciclo infinito */
1120  while (true) {
1121  KeyWords CurrKeyWord = KeyWords(HP.GetDescription());
1122 
1123  switch (CurrKeyWord) {
1124  case INITIALTIME:
1125  dInitialTime = HP.GetReal();
1126  DEBUGLCOUT(MYDEBUG_INPUT, "Initial time is "
1127  << dInitialTime << std::endl);
1128  break;
1129 
1130  case FINALTIME:
1131  dFinalTime = HP.GetReal();
1132  DEBUGLCOUT(MYDEBUG_INPUT, "Final time is "
1133  << dFinalTime << std::endl);
1134 
1135  if(dFinalTime <= dInitialTime) {
1136  silent_cerr("warning: final time " << dFinalTime
1137  << " is less than initial time "
1138  << dInitialTime << ';' << std::endl
1139  << "this will cause the simulation"
1140  " to abort" << std::endl);
1141  }
1142  break;
1143 
1144  case TIMESTEP:
1145  dInitialTimeStep = HP.GetReal();
1146  DEBUGLCOUT(MYDEBUG_INPUT, "Initial time step is "
1147  << dInitialTimeStep << std::endl);
1148 
1149  if (dInitialTimeStep == 0.) {
1150  silent_cerr("warning, null initial time step"
1151  " is not allowed" << std::endl);
1152  } else if (dInitialTimeStep < 0.) {
1154  silent_cerr("warning, negative initial time step"
1155  " is not allowed;" << std::endl
1156  << "its modulus " << dInitialTimeStep
1157  << " will be considered" << std::endl);
1158  }
1159  break;
1160 
1161  case MINTIMESTEP:
1162  dMinTimeStep = HP.GetReal();
1163  DEBUGLCOUT(MYDEBUG_INPUT, "Min time step is "
1164  << dMinTimeStep << std::endl);
1165 
1166  if (dMinTimeStep == 0.) {
1167  silent_cerr("warning, null minimum time step"
1168  " is not allowed" << std::endl);
1170  } else if (dMinTimeStep < 0.) {
1172  silent_cerr("warning, negative minimum time step"
1173  " is not allowed;" << std::endl
1174  << "its modulus " << dMinTimeStep
1175  << " will be considered" << std::endl);
1176  }
1177  break;
1178 
1179  case MAXTIMESTEP:
1181 #ifdef DEBUG
1182  if (typeid(*MaxTimeStep.pGetDriveCaller()) == typeid(PostponedDriveCaller)) {
1183  DEBUGLCOUT(MYDEBUG_INPUT, "Max time step is postponed" << std::endl);
1184 
1185  } else {
1186  DEBUGLCOUT(MYDEBUG_INPUT, "Max time step is " << MaxTimeStep.dGet() << std::endl);
1187  }
1188 #endif // DEBUG
1189 
1190  if (dGetInitialMaxTimeStep() == 0.) {
1191  silent_cout("no max time step limit will be"
1192  " considered" << std::endl);
1193  } else if (dGetInitialMaxTimeStep() < 0.) {
1194  silent_cerr("negative max time step"
1195  " is not allowed" << std::endl);
1197  }
1198  break;
1199 
1200  case ABORTAFTER: {
1201  KeyWords WhenToAbort(KeyWords(HP.GetWord()));
1202  switch (WhenToAbort) {
1203  case INPUT:
1206  "Simulation will abort after"
1207  " data input" << std::endl);
1208  break;
1209 
1210  case ASSEMBLY:
1213  "Simulation will abort after"
1214  " initial assembly" << std::endl);
1215  break;
1216 
1217  case DERIVATIVES:
1220  "Simulation will abort after"
1221  " derivatives solution" << std::endl);
1222  break;
1223 
1224  case FICTITIOUSSTEPS:
1225  case DUMMYSTEPS:
1228  "Simulation will abort after"
1229  " dummy steps solution" << std::endl);
1230  break;
1231 
1232  default:
1233  silent_cerr("Don't know when to abort,"
1234  " so I'm going to abort now" << std::endl);
1236  }
1237  break;
1238  }
1239 
1240  case OUTPUT: {
1241  unsigned OF = OUTPUT_DEFAULT;
1242  bool setOutput = false;
1243 
1244  while (HP.IsArg()) {
1245  KeyWords OutputFlag(KeyWords(HP.GetWord()));
1246  switch (OutputFlag) {
1247  case NONE:
1248  OF = OUTPUT_NONE;
1249  setOutput = true;
1250  break;
1251 
1252  case ITERATIONS:
1253  OF |= OUTPUT_ITERS;
1254  break;
1255 
1256  case RESIDUAL:
1257  OF |= OUTPUT_RES;
1258  break;
1259 
1260  case SOLUTION:
1261  OF |= OUTPUT_SOL;
1262  break;
1263 
1264  case JACOBIAN:
1265  case JACOBIANMATRIX:
1266  OF |= OUTPUT_JAC;
1267  break;
1268 
1269  case BAILOUT:
1270  OF |= OUTPUT_BAILOUT;
1271  break;
1272 
1273  case MESSAGES:
1274  OF |= OUTPUT_MSG;
1275  break;
1276 
1277  default:
1278  silent_cerr("Unknown output flag "
1279  "at line " << HP.GetLineData()
1280  << "; ignored" << std::endl);
1281  break;
1282  }
1283  }
1284 
1285  if (setOutput) {
1286  SetOutputFlags(OF);
1287  } else {
1288  AddOutputFlags(OF);
1289  }
1290 
1291  break;
1292  }
1293 
1294  case OUTPUTMETER:
1295  SetOutputMeter(HP.GetDriveCaller(true));
1296  break;
1297 
1298  case TOLERANCE: {
1299  /*
1300  * residual tolerance; can be the keyword "null",
1301  * which means that the convergence test
1302  * will be computed on the solution, or a number
1303  */
1304  if (HP.IsKeyWord("null")) {
1305  dTol = 0.;
1306 
1307  } else {
1308  dTol = HP.GetReal();
1309  if (dTol < 0.) {
1310  dTol = ::dDefaultTol;
1311  silent_cerr("warning, residual tolerance "
1312  "< 0. is illegal; "
1313  "using default value " << dTol
1314  << std::endl);
1315  }
1316  }
1317 
1318  /* safe default */
1319  if (dTol == 0.) {
1321  }
1322 
1323  if (HP.IsArg()) {
1324  if (HP.IsKeyWord("test")) {
1325  if (HP.IsKeyWord("norm")) {
1327  } else if (HP.IsKeyWord("minmax")) {
1329  } else if (HP.IsKeyWord("none")) {
1331  } else {
1332  silent_cerr("unknown test "
1333  "method at line "
1334  << HP.GetLineData()
1335  << std::endl);
1337  }
1338 
1339  if (HP.IsKeyWord("scale")) {
1341  silent_cerr("it's a nonsense "
1342  "to scale a disabled test; "
1343  "\"scale\" ignored"
1344  << std::endl);
1345  bScale = false;
1346  } else {
1347  bScale = true;
1348  }
1349  }
1350  }
1351  }
1352 
1353  if (HP.IsArg()) {
1354  if (!HP.IsKeyWord("null")) {
1355  dSolutionTol = HP.GetReal();
1356  }
1357 
1358  /* safe default */
1359  if (dSolutionTol != 0.) {
1361  }
1362 
1363  if (HP.IsArg()) {
1364  if (HP.IsKeyWord("test")) {
1365  if (HP.IsKeyWord("norm")) {
1367  } else if (HP.IsKeyWord("minmax")) {
1369  } else if (HP.IsKeyWord("none")) {
1371  } else {
1372  silent_cerr("unknown test "
1373  "method at line "
1374  << HP.GetLineData()
1375  << std::endl);
1377  }
1378  }
1379  }
1380 
1381  } else if (dTol == 0.) {
1382  silent_cerr("need solution tolerance "
1383  "with null residual tolerance"
1384  << std::endl);
1386  }
1387 
1388  if (dSolutionTol < 0.) {
1389  dSolutionTol = 0.;
1390  silent_cerr("warning, solution tolerance "
1391  "< 0. is illegal; "
1392  "solution test is disabled"
1393  << std::endl);
1394  }
1395 
1396  if (dTol == 0. && dSolutionTol == 0.) {
1397  silent_cerr("both residual and solution "
1398  "tolerances are zero" << std::endl);
1400  }
1401 
1402  DEBUGLCOUT(MYDEBUG_INPUT, "tolerance = " << dTol
1403  << ", " << dSolutionTol << std::endl);
1404  break;
1405  }
1406 
1407 
1408  case MAXITERATIONS: {
1409  iMaxIterations = HP.GetInt();
1410  if (iMaxIterations < 1) {
1412  silent_cerr("warning, max iterations "
1413  "< 1 is illegal; using default value "
1414  << iMaxIterations
1415  << std::endl);
1416  }
1418  "Max iterations = "
1419  << iMaxIterations << std::endl);
1420  if (HP.IsKeyWord("at" "most")) {
1422  }
1423  break;
1424  }
1425 
1426  case MODIFY_RES_TEST:
1427  if (bParallel) {
1428  silent_cerr("\"modify residual test\" "
1429  "not supported by schur data manager "
1430  "at line " << HP.GetLineData()
1431  << "; ignored" << std::endl);
1432  } else {
1433  bModResTest = true;
1435  "Modify residual test" << std::endl);
1436  }
1437  break;
1438 
1439  case FULL_RES_TEST:
1440  if (HP.IsArg()) {
1442 
1443  } else {
1444  bFullResTest = true;
1445  }
1446 
1448  "Full residual test: " << bFullResTest << std::endl);
1449  break;
1450 
1451  case NEWTONRAPHSON: {
1452  pedantic_cout("Newton Raphson is deprecated; use "
1453  "\"nonlinear solver: newton raphson "
1454  "[ , modified, <steps> ]\" instead"
1455  << std::endl);
1456  KeyWords NewRaph(KeyWords(HP.GetWord()));
1457  switch(NewRaph) {
1458  case MODIFIED:
1459  bTrueNewtonRaphson = 0;
1460  if (HP.IsArg()) {
1462  } else {
1464  }
1465  DEBUGLCOUT(MYDEBUG_INPUT, "Modified "
1466  "Newton-Raphson will be used; "
1467  "matrix will be assembled "
1468  "at most after "
1470  << " iterations" << std::endl);
1471  break;
1472 
1473  default:
1474  silent_cerr("warning: unknown case; "
1475  "using default" << std::endl);
1476 
1477  /* no break: fall-thru to next case */
1478  case NR_TRUE:
1479  bTrueNewtonRaphson = 1;
1481  break;
1482  }
1483  break;
1484  }
1485 
1486  case END:
1487  if (KeyWords(HP.GetWord()) != INVERSEDYNAMICS) {
1488  silent_cerr("\"end: inverse dynamics;\" expected "
1489  "at line " << HP.GetLineData()
1490  << "; aborting..." << std::endl);
1492  }
1493  goto EndOfCycle;
1494 
1495  case STRATEGY:
1496  pTSC = ReadTimeStepData(this,HP);
1497 
1498  break;
1499 
1500  case SOLVER:
1501  silent_cerr("\"solver\" keyword at line "
1502  << HP.GetLineData()
1503  << " is deprecated; "
1504  "use \"linear solver\" instead"
1505  << std::endl);
1506  case LINEARSOLVER:
1508  break;
1509 
1510  case INTERFACESOLVER:
1511  silent_cerr("\"interface solver\" keyword at line "
1512  << HP.GetLineData()
1513  << " is deprecated; "
1514  "use \"interface linear solver\" "
1515  "instead" << std::endl);
1516  case INTERFACELINEARSOLVER:
1517  ReadLinSol(CurrIntSolver, HP, true);
1518 
1519 #ifndef USE_MPI
1520  silent_cerr("Interface solver only allowed "
1521  "when compiled with MPI support" << std::endl);
1522 #endif /* ! USE_MPI */
1523  break;
1524 
1525  case NONLINEARSOLVER:
1526  switch (KeyWords(HP.GetWord())) {
1527  case DEFAULT:
1529  break;
1530 
1531  case NEWTONRAPHSON:
1533  break;
1534 
1535  case MATRIXFREE:
1537  break;
1538 
1539  default:
1540  silent_cerr("unknown nonlinear solver "
1541  "at line " << HP.GetLineData()
1542  << std::endl);
1544  break;
1545  }
1546 
1547  switch (NonlinearSolverType) {
1549  bTrueNewtonRaphson = true;
1550  bKeepJac = false;
1552 
1553  if (HP.IsKeyWord("modified")) {
1554  bTrueNewtonRaphson = false;
1556 
1557  if (HP.IsKeyWord("keep" "jacobian")) {
1558  pedantic_cout("Use of deprecated \"keep jacobian\" at line " << HP.GetLineData() << std::endl);
1559  bKeepJac = true;
1560 
1561  } else if (HP.IsKeyWord("keep" "jacobian" "matrix")) {
1562  bKeepJac = true;
1563  }
1564 
1565  DEBUGLCOUT(MYDEBUG_INPUT, "modified "
1566  "Newton-Raphson "
1567  "will be used; "
1568  "matrix will be "
1569  "assembled at most "
1570  "after "
1572  << " iterations"
1573  << std::endl);
1574  if (HP.IsKeyWord("honor" "element" "requests")) {
1575  bHonorJacRequest = true;
1577  "honor elements' "
1578  "request to update "
1579  "the preconditioner"
1580  << std::endl);
1581  }
1582  }
1583  break;
1584 
1586  switch (KeyWords(HP.GetWord())) {
1587  case DEFAULT:
1589  break;
1590 
1591 
1592  case BICGSTAB:
1594  break;
1595 
1596  case GMRES:
1598  break;
1599 
1600  default:
1601  silent_cerr("unknown iterative "
1602  "solver at line "
1603  << HP.GetLineData()
1604  << std::endl);
1606  }
1607 
1608  if (HP.IsKeyWord("tolerance")) {
1609  dIterTol = HP.GetReal();
1610  DEBUGLCOUT(MYDEBUG_INPUT,"inner "
1611  "iterative solver "
1612  "tolerance: "
1613  << dIterTol
1614  << std::endl);
1615  }
1616 
1617  if (HP.IsKeyWord("steps")) {
1618  iIterativeMaxSteps = HP.GetInt();
1619  DEBUGLCOUT(MYDEBUG_INPUT, "maximum "
1620  "number of inner "
1621  "steps for iterative "
1622  "solver: "
1624  << std::endl);
1625  }
1626 
1627  if (HP.IsKeyWord("tau")) {
1628  dIterertiveTau = HP.GetReal();
1630  "tau scaling "
1631  "coefficient "
1632  "for iterative "
1633  "solver: "
1634  << dIterertiveTau
1635  << std::endl);
1636  }
1637 
1638  if (HP.IsKeyWord("eta")) {
1639  dIterertiveEtaMax = HP.GetReal();
1640  DEBUGLCOUT(MYDEBUG_INPUT, "maximum "
1641  "eta coefficient "
1642  "for iterative "
1643  "solver: "
1645  << std::endl);
1646  }
1647 
1648  if (HP.IsKeyWord("preconditioner")) {
1649  KeyWords KPrecond = KeyWords(HP.GetWord());
1650  switch (KPrecond) {
1651  case FULLJACOBIAN:
1652  case FULLJACOBIANMATRIX:
1654  if (HP.IsKeyWord("steps")) {
1655  iPrecondSteps = HP.GetInt();
1657  "number of steps "
1658  "before recomputing "
1659  "the preconditioner: "
1660  << iPrecondSteps
1661  << std::endl);
1662  }
1663  if (HP.IsKeyWord("honor" "element" "requests")) {
1664  bHonorJacRequest = true;
1666  "honor elements' "
1667  "request to update "
1668  "the preconditioner"
1669  << std::endl);
1670  }
1671  break;
1672 
1673  /* add other preconditioners
1674  * here */
1675 
1676  default:
1677  silent_cerr("unknown "
1678  "preconditioner "
1679  "at line "
1680  << HP.GetLineData()
1681  << std::endl);
1683  }
1684  break;
1685  }
1686  break;
1687 
1688  default:
1689  ASSERT(0);
1691  }
1692  break;
1693 
1694  case REALTIME:
1695  pRTSolver = ReadRTSolver(this, HP);
1696  break;
1697 
1698  case THREADS:
1699  if (HP.IsKeyWord("auto")) {
1700 #ifdef USE_MULTITHREAD
1701  int n = get_nprocs();
1702  /* sanity checks ... */
1703  if (n <= 0) {
1704  silent_cerr("got " << n << " CPUs "
1705  "at line "
1706  << HP.GetLineData()
1707  << std::endl);
1708  nThreads = 1;
1709  } else {
1710  nThreads = n;
1711  }
1712 #else /* ! USE_MULTITHREAD */
1713  silent_cerr("configure with "
1714  "--enable-multithread "
1715  "for multithreaded assembly"
1716  << std::endl);
1717 #endif /* ! USE_MULTITHREAD */
1718 
1719  } else if (HP.IsKeyWord("disable")) {
1720 #ifdef USE_MULTITHREAD
1721  nThreads = 1;
1722 #endif /* USE_MULTITHREAD */
1723 
1724  } else {
1725 #ifdef USE_MULTITHREAD
1726  bool bAssembly = false;
1727  bool bSolver = false;
1728  bool bAll = true;
1729  unsigned nt;
1730 #endif // USE_MULTITHREAD
1731 
1732  if (HP.IsKeyWord("assembly")) {
1733 #ifdef USE_MULTITHREAD
1734  bAll = false;
1735  bAssembly = true;
1736 #endif // USE_MULTITHREAD
1737 
1738  } else if (HP.IsKeyWord("solver")) {
1739 #ifdef USE_MULTITHREAD
1740  bAll = false;
1741  bSolver = true;
1742 #endif // USE_MULTITHREAD
1743  }
1744 
1745 #ifdef USE_MULTITHREAD
1746  nt =
1747 #endif // USE_MULTITHREAD
1748  HP.GetInt();
1749 
1750 #ifdef USE_MULTITHREAD
1751  if (bAll || bAssembly) {
1752  nThreads = nt;
1753  }
1754 
1755  if (bAll || bSolver) {
1756  bSolverThreads = true;
1757  nSolverThreads = nt;
1758  }
1759 #else /* ! USE_MULTITHREAD */
1760  silent_cerr("configure with "
1761  "--enable-multithread "
1762  "for multithreaded assembly"
1763  << std::endl);
1764 #endif /* ! USE_MULTITHREAD */
1765  }
1766  break;
1767 
1768  case PROBLEMTYPE:
1769  switch (KeyWords(HP.GetWord())) {
1770  case FAC:
1771  // default
1773  break;
1774 
1775  case FANC:
1777  break;
1778 
1779  case UDUAC:
1781  if (HP.IsKeyWord("weights")) {
1782  for (int iCnt = 0; iCnt <= InverseDynamics::VELOCITY; iCnt++) {
1783  dw1[iCnt] = HP.GetReal();
1784  dw2[iCnt] = HP.GetReal();
1785  }
1788 
1789  } else {
1790  for (int iCnt = 0; iCnt <= InverseDynamics::ACCELERATION; iCnt++) {
1791  dw1[iCnt] = 1.;
1792  dw2[iCnt] = 0.;
1793  }
1794  }
1795 
1796  break;
1797 
1798  case UDFA:
1800  if (HP.IsKeyWord("weights")) {
1801  for (int iCnt = 0; iCnt <= InverseDynamics::ACCELERATION; iCnt++) {
1802  dw1[iCnt] = HP.GetReal();
1803  dw2[iCnt] = HP.GetReal();
1804  }
1805 
1806  } else {
1807  for (int iCnt = 0; iCnt <= InverseDynamics::ACCELERATION; iCnt++) {
1808  dw1[iCnt] = 1.;
1809  dw2[iCnt] = 0.;
1810  }
1811  }
1812  break;
1813 
1814  default:
1815  silent_cerr("inverse dynamics: unrecognized problem type at line " << HP.GetLineData() << std::endl);
1817  }
1818 
1819  if (HP.IsArg()) {
1820  silent_cerr("inverse dynamics: semicolon expected at line " << HP.GetLineData() << std::endl);
1822  }
1823  break;
1824 
1825  default:
1826  silent_cerr("unknown description at line "
1827  << HP.GetLineData() << "; aborting..."
1828  << std::endl);
1830  }
1831  }
1832 
1833 EndOfCycle: /* esce dal ciclo di lettura */
1834 
1835  if (MaxTimeStep.pGetDriveCaller() == 0) {
1836  DriveCaller *pDC = 0;
1838  MaxTimeStep.Set(pDC);
1839  }
1840 
1841  if (typeid(*MaxTimeStep.pGetDriveCaller()) == typeid(ConstDriveCaller)) {
1843  }
1844 
1845  if (pTSC == 0) {
1846  pedantic_cout("No time step control strategy defined; defaulting to NoChange time step control" );
1847  pTSC = new NoChange(this);
1848  }
1849 
1850  if (dFinalTime < dInitialTime) {
1852  }
1853 
1854  if (dFinalTime == dInitialTime) {
1856  }
1857 
1861  dTol,
1862  dSolutionTol,
1863  0, /* Previous Steps to be used */
1864  1, /* Unknown States: FIXME: 2? 3? */
1865  bModResTest));
1866 
1867 #ifdef USE_MULTITHREAD
1868  if (bSolverThreads) {
1869  if (CurrLinearSolver.SetNumThreads(nSolverThreads)) {
1870  silent_cerr("linear solver "
1872  << " does not support "
1873  "threaded solution" << std::endl);
1874  }
1875  }
1876 #endif /* USE_MULTITHREAD */
1877 }
1878 
bool bFullResTest
Definition: invsolver.h:75
virtual bool Start(void)
Definition: invsolver.cc:517
void SetScale(VectorHandler &XScale) const
Definition: dofman.cc:159
integer iUnkStates
Definition: solver.h:229
void ReadLinSol(LinSol &cs, HighParser &HP, bool bAllowEmpty)
Definition: readlinsol.cc:39
integer * GetDofsList(DofType who) const
MyVectorHandler * pX
Definition: solver.h:233
virtual void IDSetTest(NonlinearSolverTestRange *pResTest, NonlinearSolverTestRange *pSolTest, bool bFullResTest)
Definition: invdataman.cc:1111
integer iNumIntDofs
Definition: solver.h:319
#define DEBUG_LEVEL_MATCH(level)
Definition: myassert.h:241
std::string outputCounterPrefix
Definition: solver.h:362
std::ostream & RestartLinSol(std::ostream &out, const LinSol &cs)
Definition: readlinsol.cc:611
RTSolverBase * ReadRTSolver(Solver *pS, MBDynParser &HP)
Definition: rtsolver.cc:206
unsigned GetSolverFlags(void) const
Definition: linsol.cc:257
StepIntegrator * pCurrStepIntegrator
Definition: solver.h:283
bool bParallel
Definition: solver.h:315
virtual void SetScale(const VectorHandler *pScl)
Definition: nonlin.cc:278
bool outputCounter(void) const
bool bOutputCounter
Definition: solver.h:361
#define MBDYN_EXCEPT_ARGS
Definition: except.h:63
integer * pLocDofs
Definition: solver.h:320
void SetOutputMeter(DriveCaller *pOM)
#define DEBUGCOUTFNAME(fname)
Definition: myassert.h:256
virtual integer GetInt(integer iDefval=0)
Definition: parser.cc:1050
bool bScale
Definition: solver.h:297
void SetOutputFlags(unsigned OF)
virtual bool Advance(void)
Definition: invsolver.cc:538
virtual void OutputTypes(const bool fpred)
Definition: stepsol.cc:111
virtual void Init(void)
Definition: rtsolver.cc:69
doublereal dRefTimeStep
Definition: solver.h:240
doublereal dGetTime(void) const
Definition: dataman2.cc:165
doublereal dw2[3]
Definition: invsolver.h:71
doublereal dMinTimeStep
Definition: solver.h:111
virtual bool Output(long lStep, const doublereal &dTime, const doublereal &dTimeStep, bool force=false) const
Definition: dataman2.cc:2355
const char * what(void) const
Definition: except.cc:54
std::string sOutputFileName
Definition: solver.h:114
virtual void Wait(void)=0
integer iTotIter
Definition: solver.h:364
bool outputMsg(void) const
#define SAFEDELETEARR(pnt)
Definition: mynewmem.h:713
StepIntegrator::StepChange CurrStep
Definition: solver.h:112
integer iDummyStepsNumber
Definition: solver.h:253
doublereal dDummyStepsRatio
Definition: solver.h:254
virtual void SetDriveHandler(const DriveHandler *driveHandler)=0
bool SetNumThreads(unsigned nt)
Definition: linsol.cc:314
StepIntegrator * pDerivativeSteps
Definition: solver.h:278
virtual integer GetIntegratorMaxIters(void) const
Definition: stepsol.cc:93
integer iMaxIterations
Definition: solver.h:116
const DriveHandler * pGetDrvHdl(void) const
Definition: dataman.h:340
doublereal dTime
Definition: solver.h:109
doublereal dw1[3]
Definition: invsolver.h:70
integer iNumPreviousVectors
Definition: solver.h:228
integer iStIter
Definition: solver.h:108
void SetDataManager(DataManager *pDatMan)
Definition: stepsol.cc:74
#define NO_OP
Definition: myassert.h:74
NonlinearSolver::Type NonlinearSolverType
Definition: solver.h:304
InverseSolver(MBDynParser &HP, const std::string &sInputFileName, const std::string &sOutputFileName, unsigned int nThreads, bool bParallel=false)
Definition: invsolver.cc:79
virtual ~InverseSolver(void)
Definition: invsolver.cc:759
enum Solver::@38 eStatus
integer iIterationsBeforeAssembly
Definition: solver.h:303
InverseDynamics::Type GetProblemType(void) const
Definition: invsolver.cc:889
bool outputStep(void) const
Definition: dataman4.cc:96
virtual void StopCommanded(void)=0
std::deque< MyVectorHandler * > qXPrime
Definition: solver.h:232
virtual bool GetYesNoOrBool(bool bDefval=false)
Definition: parser.cc:1038
DriveCaller * pOutputMeter
std::string sInputFileName
Definition: solver.h:113
virtual integer TotalAssembledJacobian(void)
Definition: nonlin.cc:468
LinSol CurrLinearSolver
Definition: solver.h:292
void SetupSolmans(integer iStates, bool bCanBeParallel=false)
Definition: solver.cc:5135
void mbdyn_signal_init(int pre)
Definition: solver.cc:208
bool SetSolver(SolverType t, unsigned f=SOLVER_FLAGS_NONE)
Definition: linsol.cc:181
MatrixFreeSolver::SolverType MFSolverType
Definition: solver.h:305
static const integer iDefaultMaxIterations
Definition: solver_impl.h:99
void SetValue(VectorHandler &X, VectorHandler &XP)
Definition: dataman2.cc:1404
SolutionManager * pSM
Definition: solver.h:332
void LinkToSolution(VectorHandler &XCurr, VectorHandler &XPrimeCurr)
Definition: dataman2.cc:172
virtual void SetDrvHdl(const DriveHandler *pDH)
Definition: drive.cc:487
NonlinearSolverTest::Type SolTest
Definition: solver.h:296
std::ostream & GetOutFile(void) const
Definition: dataman.h:325
std::ostream & Restart(std::ostream &out, DataManager::eRestart type) const
Definition: invsolver.cc:791
doublereal dSolTest
Definition: solver.h:368
void AddOutputFlags(unsigned OF)
#define SAFENEW(pnt, item)
Definition: mynewmem.h:695
static const integer iDefaultIterationsBeforeAssembly
Definition: solver_impl.h:92
DataManager * pDM
Definition: solver.h:327
bool AddSolverFlags(unsigned f)
Definition: linsol.cc:292
Definition: mbdyn.h:76
NonlinearSolver *const AllocateNonlinearSolver()
Definition: solver.cc:5073
virtual bool IsKeyWord(const char *sKeyWord)
Definition: parser.cc:910
bool bTrueNewtonRaphson
Definition: solver.h:301
StepIntegrator * pRegularSteps
Definition: solver.h:282
bool bKeepJac
Definition: solver.h:302
doublereal dInitialTimeStep
Definition: solver.h:241
static const doublereal dDefaultMaxTimeStep
Definition: solver_impl.h:101
doublereal dCurrTimeStep
Definition: solver.h:107
int GetDescription(void)
Definition: parser.cc:730
virtual bool Prepare(void)
Definition: invsolver.cc:93
virtual doublereal GetIntegratorDTol(void) const
Definition: stepsol.cc:99
RTSolverBase * pRTSolver
Definition: solver.h:225
#define DEBUGCOUT(msg)
Definition: myassert.h:232
virtual void OutputPrepare(void)
Definition: dataman2.cc:1451
doublereal dFinalTime
Definition: solver.h:239
void ReadData(MBDynParser &HP)
Definition: invsolver.cc:924
MyVectorHandler * pLambda
Definition: invsolver.h:73
virtual void SetTest(NonlinearSolverTest *pr, NonlinearSolverTest *ps)
Definition: nonlin.cc:455
virtual bool IsStopCommanded(void)
Definition: rtsolver.cc:81
Dof * pDofs
Definition: solver.h:322
doublereal dIterertiveTau
Definition: solver.h:311
virtual MathParser & GetMathParser(void)
Definition: parser.cc:668
doublereal dIterertiveEtaMax
Definition: solver.h:310
virtual void Log(void)=0
int mbdyn_stop_at_end_of_time_step(void)
Definition: solver.cc:182
#define ASSERT(expression)
Definition: colamd.c:977
StepIntegrator * pDummySteps
Definition: solver.h:280
KeyWords
Definition: dataman4.cc:94
virtual void Reset(void)
Definition: vh.cc:459
#define SAFENEWWITHCONSTRUCTOR(pnt, item, constructor)
Definition: mynewmem.h:698
DriveCaller * pGetDriveCaller(void) const
Definition: drive.cc:658
InverseDynamics::Type ProblemType
Definition: invsolver.h:69
void SetTime(const doublereal &dTime, const doublereal &dTimeStep=-1., const integer &iStep=-1, bool bServePending=true)
Definition: dataman2.cc:139
virtual doublereal dGetNewStepTime(StepIntegrator::StepChange currStep, doublereal iPerformedIters)=0
DriveOwner MaxTimeStep
Definition: solver.h:110
virtual DriveCaller * pCopy(void) const =0
Definition: solver.h:78
virtual void SetDriveHandler(const DriveHandler *pDH)
Definition: stepsol.cc:117
static const doublereal dDefaultTol
Definition: solver_impl.h:95
integer iPrecondSteps
Definition: solver.h:308
Dof * pGetDofsList(void) const
virtual bool IsArg(void)
Definition: parser.cc:807
std::deque< MyVectorHandler * > qX
Definition: solver.h:231
doublereal dTotErr
Definition: solver.h:366
integer HowManyDofs(DofType who) const
virtual int GetWord(void)
Definition: parser.cc:1083
NonlinearSolverTest::Type ResTest
Definition: solver.h:295
integer * pIntDofs
Definition: solver.h:321
virtual integer GetIntegratorNumPreviousStates(void) const
Definition: stepsol.cc:81
TimeStepControl * pTSC
Definition: solver.h:106
std::ostream & GetLogFile(void) const
Definition: dataman.h:326
std::string outputCounterPostfix
Definition: solver.h:363
Preconditioner::PrecondType PcType
Definition: solver.h:307
void Set(const DriveCaller *pDC)
Definition: drive.cc:647
doublereal dGet(const doublereal &dVar) const
Definition: drive.cc:664
MyVectorHandler * pXPrime
Definition: solver.h:234
integer iNumLocDofs
Definition: solver.h:318
NonlinearSolver * pNLS
Definition: solver.h:333
int get_nprocs(void)
Definition: get_nprocs.c:44
doublereal dInitialTime
Definition: solver.h:238
TimeStepControl * ReadTimeStepData(Solver *s, MBDynParser &HP)
#define SAFENEWARR(pnt, item, sz)
Definition: mynewmem.h:701
bool bOut
Definition: solver.h:370
SchurDataManager * pSDM
Definition: solver.h:316
DriveCaller * GetDriveCaller(bool bDeferred=false)
Definition: mbpar.cc:2033
#define EMPTY
Definition: colamd.c:768
doublereal dGetInitialMaxTimeStep() const
Definition: solver.cc:4990
MyVectorHandler Scale
Definition: solver.h:298
virtual void Init(integer iMaxIterations, doublereal dMinTimeStep, const DriveOwner &MaxTimeStep, doublereal dInitialTimeStep)=0
AbortAfter eAbortAfter
Definition: solver.h:264
const char *const GetSolverName(void) const
Definition: linsol.cc:269
double doublereal
Definition: colamd.c:52
bool bSolConv
Definition: solver.h:369
void GetWeight(InverseDynamics::Order iOrder, doublereal &dw1, doublereal &dw2) const
Definition: invsolver.cc:895
integer iIterativeMaxSteps
Definition: solver.h:309
long int integer
Definition: colamd.c:51
virtual HighParser::ErrOut GetLineData(void) const
Definition: parsinc.cc:697
doublereal * pdWorkSpace
Definition: solver.h:230
integer iNumDofs
Definition: solver.h:329
Table & GetSymbolTable(void) const
Definition: mathp.cc:1931
doublereal dTest
Definition: solver.h:367
virtual void Close(void)
Definition: parsinc.cc:155
virtual void Setup(void)=0
doublereal dDerivativesCoef
Definition: solver.h:290
void CreatePartition(void)
#define DEBUGLCOUT(level, msg)
Definition: myassert.h:244
long lStep
Definition: solver.h:371
MBDynParser & HP
Definition: solver.h:115
void mbdyn_set_stop_at_end_of_time_step(void)
Definition: solver.cc:200
MyVectorHandler * pXPrimePrime
Definition: invsolver.h:72
doublereal dIterTol
Definition: solver.h:306
integer iGetNumDofs(void) const
Definition: dataman.h:809
virtual integer GetIntegratorNumUnknownStates(void) const
Definition: stepsol.cc:87
virtual doublereal GetIntegratorDSolTol(void) const
Definition: stepsol.cc:105
#define SAFEDELETE(pnt)
Definition: mynewmem.h:710
LinSol CurrIntSolver
Definition: solver.h:317
virtual doublereal GetReal(const doublereal &dDefval=0.0)
Definition: parser.cc:1056