AdaptInfo.h 17.4 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  crystal growth group                                                  ==
// ==                                                                        ==
// ==  Stiftung caesar                                                       ==
// ==  Ludwig-Erhard-Allee 2                                                 ==
// ==  53175 Bonn                                                            ==
// ==  germany                                                               ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  http://www.caesar.de/cg/AMDiS                                         ==
// ==                                                                        ==
// ============================================================================

/** \file AdaptInfo.h */

#ifndef AMDIS_ADAPTINFO_H
#define AMDIS_ADAPTINFO_H

#include "MatrixVector.h"
#include "Parameters.h"
#include "Serializable.h"

namespace AMDiS {

  /**
   * \ingroup Adaption
   * 
   * \brief
   * Holds adapt parameters and infos about the problem. Base class
   * for AdaptInfoScal and AdaptInfoVec.
   */
  class AdaptInfo : public Serializable
  {
  protected:
    /** \brief
     * Stores adapt infos for a scalar problem or for one component of a 
     * vector valued problem.
     */
    class ScalContent {
    public:
47
      /// Constructor.
Thomas Witkowski's avatar
Thomas Witkowski committed
48
      ScalContent(std::string prefix) 
49
50
51
52
53
54
55
56
57
58
	: est_sum(0.0),
	  est_t_sum(0.0),
	  est_max(0.0),
	  est_t_max(0.0),
	  spaceTolerance(1.0),
	  timeTolerance(1.0),
	  timeErrLow(1.0),
	  coarsenAllowed(0),
	  refinementAllowed(1),
	  refineBisections(1),
59
  	  coarseBisections(1)
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
      {
	double totalTol = 1.0;
	double relSpaceErr = 1.0;
	double relTimeErr = 0.5;
	double timeTheta1 = 1.0;
	double timeTheta2 = 0.3;

	GET_PARAMETER(0, prefix + "->tolerance", "%f", &totalTol);
	GET_PARAMETER(0, prefix + "->rel space error", "%f", &relSpaceErr);
	GET_PARAMETER(0, prefix + "->rel time error", "%f", &relTimeErr);
	GET_PARAMETER(0, prefix + "->time theta 1", "%f", &timeTheta1);
	GET_PARAMETER(0, prefix + "->time theta 2", "%f", &timeTheta2);
	GET_PARAMETER(0, prefix + "->coarsen allowed", "%d", &coarsenAllowed);
	GET_PARAMETER(0, prefix + "->refinement allowed", "%d", &refinementAllowed);
	GET_PARAMETER(0, prefix + "->refine bisections", "%d", &refineBisections);
	GET_PARAMETER(0, prefix + "->coarsen bisections", "%d", &coarseBisections);

	spaceTolerance = totalTol * relSpaceErr;
	timeTolerance = totalTol * relTimeErr * timeTheta1;
	timeErrLow = totalTol * relTimeErr * timeTheta2;
80
      }
81

82
      /// Sum of all error estimates
83
84
      double est_sum;

85
      /// Sum of all time error estimates
86
87
      double est_t_sum;

88
      /// maximal local error estimate.
89
90
      double est_max;

91
      /// Maximum of all time error estimates
92
93
      double est_t_max;

94
      /// Tolerance for the (absolute or relative) error
95
96
      double spaceTolerance;

97
      /// Time tolerance.
98
99
      double timeTolerance;

100
      /// Lower bound for the time error.
101
102
      double timeErrLow;

103
      /// true if coarsening is allowed, false otherwise.
104
105
      int coarsenAllowed;

106
      /// true if refinement is allowed, false otherwise.
107
108
109
110
111
112
113
      int refinementAllowed;

      /** \brief
       * parameter to tell the marking strategy how many bisections should be 
       * performed when an element is marked for refinement; usually the value is
       * 1 or DIM
       */
114
      int refineBisections;
115
116
117
118
119
120

      /** \brief
       * parameter to tell the marking strategy how many bisections should
       * be undone when an element is marked for coarsening; usually the value is 
       * 1 or DIM
       */                          
121
      int coarseBisections;    
122
123
124
    };

