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 
114 namespace Lorene {
115 class Tbl ;
116 
117 void write_formatted(const double&, ostream& ) ;
118 void write_formatted(const Tbl&, ostream& ) ;
119 
120 
121  //-------------------------//
122  // Constructors //
123  //-------------------------//
124 
125 
126 template<typename TyT>
127 Evolution<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 }
151 
152 
153 template<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  }
162 
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 }
174 
175 
176 
177 template<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  }
186 
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]) ) ;
196  }
197  else {
198  val[j] = 0x0 ;
199  }
200  }
201 
202 
203 }
204 
205 
206 
207 
208  //-----------------------//
209  // Destructor //
210  //-----------------------//
211 
212 
213 template<typename TyT>
215 
216  delete [] step ;
217  delete [] the_time ;
218 
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 
234 template<typename TyT>
236 
237  cerr << "void Evolution<TyT>::operator= : not implemented yet ! \n" ;
238  abort() ;
239 
240 }
241 
242 
243 template<typename TyT>
245 
246  if ( !(is_known(j)) ) return ; // a never updated step cannot
247  // be downdated
248 
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 
272 template<typename TyT>
273 int 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 
301 template<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 
317 template<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 
334 template<typename TyT>
335 bool 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 
359 template<typename TyT>
360 const 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 
369 template<typename TyT>
370 TyT 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 
462 template<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 
481 template<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 
503 template<typename TyT>
504 TyT 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 
588 template<typename TyT>
589 void 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 }
Lorene prototypes.
Definition: app_hor.h:67
TyT ** val
Array of pointers onto the values (size = size).
Definition: evolution.h:137
Evolution(const TyT &initial_value, int initial_j, double initial_time, int initial_size)
Constructor from some initial value.
Definition: evolution.C:127
int * step
Array of time step indices (size = size).
Definition: evolution.h:131
double * the_time
Array of values of t at the various time steps (size = size).
Definition: evolution.h:134
Time evolution (*** under development ***).
Definition: evolution.h:120