1 /* Siconos is a program dedicated to modeling, simulation and control 2 * of non smooth dynamical systems. 3 * 4 * Copyright 2021 INRIA. 5 * 6 * Licensed under the Apache License, Version 2.0 (the "License"); 7 * you may not use this file except in compliance with the License. 8 * You may obtain a copy of the License at 9 * 10 * http://www.apache.org/licenses/LICENSE-2.0 11 * 12 * Unless required by applicable law or agreed to in writing, software 13 * distributed under the License is distributed on an "AS IS" BASIS, 14 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 15 * See the License for the specific language governing permissions and 16 * limitations under the License. 17 */ 18 /*! \file TimeStepping.hpp 19 * \brief Time-Stepping simulation 20 */ 21 #ifndef TimeStepping_H 22 #define TimeStepping_H 23 24 #include "Simulation.hpp" 25 26 /** type of function used to post-treat output info from solver. */ 27 typedef void (*CheckSolverFPtr)(int, Simulation*); 28 29 /** \brief Event-capturing Time-Stepping simulation 30 * 31 * This class implements the basic algorithm for Event-capturing Time-Stepping 32 * simulations. 33 * 34 * References : 35 * 36 * V. Acary and B. Brogliato. Numerical Methods for Nonsmooth Dynamical Systems: 37 * Applications in Mechanics and Electronics, volume 35 of Lecture Notes in 38 * Applied and Computational Mechanics. Springer Verlag, 2008. 39 * 40 */ 41 42 #define SICONOS_TS_LINEAR 1 43 #define SICONOS_TS_LINEAR_IMPLICIT 2 44 #define SICONOS_TS_NONLINEAR 3 45 46 class TimeStepping : public Simulation 47 { 48 protected: 49 /** serialization hooks 50 */ 51 ACCEPT_SERIALIZATION(TimeStepping); 52 53 /** Default Newton tolerance used in call of run() of ComputeOneStep() */ 54 double _newtonTolerance; 55 56 /** Default maximum number of Newton iteration*/ 57 unsigned int _newtonMaxIteration; 58 59 /** Number of steps performed in the Newton Loop */ 60 unsigned int _newtonNbIterations; 61 62 /** Cumulative number of steps performed in the Newton Loops */ 63 unsigned int _newtonCumulativeNbIterations; 64 65 /** unsigned int _newtonOptions 66 * option in the Newon iteration 67 * SICONOS_TS_LINEAR or SICONOS_TS_LINEAR_IMPLICIT SICONOS_TS_NONLINEAR will force a single iteration of the Newton Solver 68 * SICONOS_TS_NONLINEAR (default) will perform the newton iteration up to convergence 69 */ 70 unsigned int _newtonOptions; 71 72 73 /** Maximum Residual for the Dynamical system */ 74 double _newtonResiduDSMax; 75 76 /** Maximum Residual for the output of the relation */ 77 double _newtonResiduYMax; 78 79 /** Maximum Residual for the input of the relation */ 80 double _newtonResiduRMax; 81 82 /** boolean variable to know whether the ResiduY has to be computed or not 83 * if true, the ResiduY is computed and the convergence is checked 84 */ 85 bool _computeResiduY; 86 87 /** boolean variable to know whether the ResiduR has to be computed or not 88 * if true, the ResiduR is computed and the convergence is checked 89 */ 90 bool _computeResiduR; 91 92 /** boolean variable to know whether Newton iterations converge or not 93 */ 94 bool _isNewtonConverge; 95 96 /** boolean variable indicating whether interactions should be 97 * updated within the Newton loop. 98 */ 99 bool _newtonUpdateInteractionsPerIteration; 100 101 /** boolean variable to display Newton info 102 */ 103 bool _displayNewtonConvergence; 104 105 /** boolean variable to display warning on non-convergence 106 */ 107 bool _warnOnNonConvergence; 108 109 /** boolean variable to resetAllLamda at each step (default true) 110 */ 111 bool _resetAllLambda; 112 113 /** Default Constructor 114 */ TimeStepping()115 TimeStepping() : 116 _computeResiduY(false), 117 _computeResiduR(false), 118 _isNewtonConverge(false) {}; 119 120 121 /** newton algorithm 122 * \param criterion convergence criterion 123 * \param maxStep maximum number of Newton steps 124 */ 125 virtual void newtonSolve(double criterion, unsigned int maxStep); 126 127 128 public: 129 130 /** initialisation specific to TimeStepping for OneStepNSProblem. 131 */ 132 virtual void initOSNS(); 133 134 /** Standard constructor 135 * \param nsds NonSmoothDynamicalSystem to be simulated 136 * \param td pointer to a timeDiscretisation used in the integration 137 * \param osi one step integrator (default none) 138 * \param osnspb one step non smooth problem (default none) 139 */ 140 TimeStepping(SP::NonSmoothDynamicalSystem nsds, SP::TimeDiscretisation td, 141 SP::OneStepIntegrator osi, 142 SP::OneStepNSProblem osnspb); 143 144 /** Constructor with the time-discretisation. 145 * \param nsds NonSmoothDynamicalSystem to be simulated 146 * \param td pointer to a timeDiscretisation used in the integration 147 * \param nb number of non smooth problem 148 */ 149 TimeStepping(SP::NonSmoothDynamicalSystem nsds, SP::TimeDiscretisation td, int nb =0); 150 151 /** insert an Integrator into the simulation list of integrators 152 * \param osi the OneStepIntegrator to add 153 */ 154 virtual void insertIntegrator(SP::OneStepIntegrator osi); 155 156 /** Destructor. 157 */ 158 virtual ~TimeStepping(); 159 160 /** update indexSets[i] of the topology, using current y and lambda values of Interactions. 161 * \param i the number of the set to be updated 162 */ 163 virtual void updateIndexSet(unsigned int i); 164 165 // /** Used by the updateIndexSet function in order to deactivate SP::Interaction. 166 // */ 167 // virtual bool predictorDeactivate(SP::Interaction inter, unsigned int i); 168 169 // /** Used by the updateIndexSet function in order to activate SP::Interaction. 170 // */ 171 // virtual bool predictorActivate(SP::Interaction inter, unsigned int i); 172 173 /** increment model current time according to User TimeDiscretisation and call SaveInMemory. */ 174 virtual void nextStep(); 175 176 /** integrates all the DynamicalSystems taking not into account nslaw, reactions (ie non-smooth part) ... 177 */ 178 void computeFreeState(); 179 180 /** Reset all lambdas of all interactions */ 181 void resetLambdas(); 182 183 /** step from current event to next event of EventsManager 184 */ 185 void advanceToEvent(); 186 187 /** run one time--step of the simulation 188 */ 189 void computeOneStep(); 190 191 192 193 /** To known the number of steps performed by the Newton algorithm. 194 * \return the number of steps performed by the Newton algorithm 195 */ getNewtonNbIterations()196 unsigned int getNewtonNbIterations() 197 { 198 return _newtonNbIterations; 199 } 200 201 /** To known the number of steps performed by the Newton algorithm. 202 * \return the cumulative number of steps performed by the Newton algorithm 203 */ getNewtonCumulativeNbIterations()204 unsigned int getNewtonCumulativeNbIterations() 205 { 206 return _newtonCumulativeNbIterations; 207 } 208 209 /** initialize the Newton 210 * It computes the initial residu and set the, if needed to Newton variable 211 * to start the newton algorithm. 212 */ 213 void initializeNewtonLoop(); 214 215 216 void prepareNewtonIteration(); 217 218 /** check the convergence of Newton algorithm according to criterion 219 * \param criterion convergence criterion 220 * \return bool = true if Newton method has converged 221 */ 222 bool newtonCheckConvergence(double criterion); 223 224 /*save y_k^p, the current Newton iteration*/ 225 void saveYandLambdaInOldVariables(); 226 227 /** run the simulation, from t0 to T 228 * with default parameters if any setting has been done 229 */ 230 void run(); 231 232 /** check returning value from computeOneStepNSProblem and process 233 * \param info solver-specific error code return by the nonsmooth solver 234 */ 235 void DefaultCheckSolverOutput(int info); 236 237 /** Set CheckSolverOutput function 238 * \param newF pointer to function steering the behavior of simulation when 239 * nonsmooth solver failed 240 */ 241 void setCheckSolverFunction(CheckSolverFPtr newF); 242 243 /** */ isNewtonConverge()244 bool isNewtonConverge() 245 { 246 return _isNewtonConverge; 247 }; 248 displayNewtonConvergence()249 bool displayNewtonConvergence() 250 { 251 return _displayNewtonConvergence; 252 }; setDisplayNewtonConvergence(bool newval)253 void setDisplayNewtonConvergence(bool newval) 254 { 255 _displayNewtonConvergence = newval; 256 }; 257 setWarnOnNonConvergence(bool newval)258 void setWarnOnNonConvergence(bool newval) 259 { 260 _warnOnNonConvergence = newval; 261 }; warnOnNonConvergence()262 bool warnOnNonConvergence() 263 { 264 return _warnOnNonConvergence; 265 }; 266 void displayNewtonConvergenceAtTheEnd(int info, unsigned int maxStep); 267 268 void displayNewtonConvergenceInTheLoop(); 269 setResetAllLambda(bool newval)270 void setResetAllLambda(bool newval) 271 { 272 _resetAllLambda = newval; 273 }; 274 275 276 /** To specify if the output interaction residu must be computed. 277 * \param v set to true when the output interaction residu must be computed 278 */ setComputeResiduY(bool v)279 void setComputeResiduY(bool v) 280 { 281 _computeResiduY = v; 282 }; 283 284 /** To know if the output interaction residu must be computed. 285 * \return bool _computeResiduY 286 */ computeResiduY()287 virtual bool computeResiduY() 288 { 289 return _computeResiduY; 290 }; 291 292 293 /** To specify if the input interaction residu must be computed. 294 * \param v set to true when the input interaction residu must be computed 295 */ setComputeResiduR(bool v)296 void setComputeResiduR(bool v) 297 { 298 _computeResiduR = v; 299 }; 300 301 /** To known if the input interaction residu must be computed. 302 * \return bool _computeResiduR 303 */ computeResiduR()304 virtual bool computeResiduR() 305 { 306 return _computeResiduR; 307 }; 308 309 310 /** set the Default Newton tolerance 311 * \param tol Newton solver tolerance 312 */ setNewtonTolerance(double tol)313 void setNewtonTolerance(double tol) 314 { 315 _newtonTolerance = tol; 316 }; 317 318 /** get the Newton tolerance 319 * \return default Newton solver tolerance 320 */ newtonTolerance()321 double newtonTolerance() 322 { 323 return _newtonTolerance; 324 }; 325 326 /** set the maximum number of Newton iteration 327 * \param maxStep maximum number of Newton solver iterations 328 */ setNewtonMaxIteration(unsigned int maxStep)329 void setNewtonMaxIteration(unsigned int maxStep) 330 { 331 _newtonMaxIteration = maxStep; 332 }; 333 334 /** get the maximum number of Newton iteration 335 * \return maximum number of Newton solver iterations 336 */ newtonMaxIteration()337 unsigned int newtonMaxIteration() 338 { 339 340 return _newtonMaxIteration; 341 }; 342 343 /** set the NewtonOptions 344 * \param v Newton solver options 345 */ setNewtonOptions(unsigned int v)346 void setNewtonOptions(unsigned int v) 347 { 348 _newtonOptions = v; 349 }; 350 351 /** get the NewtonOptions 352 * \return Newton solver options - SICONOS_TS_LINEAR 1, 353 * SICONOS_TS_LINEAR_IMPLICIT 2, SICONOS_TS_NONLINEAR 3 354 */ newtonOptions()355 unsigned int newtonOptions() 356 { 357 return _newtonOptions; 358 }; 359 360 361 /** accessor to _newtonResiduDSMax 362 * \return _newtonResiduDSMax 363 */ newtonResiduDSMax()364 double newtonResiduDSMax() 365 { 366 return _newtonResiduDSMax; 367 }; 368 369 /** accessor to _newtonResiduYMax 370 * \return _newtonResiduYMax 371 */ newtonResiduYMax()372 double newtonResiduYMax() 373 { 374 return _newtonResiduYMax; 375 }; 376 377 /** accessor to _newtonResiduRMax 378 * \return _newtonResiduRMax 379 */ newtonResiduRMax()380 double newtonResiduRMax() 381 { 382 return _newtonResiduRMax; 383 }; 384 385 /** visitors hook 386 */ 387 ACCEPT_STD_VISITORS(); 388 389 }; 390 391 #endif // TimeStepping_H 392