00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 #ifndef __DIFFERENTIALSTEPPER_HPP
00032 #define __DIFFERENTIALSTEPPER_HPP
00033
00034 #include "libecs.hpp"
00035 #include "Stepper.hpp"
00036
00037 #include "boost/multi_array.hpp"
00038
00039 namespace libecs
00040 {
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056 typedef boost::multi_array<Real, 2> RealMatrix_;
00057 DECLARE_TYPE( RealMatrix_, RealMatrix );
00058
00059
00060 DECLARE_CLASS( DifferentialStepper );
00061
00062 LIBECS_DM_CLASS( DifferentialStepper, Stepper )
00063 {
00064
00065 public:
00066
00067 LIBECS_DM_OBJECT_ABSTRACT( DifferentialStepper )
00068 {
00069 INHERIT_PROPERTIES( Stepper );
00070
00071
00072 PROPERTYSLOT( Real, StepInterval,
00073 &DifferentialStepper::initializeStepInterval,
00074 &DifferentialStepper::getStepInterval );
00075
00076 PROPERTYSLOT_GET_NO_LOAD_SAVE( Real, NextStepInterval );
00077 PROPERTYSLOT_SET_GET_NO_LOAD_SAVE( Real, TolerableStepInterval );
00078 PROPERTYSLOT_GET_NO_LOAD_SAVE( Integer, Stage );
00079 PROPERTYSLOT_GET_NO_LOAD_SAVE( Integer, Order );
00080 }
00081
00082 class Interpolant
00083 :
00084 public libecs::Interpolant
00085 {
00086
00087 public:
00088
00089 Interpolant( DifferentialStepperRef aStepper,
00090 VariablePtr const aVariablePtr )
00091 :
00092 libecs::Interpolant( aVariablePtr ),
00093 theStepper( aStepper ),
00094 theIndex( theStepper.getVariableIndex( aVariablePtr ) )
00095 {
00096 ;
00097 }
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149 virtual const Real getDifference( RealParam aTime,
00150 RealParam anInterval ) const
00151 {
00152 if ( !theStepper.theStateFlag )
00153 {
00154 return 0.0;
00155 }
00156
00157 const Real aTimeInterval1( aTime - theStepper.getCurrentTime() );
00158 const Real aTimeInterval2( aTimeInterval1 - anInterval );
00159
00160 RealMatrixCref aTaylorSeries( theStepper.getTaylorSeries() );
00161 RealCptr aTaylorCoefficientPtr( aTaylorSeries.origin() + theIndex );
00162
00163
00164
00165
00166
00167 Real aValue1( *aTaylorCoefficientPtr * aTimeInterval1 );
00168 Real aValue2( *aTaylorCoefficientPtr * aTimeInterval2 );
00169
00170
00171
00172
00173 const RealMatrix::size_type aTaylorSize( theStepper.getOrder() );
00174 if( aTaylorSize >= 2)
00175 {
00176 const Real
00177 aStepIntervalInv( 1.0 / theStepper.getTolerableStepInterval() );
00178
00179 const RealMatrix::size_type aStride( aTaylorSeries.strides()[0] );
00180
00181 Real aFactorialInv1( aTimeInterval1 );
00182 Real aFactorialInv2( aTimeInterval2 );
00183
00184 RealMatrix::size_type s( aTaylorSize - 1 );
00185
00186 const Real theta1( aTimeInterval1 * aStepIntervalInv );
00187 const Real theta2( aTimeInterval2 * aStepIntervalInv );
00188
00189 do
00190 {
00191
00192
00193
00194 aTaylorCoefficientPtr += aStride;
00195 const Real aTaylorCoefficient( *aTaylorCoefficientPtr );
00196
00197 aFactorialInv1 *= theta1;
00198 aFactorialInv2 *= theta2;
00199
00200 aValue1 += aTaylorCoefficient * aFactorialInv1;
00201 aValue2 += aTaylorCoefficient * aFactorialInv2;
00202
00203
00204 --s;
00205 } while( s != 0 );
00206 }
00207
00208 return aValue1 - aValue2;
00209 }
00210
00211 virtual const Real getVelocity( RealParam aTime ) const
00212 {
00213 if ( !theStepper.theStateFlag )
00214 {
00215 return 0.0;
00216 }
00217
00218 const Real aTimeInterval( aTime - theStepper.getCurrentTime() );
00219
00220 RealMatrixCref aTaylorSeries( theStepper.getTaylorSeries() );
00221 RealCptr aTaylorCoefficientPtr( aTaylorSeries.origin() + theIndex );
00222
00223
00224
00225
00226
00227 Real aValue( *aTaylorCoefficientPtr );
00228
00229
00230
00231
00232 const RealMatrix::size_type aTaylorSize( theStepper.getStage() );
00233 if( aTaylorSize >= 2 && aTimeInterval != 0.0 )
00234 {
00235 const RealMatrix::size_type aStride( aTaylorSeries.strides()[0] );
00236
00237 Real aFactorialInv( 1.0 );
00238
00239 RealMatrix::size_type s( 1 );
00240
00241 const Real theta( aTimeInterval
00242 / theStepper.getTolerableStepInterval() );
00243
00244 do
00245 {
00246
00247 ++s;
00248
00249 aTaylorCoefficientPtr += aStride;
00250 const Real aTaylorCoefficient( *aTaylorCoefficientPtr );
00251
00252 aFactorialInv *= theta * s;
00253
00254 aValue += aTaylorCoefficient * aFactorialInv;
00255
00256
00257 } while( s != aTaylorSize );
00258 }
00259
00260 return aValue;
00261 }
00262
00263 protected:
00264
00265 DifferentialStepperRef theStepper;
00266 VariableVector::size_type theIndex;
00267
00268 };
00269
00270 public:
00271
00272 ECELL_API DifferentialStepper();
00273
00274 ECELL_API virtual ~DifferentialStepper();
00275
00276 SET_METHOD( Real, NextStepInterval )
00277 {
00278 theNextStepInterval = value;
00279 }
00280
00281 GET_METHOD( Real, NextStepInterval )
00282 {
00283 return theNextStepInterval;
00284 }
00285
00286 SET_METHOD( Real, TolerableStepInterval )
00287 {
00288 theTolerableStepInterval = value;
00289 }
00290
00291 GET_METHOD( Real, TolerableStepInterval )
00292 {
00293 return theTolerableStepInterval;
00294 }
00295
00296 void initializeStepInterval( RealParam aStepInterval )
00297 {
00298 setStepInterval( aStepInterval );
00299 setTolerableStepInterval( aStepInterval );
00300 setNextStepInterval( aStepInterval );
00301 }
00302
00303 ECELL_API void resetAll();
00304
00305 ECELL_API void interIntegrate();
00306
00307 void initializeVariableReferenceList();
00308
00309 ECELL_API void setVariableVelocity( boost::detail::multi_array::sub_array<Real, 1> aVelocityBuffer );
00310
00311 ECELL_API virtual void initialize();
00312
00313 ECELL_API virtual void reset();
00314
00315 ECELL_API virtual void interrupt( TimeParam aTime );
00316
00317 virtual InterpolantPtr createInterpolant( VariablePtr aVariable )
00318 {
00319 return new DifferentialStepper::Interpolant( *this, aVariable );
00320 }
00321
00322 virtual GET_METHOD( Integer, Stage )
00323 {
00324 return 1;
00325 }
00326
00327 virtual GET_METHOD( Integer, Order )
00328 {
00329 return getStage();
00330 }
00331
00332 RealMatrixCref getTaylorSeries() const
00333 {
00334 return theTaylorSeries;
00335 }
00336
00337 protected:
00338
00339 const bool isExternalErrorTolerable() const;
00340
00341 RealMatrix theTaylorSeries;
00342
00343 std::vector< std::vector<Integer> > theVariableReferenceListVector;
00344
00345 bool theStateFlag;
00346
00347 private:
00348
00349 Real theNextStepInterval;
00350 Real theTolerableStepInterval;
00351 };
00352
00353
00354
00355
00356
00357
00358
00359
00360 DECLARE_CLASS( AdaptiveDifferentialStepper );
00361
00362 LIBECS_DM_CLASS( AdaptiveDifferentialStepper, DifferentialStepper )
00363 {
00364
00365 public:
00366
00367 LIBECS_DM_OBJECT_ABSTRACT( AdaptiveDifferentialStepper )
00368 {
00369 INHERIT_PROPERTIES( DifferentialStepper );
00370
00371 PROPERTYSLOT_SET_GET( Real, Tolerance );
00372 PROPERTYSLOT_SET_GET( Real, AbsoluteToleranceFactor );
00373 PROPERTYSLOT_SET_GET( Real, StateToleranceFactor );
00374 PROPERTYSLOT_SET_GET( Real, DerivativeToleranceFactor );
00375
00376 PROPERTYSLOT( Integer, IsEpsilonChecked,
00377 &AdaptiveDifferentialStepper::setEpsilonChecked,
00378 &AdaptiveDifferentialStepper::isEpsilonChecked );
00379 PROPERTYSLOT_SET_GET( Real, AbsoluteEpsilon );
00380 PROPERTYSLOT_SET_GET( Real, RelativeEpsilon );
00381
00382 PROPERTYSLOT_GET_NO_LOAD_SAVE( Real, MaxErrorRatio );
00383 }
00384
00385 public:
00386
00387 ECELL_API AdaptiveDifferentialStepper();
00388
00389 ECELL_API virtual ~AdaptiveDifferentialStepper();
00390
00391
00392
00393
00394
00395
00396
00397 SET_METHOD( Real, Tolerance )
00398 {
00399 theTolerance = value;
00400 }
00401
00402 GET_METHOD( Real, Tolerance )
00403 {
00404 return theTolerance;
00405 }
00406
00407 SET_METHOD( Real, AbsoluteToleranceFactor )
00408 {
00409 theAbsoluteToleranceFactor = value;
00410 }
00411
00412 GET_METHOD( Real, AbsoluteToleranceFactor )
00413 {
00414 return theAbsoluteToleranceFactor;
00415 }
00416
00417 SET_METHOD( Real, StateToleranceFactor )
00418 {
00419 theStateToleranceFactor = value;
00420 }
00421
00422 GET_METHOD( Real, StateToleranceFactor )
00423 {
00424 return theStateToleranceFactor;
00425 }
00426
00427 SET_METHOD( Real, DerivativeToleranceFactor )
00428 {
00429 theDerivativeToleranceFactor = value;
00430 }
00431
00432 GET_METHOD( Real, DerivativeToleranceFactor )
00433 {
00434 return theDerivativeToleranceFactor;
00435 }
00436
00437 SET_METHOD( Real, MaxErrorRatio )
00438 {
00439 theMaxErrorRatio = value;
00440 }
00441
00442 GET_METHOD( Real, MaxErrorRatio )
00443 {
00444 return theMaxErrorRatio;
00445 }
00446
00447
00448
00449
00450
00451 SET_METHOD( Integer, EpsilonChecked )
00452 {
00453 if ( value > 0 ) {
00454 theEpsilonChecked = true;
00455 }
00456 else {
00457 theEpsilonChecked = false;
00458 }
00459 }
00460
00461 const Integer isEpsilonChecked() const
00462 {
00463 return theEpsilonChecked;
00464 }
00465
00466 SET_METHOD( Real, AbsoluteEpsilon )
00467 {
00468 theAbsoluteEpsilon = value;
00469 }
00470
00471 GET_METHOD( Real, AbsoluteEpsilon )
00472 {
00473 return theAbsoluteEpsilon;
00474 }
00475
00476 SET_METHOD( Real, RelativeEpsilon )
00477 {
00478 theRelativeEpsilon = value;
00479 }
00480
00481 GET_METHOD( Real, RelativeEpsilon )
00482 {
00483 return theRelativeEpsilon;
00484 }
00485
00486 ECELL_API virtual void initialize();
00487
00488 ECELL_API virtual void step();
00489
00490 virtual bool calculate() = 0;
00491
00492 virtual GET_METHOD( Integer, Stage )
00493 {
00494 return 2;
00495 }
00496
00497 private:
00498
00499 Real safety;
00500 Real theTolerance;
00501 Real theAbsoluteToleranceFactor;
00502 Real theStateToleranceFactor;
00503 Real theDerivativeToleranceFactor;
00504
00505 bool theEpsilonChecked;
00506 Real theAbsoluteEpsilon;
00507 Real theRelativeEpsilon;
00508
00509 Real theMaxErrorRatio;
00510 };
00511
00512
00513 }
00514
00515 #endif
00516
00517
00518
00519
00520
00521
00522
00523
00524