LORENE
evolution.C
1/*
2 * Methods of template class Evolution
3 *
4 * (see file evolution.h for documentation).
5 *
6 */
7
8/*
9 * Copyright (c) 2004 Eric Gourgoulhon & Jerome Novak
10 *
11 * This file is part of LORENE.
12 *
13 * LORENE is free software; you can redistribute it and/or modify
14 * it under the terms of the GNU General Public License version 2
15 * as published by the Free Software Foundation.
16 *
17 * LORENE is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU General Public License for more details.
21 *
22 * You should have received a copy of the GNU General Public License
23 * along with LORENE; if not, write to the Free Software
24 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25 *
26 */
27
28/*
29 * $Id: evolution.C,v 1.18 2014/10/13 08:52:38 j_novak Exp $
30 * $Log: evolution.C,v $
31 * Revision 1.18 2014/10/13 08:52:38 j_novak
32 * Lorene classes and functions now belong to the namespace Lorene.
33 *
34 * Revision 1.17 2014/03/27 16:59:41 j_novak
35 * Added methods next_position(int) and previous_position(int). Changed (corrected + simplified) the interpolation method.
36 *
37 * Revision 1.16 2013/07/19 15:50:24 j_novak
38 * Implementation of the interpolation function for Evolution, with order=0, 1 or 2.
39 *
40 * Revision 1.15 2008/09/19 13:24:29 j_novak
41 * Added the third-order scheme for the time derivativre computation.
42 *
43 * Revision 1.14 2004/05/14 08:51:01 p_grandclement
44 * *** empty log message ***
45 *
46 * Revision 1.13 2004/05/14 08:41:05 e_gourgoulhon
47 * Added declaration "class Tbl ;" before the declarations of
48 * write_formatted.
49 *
50 * Revision 1.12 2004/05/13 21:30:32 e_gourgoulhon
51 * Use of function write_formatted in method save( ).
52 *
53 * Revision 1.11 2004/05/11 20:12:49 e_gourgoulhon
54 * Added methods j_min, j_max and save.
55 *
56 * Revision 1.10 2004/05/03 15:23:22 e_gourgoulhon
57 * Method downdate: changed the order of conditions (pos_jtop>=0)
58 * and (val[pos_jtop] == 0x0) in the while(...) test.
59 *
60 * Revision 1.9 2004/03/31 20:24:04 e_gourgoulhon
61 * Method time_derive: result object created by the copy constructor of
62 * class TyT, since the arithmetics may not return an object of exactly
63 * class TyT.
64 *
65 * Revision 1.8 2004/03/26 13:31:09 j_novak
66 * Definition of the macro UNDEF_STEP for non-defined time-steps.
67 * Changes in the way the time derivative is calculated.
68 *
69 * Revision 1.7 2004/03/26 08:22:13 e_gourgoulhon
70 * *** Full reorganization of class Evolution ***
71 * Introduction of the notion of absoluteuniversal time steps,
72 * stored in the new array 'step'.
73 * The new function position(int j) makes a correspondence
74 * between a universal time step j and the position in the
75 * arrays step, the_time and val.
76 * Only method update is now virtual.
77 * Methods operator[], position, is_known, downdate belong to
78 * the base class.
79 *
80 * Revision 1.6 2004/03/24 14:55:47 e_gourgoulhon
81 * Added method last_value().
82 *
83 * Revision 1.5 2004/03/23 14:50:41 e_gourgoulhon
84 * Added methods is_updated, downdate, get_jlast, get_size,
85 * as well as constructors without any initial value.
86 * Formatted documentation for Doxygen.
87 *
88 * Revision 1.4 2004/03/06 21:13:15 e_gourgoulhon
89 * Added time derivation (method time_derive).
90 *
91 * Revision 1.3 2004/02/17 22:13:34 e_gourgoulhon
92 * Suppressed declaration of global char[] evolution_C = ...
93 *
94 * Revision 1.2 2004/02/15 21:55:33 e_gourgoulhon
95 * Introduced derived classes Evolution_full and Evolution_std.
96 * Evolution is now an abstract base class.
97 *
98 * Revision 1.1 2004/02/13 15:53:20 e_gourgoulhon
99 * New (template) class for time evolution.
100 *
101 *
102 * $Header: /cvsroot/Lorene/C++/Include/Template/evolution.C,v 1.18 2014/10/13 08:52:38 j_novak Exp $
103 *
104 */
105
106// C++ headers
107#include "headcpp.h"
108
109// C headers
110#include <cstdlib>
111#include <cassert>
112#include <cmath>
113
114namespace Lorene {
115class Tbl ;
116
117void write_formatted(const double&, ostream& ) ;
118void write_formatted(const Tbl&, ostream& ) ;
119
120
121 //-------------------------//
122 // Constructors //
123 //-------------------------//
124
125
126template<typename TyT>
127Evolution<TyT>::Evolution(const TyT& initial_value, int initial_step,
128 double initial_time, int size_i)
129 : size(size_i),
130 pos_jtop(0) {
131
132 step = new int[size] ;
133 step[0] = initial_step ;
134 for (int j=1; j<size; j++) {
135 step[j] = UNDEF_STEP ;
136 }
137
138 the_time = new double[size] ;
139 the_time[0] = initial_time ;
140 for (int j=1; j<size; j++) {
141 the_time[j] = -1e20 ;
142 }
143
144 val = new TyT*[size] ;
145 val[0] = new TyT(initial_value) ;
146 for (int j=1; j<size; j++) {
147 val[j] = 0x0 ;
148 }
149
150}
152
153template<typename TyT>
155 : size(size_i),
156 pos_jtop(-1) {
157
158 step = new int[size] ;
159 for (int j=0; j<size; j++) {
160 step[j] = UNDEF_STEP ;
161 }
163 the_time = new double[size] ;
164 for (int j=0; j<size; j++) {
165 the_time[j] = -1e20 ;
166 }
167
168 val = new TyT*[size] ;
169 for (int j=0; j<size; j++) {
170 val[j] = 0x0 ;
171 }
172
173}
175
176
177template<typename TyT>
179 : size(evo.size),
180 pos_jtop(evo.pos_jtop) {
181
182 step = new int[size] ;
183 for (int j=0; j<size; j++) {
184 step[j] = evo.step[j] ;
185 }
187 the_time = new double[size] ;
188 for (int j=0; j<size; j++) {
189 the_time[j] = evo.the_time[j] ;
190 }
191
192 val = new TyT*[size] ;
193 for (int j=0; j<size; j++) {
194 if (evo.val[j] != 0x0) {
195 val[j] = new TyT( *(evo.val[j]) ) ;
197 else {
198 val[j] = 0x0 ;
199 }
200 }
201
203}
204
205
206
207
208 //-----------------------//
209 // Destructor //
210 //-----------------------//
212
213template<typename TyT>
215
216 delete [] step ;
217 delete [] the_time ;
219 for (int j=0; j<size; j++) {
220 if (val[j] != 0x0) delete val[j] ;
221 }
222
223 delete [] val ;
224
225}
226
227
228
229 //-----------------------//
230 // Mutators //
231 //-----------------------//
232
233
234template<typename TyT>
237 cerr << "void Evolution<TyT>::operator= : not implemented yet ! \n" ;
238 abort() ;
239
240}
241
242
243template<typename TyT>
245
246 if ( !(is_known(j)) ) return ; // a never updated step cannot
247 // be downdated
249 int pos = position(j) ;
250
251 assert( val[pos] != 0x0) ;
252
253 delete val[pos] ;
254 val[pos] = 0x0 ;
255 step[pos] = UNDEF_STEP ;
256 the_time[pos] = -1e20 ;
257
258 if (pos == pos_jtop) { // pos_jtop must be decreased
259 pos_jtop-- ;
260 while ( (pos_jtop>=0) && (val[pos_jtop] == 0x0) ) pos_jtop-- ;
261 }
262
263}
264
265
266
267
268 //------------//
269 // Accessors //
270 //------------//
271
272template<typename TyT>
273int Evolution<TyT>::position(int j) const {
274
275 assert(pos_jtop >= 0) ;
276 int jmax = step[pos_jtop] ;
277
278 if (j == jmax) return pos_jtop ; // for efficiency purpose
279
280 int pos = - 1 ;
281
282 if ( (j>=step[0]) && (j<jmax) ) {
283
284 for (int i=pos_jtop-1; i>=0; i--) { // cas i=pos_jtop treated above
285 if (step[i] == j) {
286 pos = i ;
287 break ;
288 }
289 }
290 }
291
292 if (pos == -1) {
293 cerr << "Evolution<TyT>::position: time step j = " <<
294 j << " not found !" << endl ;
295 abort() ;
296 }
297
298 return pos ;
299}
300
301template<typename TyT>
303
304 assert( (j>=0) && (j<=pos_jtop) ) ;
305 assert( step[j] != UNDEF_STEP ) ;
306
307 int n_pos = -1 ;
308
309 while ( (n_pos == -1) && ( j < pos_jtop ) ) {
310
311 j++ ;
312 if (step[j] != UNDEF_STEP) n_pos = j ;
313 }
314 return n_pos ;
315}
316
317template<typename TyT>
319
320 assert( (j>=0) && (j<=pos_jtop) ) ;
321 assert( step[j] != UNDEF_STEP ) ;
322
323 int n_pos = -1 ;
324
325 while ( (n_pos == -1) && ( j > 0 ) ) {
326
327 j-- ;
328 if (step[j] != UNDEF_STEP) n_pos = j ;
329 }
330 return n_pos ;
331}
332
333
334template<typename TyT>
335bool Evolution<TyT>::is_known(int j) const {
336
337 if (pos_jtop == -1) return false ;
338
339 assert(pos_jtop >= 0) ;
340
341 int jmax = step[pos_jtop] ;
342
343 if (j == jmax) {
344 return ( val[pos_jtop] != 0x0 ) ;
345 }
346
347 if ( (j>=step[0]) && (j<jmax) ) {
348
349 for (int i=pos_jtop-1; i>=0; i--) { // cas i=pos_jtop treated above
350
351 if (step[i] == j) return ( val[i] != 0x0 ) ;
352 }
353 }
354
355 return false ;
356}
357
358
359template<typename TyT>
360const TyT& Evolution<TyT>::operator[](int j) const {
361
362 TyT* pval = val[position(j)] ;
363 assert(pval != 0x0) ;
364 return *pval ;
365
366}
367
368
369template<typename TyT>
370TyT Evolution<TyT>::operator()(double t, int order) const {
371
372 int imin = position( j_min() ) ;
373
374 if ( ( t < the_time[ imin ] ) || ( t > the_time[pos_jtop] ) ) {
375 cerr << "Evolution<TyT>::operator()(double t, int order) : \n"
376 << "Requested time outside stored range!" << endl ;
377 abort() ;
378 }
379 assert ( order <= 2 ) ;
380 assert ( pos_jtop > order ) ;
381
382 int j0 = imin ;
383
384 while ((the_time[j0] < t) && ( j0 < pos_jtop )) {
385 j0 = next_position(j0) ;
386 assert( j0 != -1 ) ;
387 }
388
389 switch (order) {
390
391 case 0: {
392
393 return (*val[j0]) ;
394
395 break ;
396 }
397
398 case 1: {
399
400 int j1 = ( (j0 > imin) ? previous_position(j0) : next_position(j0) ) ;
401 assert( j1 != -1) ;
402
403 double x0 = the_time[j1] - t ;
404 double x1 = the_time[j0] - t ;
405 double dx = the_time[j0] - the_time[j1] ;
406
407 double a0 = x1 / dx ;
408 double a1 = x0 / dx ;
409
410 return (a0*(*val[j0]) - a1*(*val[j1])) ;
411
412 break ;
413 }
414
415 case 2: {
416
417 int j1 = ( (j0 == imin) ? next_position(j0) : previous_position(j0) ) ;
418 assert( j1 != -1) ;
419
420 int j2 = -1 ;
421 if (j0 == imin)
422 j2 = next_position(j1) ;
423 else {
424 if (j1 == imin)
425 j2 = next_position(j0) ;
426 else
427 j2 = previous_position(j1) ;
428 }
429 assert (j2 != -1 ) ;
430
431 double x0 = t - the_time[j0] ;
432 double x1 = t - the_time[j1] ;
433 double dx = the_time[j0] - the_time[j1] ;
434
435 double x2 = t - the_time[j2] ;
436
437 double dx0 = the_time[j1] - the_time[j2] ;
438 double dx1 = the_time[j0] - the_time[j2] ;
439 double dx2 = dx ;
440
441 double a0 = ( x2*x1 ) / ( dx2*dx1 ) ;
442 double a1 = ( x0*x2 ) / ( dx0*dx2 ) ;
443 double a2 = ( x0*x1 ) / ( dx0*dx1 ) ;
444
445 return ( a0*(*val[j0]) - a1*(*val[j1]) + a2*(*val[j2]) ) ;
446
447 break ;
448 }
449
450 default: {
451 cerr << " Evolution<TyT>::operator()(double t, int order) : \n" << endl ;
452 cerr << " The case order = " << order << " is not implemented!" << endl ;
453 abort() ;
454 break ;
455 }
456 }
457
458 return *val[j0] ;
459
460}
461
462template<typename TyT>
464
465 int resu = UNDEF_STEP ;
466 for (int i=0; i<=pos_jtop; i++) {
467 if (step[i] != UNDEF_STEP) {
468 resu = step[i] ;
469 break ;
470 }
471 }
472
473 if (resu == UNDEF_STEP) {
474 cerr << "Evolution<TyT>::j_min() : no valid time step found !" << endl ;
475 abort() ;
476 }
477
478 return resu ;
479}
480
481template<typename TyT>
483
484 if (pos_jtop == -1) {
485 cerr << "Evolution<TyT>::j_max() : no valid time step found !" << endl ;
486 abort() ;
487 }
488
489 assert(pos_jtop >=0) ;
490 int jmax = step[pos_jtop] ;
491 assert(jmax != UNDEF_STEP) ;
492 return jmax ;
493}
494
495
496
497
498
499 //-----------------------//
500 // Time derivative //
501 //-----------------------//
502
503template<typename TyT>
504TyT Evolution<TyT>::time_derive(int j, int n) const {
505
506 if (n == 0) {
507 TyT resu ( operator[](j) ) ;
508 resu = 0 * resu ;
509
510 return resu ;
511
512 }
513
514 else {
515
516 int pos = position(j) ;
517 assert ( pos > 0 ) ;
518 assert ( step[pos-1] != UNDEF_STEP ) ;
519
520 switch (n) {
521
522 case 1 : {
523
524 double dt = the_time[pos] - the_time[pos-1] ;
525
526 return ( (*val[pos]) - (*val[pos-1]) ) / dt ;
527 break ;
528 }
529
530 case 2 : {
531
532 assert ( pos > 1 ) ;
533 assert ( step[pos-2] != UNDEF_STEP ) ;
534 double dt = the_time[pos] - the_time[pos-1] ;
535#ifndef NDEBUG
536 double dt2 = the_time[pos-1] - the_time[pos-2] ;
537 if (fabs(dt2 -dt) > 1.e-13) {
538 cerr <<
539 "Evolution<TyT>::time_derive: the current version is"
540 << " valid only for \n"
541 << " a constant time step !" << endl ;
542 abort() ;
543 }
544#endif
545 return ( 1.5* (*val[pos]) - 2.* (*val[pos-1]) + 0.5* (*val[pos-2]) ) / dt ;
546 break ;
547 }
548
549 case 3 : {
550
551 assert ( pos > 2 ) ;
552 assert ( step[pos-2] != UNDEF_STEP ) ;
553 assert ( step[pos-3] != UNDEF_STEP ) ;
554 double dt = the_time[pos] - the_time[pos-1] ;
555#ifndef NDEBUG
556 double dt2 = the_time[pos-1] - the_time[pos-2] ;
557 double dt3 = the_time[pos-2] - the_time[pos-3] ;
558 if ((fabs(dt2 -dt) > 1.e-13)||(fabs(dt3 -dt2) > 1.e-13)) {
559 cerr <<
560 "Evolution<TyT>::time_derive: the current version is valid only for \n"
561 << " a constant time step !" << endl ;
562 abort() ;
563 }
564#endif
565 return ( 11.*(*val[pos]) - 18.*(*val[pos-1]) + 9.*(*val[pos-2])
566 - 2.*(*val[pos-3])) / (6.*dt) ;
567 break ;
568 }
569 default : {
570 cerr << "Evolution<TyT>::time_derive: the case n = " << n
571 << " is not implemented !" << endl ;
572 abort() ;
573 break ;
574 }
575 }
576
577 }
578 return operator[](j) ;
579}
580
581
582
583 //-----------------------//
584 // Outputs //
585 //-----------------------//
586
587
588template<typename TyT>
589void Evolution<TyT>::save(const char* filename) const {
590
591 ofstream fich(filename) ;
592
593 time_t temps = time(0x0) ;
594
595 fich << "# " << filename << " " << ctime(&temps) ;
596 fich << "# " << size << " size" << '\n' ;
597 fich << "# " << pos_jtop << " pos_jtop" << '\n' ;
598 fich << "# t value... \n" ;
599
600 fich.precision(14) ;
601 fich.setf(ios::scientific) ;
602 fich.width(20) ;
603 for (int i=0; i<=pos_jtop; i++) {
604 if (step[i] != UNDEF_STEP) {
605 fich << the_time[i] ; fich.width(23) ;
606 assert(val[i] != 0x0) ;
607 write_formatted(*(val[i]), fich) ;
608 fich << '\n' ;
609 }
610 }
611
612 fich.close() ;
613}
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632}
Time evolution (*** under development ***).
Definition evolution.h:120
void downdate(int j)
Suppresses a stored value.
Definition evolution.C:244
void save(const char *filename) const
Saves *this in a formatted file.
Definition evolution.C:589
Evolution(const TyT &initial_value, int initial_j, double initial_time, int initial_size)
Constructor from some initial value.
Definition evolution.C:127
const TyT & operator[](int j) const
Returns the value at time step j.
Definition evolution.C:360
int j_max() const
Returns the larger time step j stored in *this.
Definition evolution.C:482
int next_position(int i) const
Returns the next valid position (returns -1 if none is found)
Definition evolution.C:302
TyT time_derive(int j, int n=2) const
Computes the time derivative at time step j by means of a n-th order scheme, from the values at steps...
Definition evolution.C:504
TyT ** val
Array of pointers onto the values (size = size).
Definition evolution.h:137
TyT operator()(double t, int order=2) const
Returns the value at time t, with a scheme of order order.
Definition evolution.C:370
int size
Maximum number of stored time steps.
Definition evolution.h:128
virtual ~Evolution()
Destructor.
Definition evolution.C:214
int * step
Array of time step indices (size = size).
Definition evolution.h:131
int j_min() const
Returns the smaller time step j stored in *this.
Definition evolution.C:463
double * the_time
Array of values of t at the various time steps (size = size).
Definition evolution.h:134
int position(int j) const
Gives the position in the arrays step, the_time and val corresponding to the time step j.
Definition evolution.C:273
bool is_known(int j) const
Checks whether the value a given time step has been set.
Definition evolution.C:335
virtual void operator=(const Evolution< TyT > &t_in)
Assignement.
Definition evolution.C:235
int previous_position(int i) const
Returns the previous valid position (returns -1 if none is found)
Definition evolution.C:318
Lorene prototypes.
Definition app_hor.h:64