  public:
125
    /// Constructor.
Thomas Witkowski's avatar
Thomas Witkowski committed
126
    AdaptInfo(std::string name_, int size = 1) 
127
128
129
130
131
132
133
134
135
136
137
138
139
140
      : name(name_), 
	spaceIteration(-1),
	maxSpaceIteration(-1),
	timestepIteration(0),
	maxTimestepIteration(30),
	timeIteration(0),
	maxTimeIteration(30),
	time(0.0),
	startTime(0.0),
	endTime(1.0),
	timestep(0.0),
	minTimestep(0.0),
	maxTimestep(1.0),
	timestepNumber(0),
Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
141
	nTimesteps(0),
142
143
144
145
	solverIterations(0),
	maxSolverIterations(0),
	solverTolerance(1e-8),
	solverResidual(0.0),
146
147
        scalContents(size),
	isDeserialized_(false)
148
149
150
151
152
153
154
155
156
157
158
159
    {
      GET_PARAMETER(0, name_ + "->start time", "%f", &startTime);
      time = startTime;
      GET_PARAMETER(0, name_ + "->timestep", "%f", &timestep);
      GET_PARAMETER(0, name_ + "->end time", "%f", &endTime);
      GET_PARAMETER(0, name_ + "->max iteration", "%d", &maxSpaceIteration);
      GET_PARAMETER(0, name_ + "->max timestep iteration", "%d", &maxTimestepIteration);
      GET_PARAMETER(0, name_ + "->max time iteration", "%d", &maxTimeIteration);

      GET_PARAMETER(0, name_ + "->min timestep", "%f", &minTimestep);
      GET_PARAMETER(0, name_ + "->max timestep", "%f", &maxTimestep);

Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
160
161
      GET_PARAMETER(0, name_ + "->number of timesteps", "%d", &nTimesteps);

162
163
164
165
166
167
      if (size == 1) {
	scalContents[0] = new ScalContent(name); 
      } else {
	char number[5];
	for (int i = 0; i < size; i++) {
	  sprintf(number, "[%d]", i);
Thomas Witkowski's avatar
Thomas Witkowski committed
168
	  scalContents[i] = new ScalContent(name + std::string(number)); 
169
170
	}
      }
171
    }
172

173
    /// Destructor.
174
    virtual ~AdaptInfo() {
175
      for (int i = 0;  i < scalContents.getSize(); i++)
176
	delete scalContents[i];
177
    }
178
179
180
181
182
183
184
185
186
187
188

