MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
TimeStepControl.cc
Go to the documentation of this file.
1 /*
2  * MBDyn (C) is a multibody analysis code.
3  * http://www.mbdyn.org
4  *
5  * Copyright (C) 1996-2017
6  *
7  * Pierangelo Masarati <masarati@aero.polimi.it>
8  * Paolo Mantegazza <mantegazza@aero.polimi.it>
9  *
10  * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano
11  * via La Masa, 34 - 20156 Milano, Italy
12  * http://www.aero.polimi.it
13  *
14  * Changing this copyright notice is forbidden.
15  *
16  * This program is free software; you can redistribute it and/or modify
17  * it under the terms of the GNU General Public License as published by
18  * the Free Software Foundation (version 2 of the License).
19  *
20  *
21  * This program is distributed in the hope that it will be useful,
22  * but WITHOUT ANY WARRANTY; without even the implied warranty of
23  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24  * GNU General Public License for more details.
25  *
26  * You should have received a copy of the GNU General Public License
27  * along with this program; if not, write to the Free Software
28  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
29  */
30 
31  /*
32  * With the contribution of Ankit Aggarwal <ankit.ankit.aggarwal@gmail.com>
33  * during Google Summer of Code 2016
34  */
35 
36 #include "mbconfig.h"
37 
38 #include "TimeStepControl.h"
39 #include "solver_impl.h"
40 
41 typedef std::map<std::string, TimeStepRead *> TimeStepFuncMapType;
42 
44 
46  bool IsWord(const std::string& s) const {
47  return TimeStepReadMap.find(s) != TimeStepReadMap.end();
48  };
49 };
50 
52 
53 bool
54 SetTimeStepData(const char *name, TimeStepRead * sr) {
55  return TimeStepReadMap.insert(std::make_pair(name, sr)).second;
56 }
57 
59  NO_OP;
60 }
61 
64 {
65  return dCurrTimeStep;
66 }
67 
68 void
69 NoChange::Init(integer iMaxIterations, doublereal dMinTimeStep, const DriveOwner& MaxTimeStep, doublereal dInitialTimeStep)
70 {
71  dCurrTimeStep = dInitialTimeStep;
72 
73  doublereal dInitialMaxTimeStep;
74  if (typeid(*MaxTimeStep.pGetDriveCaller()) == typeid(PostponedDriveCaller)) {
75  dInitialMaxTimeStep = std::numeric_limits<doublereal>::max();
76 
77  } else{
78  dInitialMaxTimeStep = MaxTimeStep.dGet();
79  }
80 
81  if (dMinTimeStep > dInitialMaxTimeStep) {
82  silent_cerr("warning: minimum time step " << dMinTimeStep << " greater than maximum (initial) time step " << dInitialMaxTimeStep << " (ignored)" << std::endl);
83  }
84 
85  if (dInitialTimeStep < dMinTimeStep) {
86  silent_cerr("warning: (constant) time step " << dInitialTimeStep << " less than minimum time step " << dMinTimeStep << " (ignored)" << std::endl);
87  }
88 
89  if (dInitialTimeStep > dInitialMaxTimeStep) {
90  silent_cerr("warning: (constant) time step " << dInitialTimeStep << " greater than maximum (initial) time step " << dInitialMaxTimeStep << " (ignored)" << std::endl);
91  }
92 }
93 
94 struct NoChangeTSR : public TimeStepRead {
95 public:
97 };
98 
101 {
102  // FIXME: the (constant) time step could be parsed here
103  return new NoChange(s);
104 }
105 
106 ChangeStep::ChangeStep(Solver *s, DriveCaller* pStrategyChangeDrive)
107 : s(s),
108 pStrategyChangeDrive(pStrategyChangeDrive),
109 dMinTimeStep(-1.)
110 {
111  NO_OP;
112 }
113 
116 {
117  doublereal dNewStep = pStrategyChangeDrive->dGet();
118  doublereal dMaxTimeStep = MaxTimeStep.dGet();;
119 
120  // are we sure we intend to change the time step if repeat is requested?
121  if (currStep == StepIntegrator::REPEATSTEP && dNewStep == dCurrTimeStep) {
122  return (dCurrTimeStep = dMinTimeStep/2.);
123  }
124 
125  return (dCurrTimeStep = std::max(std::min(dNewStep, dMaxTimeStep), dMinTimeStep));
126 }
127 
128 void
130 {
131  pStrategyChangeDrive->SetDrvHdl(driveHandler);
132 }
133 
134 void
135 ChangeStep::Init(integer iMaxIterations, doublereal dMinTimeStep, const DriveOwner& MaxTimeStep, doublereal dInitialTimeStep)
136 {
137  this->dCurrTimeStep = dInitialTimeStep;
138  this->MaxTimeStep.Set(MaxTimeStep.pGetDriveCaller()->pCopy());
139  this->dMinTimeStep = dMinTimeStep;
140 
141  doublereal dInitialMaxTimeStep ;
142  if (typeid(*MaxTimeStep.pGetDriveCaller()) == typeid(PostponedDriveCaller)) {
143  dInitialMaxTimeStep = std::numeric_limits<doublereal>::max();
144 
145  } else{
146  dInitialMaxTimeStep = MaxTimeStep.dGet();
147  }
148 
149  if (dMinTimeStep > dInitialMaxTimeStep) {
150  silent_cerr("error: minimum time step " << dMinTimeStep << " greater than maximum (initial) time step " << dInitialMaxTimeStep << std::endl);
152  }
153 
154  if (dInitialTimeStep < dMinTimeStep) {
155  silent_cerr("error: initial time step " << dInitialTimeStep << " less than minimum time step " << dMinTimeStep << std::endl);
157  }
158 
159  if (dInitialTimeStep > dInitialMaxTimeStep) {
160  silent_cerr("error: initial time step " << dInitialTimeStep << " greater than maximum (initial) time step " << dInitialMaxTimeStep << std::endl);
162  }
163 }
164 
165 class ChangeStepTSR : public TimeStepRead {
166 public:
168 };
169 
172 {
173  return new ChangeStep(s , HP.GetDriveCaller(true));
174 }
175 
177  doublereal dReductionFactor,
178  doublereal iStepsBeforeReduction,
179  doublereal dRaiseFactor,
180  doublereal iStepsBeforeRaise,
181  doublereal iMinIters,
182  doublereal iMaxIters)
183 : s(s),
184 dReductionFactor(dReductionFactor),
185 iStepsBeforeReduction(iStepsBeforeReduction),
186 dRaiseFactor(dRaiseFactor),
187 iStepsBeforeRaise(iStepsBeforeRaise),
188 iMinIters(iMinIters),
189 iMaxIters(iMaxIters),
190 bLastChance(false),
191 iStepsAfterReduction(0),
192 iStepsAfterRaise(0),
193 iWeightedPerformedIters(0),
194 dMinTimeStep(::dDefaultMinTimeStep)
195 {
196  NO_OP;
197 }
198 
201 {
202  doublereal dMaxTimeStep = MaxTimeStep.dGet();
203 
204  switch (Why) {
207  if (bLastChance) {
208  bLastChance = false;
209  }
211  return (dCurrTimeStep = std::min(dCurrTimeStep*dReductionFactor, dMaxTimeStep));
212  }
213 
214  if (!bLastChance) {
215  bLastChance = true;
216  return (dCurrTimeStep = dMinTimeStep);
217  }
218 
219  return (dCurrTimeStep = std::min(dCurrTimeStep*dReductionFactor, dMaxTimeStep));
220 
224 
225  iWeightedPerformedIters = (10*iPerformedIters + 9*iWeightedPerformedIters)/10;
226 
227  if (iPerformedIters > iMaxIters) {
229  bLastChance = false;
230  return (dCurrTimeStep = std::min(std::max(dCurrTimeStep*dReductionFactor, dMinTimeStep), dMaxTimeStep));
231 
232  }
233 
234  if (iPerformedIters <= iMinIters
237  && dCurrTimeStep < dMaxTimeStep)
238  {
239  iStepsAfterRaise = 0;
241  return (dCurrTimeStep = std::min(dCurrTimeStep*dRaiseFactor, dMaxTimeStep));
242  }
243 
244  return dCurrTimeStep;
245 
246  default:
247  // Should Not Reach Over here
248  ASSERT(0);
250  }
251 }
252 
253 
254 void
255 Factor::Init(integer iMaxIterations, doublereal dMinTimeStep, const DriveOwner& MaxTimeStep, doublereal dInitialTimeStep)
256 {
257  this->dCurrTimeStep = dInitialTimeStep;
258  this->MaxTimeStep.Set(MaxTimeStep.pGetDriveCaller()->pCopy());
259  this->dMinTimeStep = dMinTimeStep;
260 
261  if (iMaxIters <= iMinIters) {
262  silent_cerr("error: maximum number of iterations " << iMaxIters << " less than or equal to minimum " << iMinIters << std::endl);
264  }
265 
266  doublereal dInitialMaxTimeStep = MaxTimeStep.dGet();
267 
268  if (typeid(*MaxTimeStep.pGetDriveCaller()) == typeid(ConstDriveCaller)
269  && dInitialMaxTimeStep == std::numeric_limits<doublereal>::max())
270  {
271  silent_cerr("warning: maximum time step not set; the initial time step value " << dInitialTimeStep << " will be used" << std::endl);
272  }
273 
274  if (dMinTimeStep == ::dDefaultMinTimeStep) {
275  silent_cerr("warning: minimum time step not set; the initial time step value " << dInitialTimeStep << " will be used" << std::endl);
276  dMinTimeStep = dInitialTimeStep;
277  }
278 
279  if (dMinTimeStep == dInitialMaxTimeStep) {
280  silent_cerr("warning: minimum time step " << dMinTimeStep << "is equal to (initial) maximum time step " << dInitialMaxTimeStep << "; no time step adaptation will take place" << std::endl);
281 
282  } else if (dMinTimeStep == dInitialMaxTimeStep) {
283  silent_cerr("error: minimum time step " << dMinTimeStep << "is greater than (initial) maximum time step " << dInitialMaxTimeStep << std::endl);
285  }
286 }
287 
288 class FactorTSR : public TimeStepRead {
289 public:
291 };
292 
295 {
296  integer iMinIters = 1;
297  integer iMaxIters = 0;
298  doublereal dReductionFactor = 1.;
299  try {
300  dReductionFactor = HP.GetReal(1., HighParser::range_gt_le<doublereal>(0., 1.));
301 
303  silent_cerr("error: invalid reduction factor " << e.Get() << " (must be positive and less than (or equal to) one) [" << e.what() << "] at line " << HP.GetLineData() << std::endl);
304  throw e;
305  }
306 
307  integer iStepsBeforeReduction = 1;
308  try {
309  iStepsBeforeReduction = HP.GetInt(1, HighParser::range_ge<integer>(1));
310 
312  silent_cerr("error: invalid number of steps before reduction " << e.Get() << " (must be greater than zero) [" << e.what() << "] at line " << HP.GetLineData() << std::endl);
313  throw e;
314  }
315 
316  doublereal dRaiseFactor = 1.;
317  try {
318  dRaiseFactor = HP.GetReal(1., HighParser::range_ge<doublereal>(1.));
319 
321  silent_cerr("error: invalid increase factor " << e.Get() << " (must be greater than (or equal to) one) [" << e.what() << "] at line " << HP.GetLineData() << std::endl);
322  throw e;
323  }
324 
325  integer iStepsBeforeRaise = 1;
326  try {
327  iStepsBeforeRaise = HP.GetInt(1, HighParser::range_ge<integer>(1));
328 
330  silent_cerr("error: invalid number of steps before increase " << e.Get() << " (must be greater than zero) [" << e.what() << "] at line " << HP.GetLineData() << std::endl);
331  throw e;
332  }
333 
334  try {
335  iMinIters = HP.GetInt(1, HighParser::range_ge<integer>(1));
336 
338  silent_cerr("error: invalid number of minimum iterations " << e.Get() << " (must be greater than zero) [" << e.what() << "] at line " << HP.GetLineData() << std::endl);
339  throw e;
340  }
341 
342  iMaxIters = 0;
343  if (HP.IsArg()) {
344  try {
345  iMaxIters = HP.GetInt(1, HighParser::range_ge<integer>(1));
346 
348  silent_cerr("error: invalid number of maximum iterations " << e.Get() << " (must be greater than zero) [" << e.what() << "] at line " << HP.GetLineData() << std::endl);
349  throw e;
350  }
351  }
352 
353  // RETURN THE new FACTOR
354  return new Factor(s, dReductionFactor,
355  iStepsBeforeReduction,
356  dRaiseFactor,
357  iStepsBeforeRaise,
358  iMinIters,
359  iMaxIters);
360 }
361 
362 void
364 {
365  SetTimeStepData("no" "change", new NoChangeTSR);
366  SetTimeStepData("change", new ChangeStepTSR);
367  SetTimeStepData("factor", new FactorTSR);
368 }
369 
370 void
372 {
373  for (TimeStepFuncMapType::iterator it = TimeStepReadMap.begin(); it != TimeStepReadMap.end(); ++it) {
374  delete it->second;
375  }
376  TimeStepReadMap.clear();
377 }
378 
381 {
382  const char *str = HP.IsWord(TimeStepWordSet);
383  TimeStepFuncMapType::iterator func = TimeStepReadMap.find(std::string(str));
384  if (func == TimeStepReadMap.end()) {
385  silent_cerr("unknown TimeStepControl type " << s << "; using \"no change\"" << std::endl);
386  func = TimeStepReadMap.find(std::string("no" "change"));
387  }
388  ASSERT(func != TimeStepReadMap.end());
389  return func->second->Read(s, HP);
390 }
391 
void Init(integer iMaxIterations, doublereal dMinTimeStep, const DriveOwner &MaxTimeStep, doublereal dInitialTimeStep)
static const doublereal dDefaultMinTimeStep
Definition: solver_impl.h:100
doublereal dGetNewStepTime(StepIntegrator::StepChange Why, doublereal iPerformedIters)
void InitTimeStepData(void)
doublereal dReductionFactor
#define MBDYN_EXCEPT_ARGS
Definition: except.h:63
virtual integer GetInt(integer iDefval=0)
Definition: parser.cc:1050
void Init(integer iMaxIterations, doublereal dMinTimeStep, const DriveOwner &MaxTimeStep, doublereal dInitialTimeStep)
virtual const char * IsWord(const HighParser::WordSet &ws)
Definition: parser.cc:977
doublereal iStepsAfterRaise
DriveOwner MaxTimeStep
Factor(Solver *s, doublereal dReductionFactor, doublereal iStepsBeforeReduction, doublereal dRaiseFactor, doublereal iStepsBeforeRaise, doublereal iMinIters, doublereal iMaxIters)
const char * what(void) const
Definition: except.cc:54
ChangeStep(Solver *s, DriveCaller *pStrategyChangeDrive)
TimeStepControl * Read(Solver *s, MBDynParser &HP)
doublereal iStepsAfterReduction
TimeStepFuncMapType TimeStepReadMap
doublereal dMinTimeStep
#define NO_OP
Definition: myassert.h:74
TimeStepControl * Read(Solver *s, MBDynParser &HP)
DriveCaller * pStrategyChangeDrive
bool bLastChance
void DestroyTimeStepData(void)
void func(const T &u, const T &v, const T &w, doublereal e, T &f)
void Init(integer iMaxIterations, doublereal dMinTimeStep, const DriveOwner &MaxTimeStep, doublereal dInitialTimeStep)
doublereal iStepsBeforeReduction
virtual void SetDrvHdl(const DriveHandler *pDH)
Definition: drive.cc:487
NoChange(Solver *s)
Definition: mbdyn.h:76
doublereal iMaxIters
doublereal dGetNewStepTime(StepIntegrator::StepChange currStep, doublereal iPerformedIters)
doublereal dGetNewStepTime(StepIntegrator::StepChange currStep, doublereal iPerformedIters)
doublereal dMinTimeStep
#define ASSERT(expression)
Definition: colamd.c:977
DriveCaller * pGetDriveCaller(void) const
Definition: drive.cc:658
doublereal dCurrTimeStep
virtual DriveCaller * pCopy(void) const =0
Definition: solver.h:78
virtual doublereal dGet(const doublereal &dVar) const =0
doublereal iWeightedPerformedIters
virtual bool IsArg(void)
Definition: parser.cc:807
std::map< std::string, TimeStepRead * > TimeStepFuncMapType
void Set(const DriveCaller *pDC)
Definition: drive.cc:647
doublereal dGet(const doublereal &dVar) const
Definition: drive.cc:664
DriveOwner MaxTimeStep
bool SetTimeStepData(const char *name, TimeStepRead *sr)
void SetDriveHandler(const DriveHandler *driveHandler)
TimeStepControl * ReadTimeStepData(Solver *s, MBDynParser &HP)
DriveCaller * GetDriveCaller(bool bDeferred=false)
Definition: mbpar.cc:2033
doublereal iMinIters
doublereal dRaiseFactor
doublereal iStepsBeforeRaise
bool IsWord(const std::string &s) const
double doublereal
Definition: colamd.c:52
long int integer
Definition: colamd.c:51
virtual HighParser::ErrOut GetLineData(void) const
Definition: parsinc.cc:697
TimeStepWordSetType TimeStepWordSet
virtual doublereal GetReal(const doublereal &dDefval=0.0)
Definition: parser.cc:1056
TimeStepControl * Read(Solver *s, MBDynParser &HP)