    inline void reset() 
    {
      spaceIteration = -1;
      timestepIteration = 0;
      timeIteration = 0;
      time = 0.0;
      timestep = 0.0;
      timestepNumber = 0;
      solverIterations = 0;
      solverResidual = 0.0;
189
190

      GET_PARAMETER(0, name + "->timestep", "%f", &timestep);
191
    }
192

193
    /// Returns whether space tolerance is reached.
194
    virtual bool spaceToleranceReached() {
195
      int size = scalContents.getSize();
196
197
      for (int i = 0; i < size; i++)
	if (!(scalContents[i]->est_sum < scalContents[i]->spaceTolerance))
198
	  return false;
199

200
      return true;
201
    }
202

203
    /// Returns whether space tolerance of component i is reached.
204
    virtual bool spaceToleranceReached(int i) {
205
      if (!(scalContents[i]->est_sum < scalContents[i]->spaceTolerance))
206
	return false;
207
      else
208
	return true;
209
    }
210

211
    /// Returns whether time tolerance is reached.
212
    virtual bool timeToleranceReached() {
213
      int size = scalContents.getSize();
214
215
      for (int i = 0; i < size; i++)
	if (!(scalContents[i]->est_t_sum < scalContents[i]->timeTolerance))
216
	  return false;
217

218
      return true;
219
    }
220

221
    /// Returns whether time tolerance of component i is reached.
222
    virtual bool timeToleranceReached(int i) {
223
      if (!(scalContents[i]->est_t_sum < scalContents[i]->timeTolerance))
224
	return false;
225
      else
226
	return true;
227
    }
228

229
    /// Returns whether time error is under its lower bound.
230
231
    virtual bool timeErrorLow() {
      int size = scalContents.getSize();
232
233
      for (int i = 0; i < size; i++)
	if (!(scalContents[i]->est_t_sum < scalContents[i]->timeErrLow))
234
	  return false;
235

236
      return true;
237
238
    }

239
    /// Print debug information about time error and its bound.
240
    void printTimeErrorLowInfo() {
241
      for (int i = 0; i < scalContents.getSize(); i++)
242
243
244
	std::cout << "    Time error estimate = " << scalContents[i]->est_t_sum 
		  << "    Time error bound = " << scalContents[i]->timeErrLow << "\n";
    }
245

246
    /// Returns \ref spaceIteration.
247
248
    inline int getSpaceIteration() { 
      return spaceIteration; 
249
    }
250

251
    /// Sets \ref spaceIteration.
252
253
    inline void setSpaceIteration(int it) { 
      spaceIteration = it; 
254
    }
255
  
256
    /// Returns \ref maxSpaceIteration.
257
258
    inline int getMaxSpaceIteration() { 
      return maxSpaceIteration;
259
    }
260

261
    /// Sets \ref maxSpaceIteration.
262
263
    inline void setMaxSpaceIteration(int it) { 
      maxSpaceIteration = it; 
264
    }
265
  
266
    /// Increments \ref spaceIteration by 1;
267
268
    inline void incSpaceIteration() { 
      spaceIteration++; 
269
    }
270

271
    /// Sets \ref timestepIteration.
272
273
    inline void setTimestepIteration(int it) { 
      timestepIteration = it; 
274
    }
275
  
276
    /// Returns \ref timestepIteration.
277
278
    inline int getTimestepIteration() { 
      return timestepIteration; 
279
    }
280

281
    /// Increments \ref timestepIteration by 1;
282
283
    inline void incTimestepIteration() { 
      timestepIteration++; 
284
    }
285

286
    /// Returns \ref maxTimestepIteration.
287
288
    inline int getMaxTimestepIteration() { 
      return maxTimestepIteration; 
289
    }
290

291
    /// Sets \ref maxTimestepIteration.
292
293
    inline void setMaxTimestepIteration(int it) { 
      maxTimestepIteration = it; 
294
    }
295
  
296
    /// Sets \ref timeIteration.
297
298
    inline void setTimeIteration(int it) { 
      timeIteration = it; 
299
    }
300
  
301
    /// Returns \ref timeIteration.
302
303
    inline int getTimeIteration() { 
      return timeIteration; 
304
    }
305

306
    /// Increments \ref timesIteration by 1;
307
308
    inline void incTimeIteration() { 
      timeIteration++; 
309
    }
310

311
    /// Returns \ref maxTimeIteration.
312
313
    inline int getMaxTimeIteration() { 
      return maxTimeIteration; 
314
    }
315

316
    /// Sets \ref maxTimeIteration.
317
318
    inline void setMaxTimeIteration(int it) { 
      maxTimeIteration = it; 
319
    }
320
  
321
    /// Returns \ref timestepNumber.
322
323
    inline int getTimestepNumber() { 
      return timestepNumber; 
324
    }
325

Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
326
327
328
329
330
    /// Returns \ref nTimesteps.
    inline int getNumberOfTimesteps() {
      return nTimesteps;
    }

331
    /// Increments \ref timestepNumber by 1;
332
333
    inline void incTimestepNumber() { 
      timestepNumber++; 
334
    }
335

336
    /// Sets \ref est_sum.
337
338
    inline void setEstSum(double e, int index) {
      scalContents[index]->est_sum = e;
339
    }
340

341
    /// Sets \ref est_max.
342
343
    inline void setEstMax(double e, int index) {
      scalContents[index]->est_max = e;
344
    }
345

346
    /// Sets \ref est_max.
347
348
    inline void setTimeEstMax(double e, int index) {
      scalContents[index]->est_t_max = e;
349
    }
350

351
    /// Sets \ref est_t_sum.
352
353
    inline void setTimeEstSum(double e, int index) {
      scalContents[index]->est_t_sum = e;
354
    }
355

356
    /// Returns \ref est_sum.
357
358
    inline double getEstSum(int index) { 
      return scalContents[index]->est_sum; 
359
    }
360

361
    /// Returns \ref est_t_sum.
362
363
    inline double getEstTSum(int index) { 
      return scalContents[index]->est_t_sum; 
364
    }
365

366
    /// Returns \ref est_max.
367
368
    inline double getEstMax(int index) { 
      return scalContents[index]->est_max; 
369
    }
370

371
    /// Returns \ref est_max.
372
373
    inline double getTimeEstMax(int index) { 
      return scalContents[index]->est_t_max; 
374
    }
375

376
    /// Returns \ref est_t_sum.
377
378
    inline double getTimeEstSum(int index) { 
      return scalContents[index]->est_t_sum; 
379
    }
380

381
    /// Returns \ref spaceTolerance.
382
383
    inline double getSpaceTolerance(int index) { 
      return scalContents[index]->spaceTolerance; 
384
    }  
385

386
    /// Sets \ref spaceTolerance.
387
388
    inline void setSpaceTolerance(int index, double tol) { 
      scalContents[index]->spaceTolerance = tol; 
389
    }  
390

391
    /// Returns \ref timeTolerance.
392
393
    inline double getTimeTolerance(int index) { 
      return scalContents[index]->timeTolerance; 
394
    }  
395

396
    /// Sets \ref time
397
398
    inline double setTime(double t) { 
      time = t; 
399
400
401
402
403
      if (time > endTime) 
	time = endTime;
      if (time < startTime) 
	time = startTime;

404
      return time;
405
    }
406

407
    /// Gets \ref time
408
409
    inline double getTime() { 
      return time; 
410
    }  
411

412
    /// Gets \ref &time
413
414
    inline double *getTimePtr() { 
      return &time; 
415
    }  
416

417
    /// Sets \ref timestep
418
419
    inline double setTimestep(double t) { 
      timestep = t; 
420
      if (timestep > maxTimestep) {
421
	timestep = maxTimestep;
422
423
      }
      if (timestep < minTimestep) {
424
	timestep = minTimestep;
425
426
      }
      if (time + timestep > endTime) {
427
	timestep = endTime - time;
428
      }
429
      
430
      return timestep;
431
    }
432

Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
433
434
435
436
437
    /** \brief
     * Returns true, if the end time is reached and no more timestep
     * computations must be done.
     */
    inline bool reachedEndTime() {
438
439
440
441
      if (nTimesteps > 0) 
	return !(timestepNumber < nTimesteps);

      return !(time < endTime - DBL_TOL);
Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
442
443
    }

444
    /// Gets \ref timestep
445
446
    inline double getTimestep() { 
      return timestep; 
447
    }
448

449
    /// Sets \ref minTimestep
450
451
    inline void setMinTimestep(double t) { 
      minTimestep = t; 
452
    }
453

454
    /// Gets \ref minTimestep
455
456
    inline double getMinTimestep() { 
      return minTimestep; 
457
    }  
458

459
    /// Sets \ref maxTimestep
460
461
    inline void setMaxTimestep(double t) { 
      maxTimestep = t; 
462
    }
463

464
    /// Gets \ref maxTimestep
465
466
    inline double getMaxTimestep() { 
      return maxTimestep; 
467
    }  
468

469
    /// Gets \ref &timestep
470
471
    inline double *getTimestepPtr() { 
      return &timestep; 
472
    }  
473

474
    /// Sets \ref startTime = time
475
476
    inline void setStartTime(double time) { 
      startTime = time; 
477
    }
478

479
    /// Sets \ref endTime = time
480
481
    inline void setEndTime(double time) { 
      endTime = time; 
482
    }
483

484
    /// Returns \ref startTime
485
486
    inline double getStartTime() { 
      return startTime; 
487
    }
488

489
    /// Returns \ref endTime
490
491
    inline double getEndTime() { 
      return endTime; 
492
    }
493

494
    /// Returns \ref timeErrLow.
495
496
    inline double getTimeErrLow(int index) { 
      return scalContents[index]->timeErrLow; 
497
    }  
498

499
    /// Returns whether coarsening is allowed or not.
500
501
    inline bool isCoarseningAllowed(int index) {
      return (scalContents[index]->coarsenAllowed == 1);
502
    }
503

504
    /// Returns whether coarsening is allowed or not.
505
506
    inline bool isRefinementAllowed(int index) {
      return (scalContents[index]->refinementAllowed == 1);
507
    }
508

509
    ///
510
511
    inline void allowRefinement(bool allow, int index) {
      scalContents[index]->refinementAllowed = allow;
512
    }
513

514
    ///
515
516
    inline void allowCoarsening(bool allow, int index) {
      scalContents[index]->coarsenAllowed = allow;
517
    }
518

519
    /// Returns \ref refineBisections
520
521
    inline const int getRefineBisections(int index) const {
      return scalContents[index]->refineBisections;
522
    }
523

524
    /// Returns \ref coarseBisections
525
526
    inline const int getCoarseBisections(int index) const {
      return scalContents[index]->coarseBisections;
527
    }
528
529
530

    inline int getSize() { 
      return scalContents.getSize(); 
531
    }
532
533
534

    inline void setSolverIterations(int it) {
      solverIterations = it;
535
    }
536
537
538
  
    inline int getSolverIterations() {
      return solverIterations;
539
    }
540
541
542
  
    inline void setMaxSolverIterations(int it) {
      maxSolverIterations = it;
543
    }
544
545
546
  
    inline int getMaxSolverIterations() {
      return maxSolverIterations;
547
    }
548
549
550
  
    inline void setSolverTolerance(double tol) {
      solverTolerance = tol;
551
    }
552
553
554
  
    inline double getSolverTolerance() {
      return solverTolerance;
555
    }
556
557
558
  
    inline void setSolverResidual(double res) {
      solverResidual = res;
559
    }
560
561
562
  
    inline double getSolverResidual() {
      return solverResidual;
563
    }
564

565
    /// Returns true, if the adaptive procedure was deserialized from a file.
566
567
568
569
570
571
572
573
    const bool isDeserialized() const {
      return isDeserialized_;
    }

    inline void setIsDeserialized(bool b) {
      isDeserialized_ = b;
    }

574
    /// Creates new scalContents with the given size.
575
576
    void setScalContents(int newSize);

577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
    /** \brief
     * Resets timestep, current time and time boundaries without
     * any check. Is used by the parareal algorithm.
     */
    void resetTimeValues(double newTimeStep,
			 double newStartTime,
			 double newEndTime)
    {
      time = newStartTime;
      startTime = newStartTime;
      endTime = newEndTime;
      timestep = newTimeStep;
      timestepNumber = 0;
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
592
    void serialize(std::ostream& out);
593

Thomas Witkowski's avatar
Thomas Witkowski committed
594
    void deserialize(std::istream& in);
595
596

  protected:
597
    /// Name.
Thomas Witkowski's avatar
Thomas Witkowski committed
598
    std::string name;
599

600
    /// Current space iteration
601
602
603
604
605
606
607
608
    int spaceIteration;

    /** \brief
     * maximal allowed number of iterations of the adaptive procedure; if 
     * maxIteration <= 0, no iteration bound is used
     */
    int maxSpaceIteration;

609
    /// Current timestep iteration
610
611
    int timestepIteration;

612
    /// Maximal number of iterations for choosing a timestep
613
614
    int maxTimestepIteration;

615
    /// Current time iteration
616
617
    int timeIteration;

618
    /// Maximal number of time iterations
619
620
    int maxTimeIteration;

621
    /// Actual time, end of time interval for current time step
622
623
    double time;

624
    /// Initial time
625
626
    double startTime;

627
    /// Final time
628
629
    double endTime;

630
    /// Current time step size
631
632
    double timestep;

633
    /// Minimal step size
634
635
    double minTimestep;

636
    /// Maximal step size
637
638
    double maxTimestep;

639
    /// Number of current time step
640
    int timestepNumber;
Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
641
642
643
644
645
646
647

    /** \brief
     * Per default this value is 0 and not used. If it is set to a non-zero value,
     * the computation of the stationary problem is done nTimesteps times with a
     * fixed timestep.
     */
    int nTimesteps;
648
  
649
    /// number of iterations needed of linear or nonlinear solver
650
651
    int solverIterations;

652
    /// maximal number of iterations needed of linear or nonlinear solver
653
654
    int maxSolverIterations;

655
    ///
656
657
    double solverTolerance;

658
    ///
659
660
    double solverResidual;

661
    /// Scalar adapt infos.
662
    Vector<ScalContent*> scalContents;
663

664
    /// Is true, if the adaptive procedure was deserialized from a file.
665
    bool isDeserialized_;
666
667
668
669
670
  };

}

#endif //  AMDIS_ADAPTINFO_H