AdaptInfo.h 18.7 KB
Newer Older
1
2
3
4
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
5
// ==  http://www.amdis-fem.org                                              ==
6
7
// ==                                                                        ==
// ============================================================================
8
9
10
11
12
13
14
15
16
17
18
19
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


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
47

/** \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:
48
      /// Constructor.
Thomas Witkowski's avatar
Thomas Witkowski committed
49
      ScalContent(std::string prefix) 
50
	: est_sum(0.0),
51
52
53
54
55
56
57
58
59
60
61
62
	  est_t_sum(0.0),
	  est_max(0.0),
	  est_t_max(0.0),
	  fac_max(0.0),
	  fac_sum(1.0),
	  spaceTolerance(0.0),
	  timeTolerance(0.0),
	  timeErrLow(0.0),
	  coarsenAllowed(0),
	  refinementAllowed(1),
	  refineBisections(1),
	  coarseBisections(1)	  	
63
      {
64
65
	GET_PARAMETER(0, prefix + "->tolerance", "%f", &spaceTolerance);
	GET_PARAMETER(0, prefix + "->time tolerance", "%f", &timeTolerance);
66
67
68
69
	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);
70
71
	GET_PARAMETER(0, prefix + "->sum factor", "%f", &fac_sum);
	GET_PARAMETER(0, prefix + "->max factor", "%f", &fac_max);
72

73
	timeErrLow = timeTolerance * 0.3;
74
      }
75

76
      /// Sum of all error estimates
77
78
      double est_sum;

79
      /// Sum of all time error estimates
80
81
      double est_t_sum;

82
      /// maximal local error estimate.
83
84
      double est_max;

85
      /// Maximum of all time error estimates
86
      double est_t_max;
87
88
89
      
      /// factors to combine max and integral time estimate
      double fac_max, fac_sum;
90

91
      /// Tolerance for the (absolute or relative) error
92
93
      double spaceTolerance;

94
      /// Time tolerance.
95
96
      double timeTolerance;

97
      /// Lower bound for the time error.
98
99
      double timeErrLow;

100
      /// true if coarsening is allowed, false otherwise.
101
102
      int coarsenAllowed;

103
      /// true if refinement is allowed, false otherwise.
104
105
106
107
108
109
110
      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
       */
111
      int refineBisections;
112
113
114
115
116
117

      /** \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
       */                          
118
      int coarseBisections;    
119
120
121
    };

  public:
122
    /// Constructor.
Thomas Witkowski's avatar
Thomas Witkowski committed
123
    AdaptInfo(std::string name_, int size = 1) 
124
125
126
127
128
129
130
131
132
133
134
      : 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),
135
	lastProcessedTimestep(0.0),
136
137
138
	minTimestep(0.0),
	maxTimestep(1.0),
	timestepNumber(0),
Thomas Witkowski's avatar
* Bla  
Thomas Witkowski committed
139
	nTimesteps(0),
140
141
142
143
	solverIterations(0),
	maxSolverIterations(0),
	solverTolerance(1e-8),
	solverResidual(0.0),
144
        scalContents(size),
145
146
	deserialized(false),
	rosenbrockMode(false)
147
148
149
150
151
152
153
154
155
156
157
158
    {
      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
159
160
      GET_PARAMETER(0, name_ + "->number of timesteps", "%d", &nTimesteps);

161
162
163
164
165
166
      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
167
	  scalContents[i] = new ScalContent(name + std::string(number)); 
168
169
	}
      }
170
    }
171

172
    /// Destructor.
Thomas Witkowski's avatar
Thomas Witkowski committed
173
174
175
    virtual ~AdaptInfo() 
    {
      for (unsigned int i = 0;  i < scalContents.size(); 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
      lastProcessedTimestep=timestep;
192
    }
193

194
    /// Returns whether space tolerance is reached.
195
196
    virtual bool spaceToleranceReached() 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
197
198
199
200
      for (unsigned int i = 0; i < scalContents.size(); i++) {
	std::cout << "est_sum:" <<scalContents[i]->est_sum 
		  << " spaceTol: " << scalContents[i]->spaceTolerance 
		  << std::endl;
201
	if (!(scalContents[i]->est_sum < scalContents[i]->spaceTolerance))
202
	  return false;
203
      }
204

205
      return true;
206
    }
207

208
    /// Returns whether space tolerance of component i is reached.
209
210
    virtual bool spaceToleranceReached(int i) 
    {
211
      if (!(scalContents[i]->est_sum < scalContents[i]->spaceTolerance))
212
	return false;
213
      else
214
	return true;
215
    }
216

217
    /// Returns whether time tolerance is reached.
218
219
    virtual bool timeToleranceReached() 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
220
      for (unsigned int i = 0; i < scalContents.size(); i++)
221
	if (!(getTimeEstCombined(i) < scalContents[i]->timeTolerance))
222
	  return false;
223

224
      return true;
225
    }
226

227
    /// Returns whether time tolerance of component i is reached.
228
229
    virtual bool timeToleranceReached(int i) 
    {
230
      if (!(getTimeEstCombined(i) < scalContents[i]->timeTolerance))
231
	return false;
232
      else
233
	return true;
234
    }
235

236
    /// Returns whether time error is under its lower bound.
237
238
    virtual bool timeErrorLow() 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
239
      for (unsigned int i = 0; i < scalContents.size(); i++)
240
	if (!(getTimeEstCombined(i) < scalContents[i]->timeErrLow))
241
	  return false;
242

243
      return true;
244
    }
245
246
247
248
    /// Returns the time estimation as a combination 
    /// of maximal and integral time error 
    double getTimeEstCombined(unsigned i) const 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
249
250
251
      return 
	scalContents[i]->est_t_max * scalContents[i]->fac_max +
	scalContents[i]->est_t_sum * scalContents[i]->fac_sum; 
252
253
    }

254

255
    /// Print debug information about time error and its bound.
256
257
    void printTimeErrorLowInfo() 
    {
258
259
260
261
      for (unsigned int i = 0; i < scalContents.size(); i++){
	std::cout << "    Time error estimate  = " << getTimeEstCombined(i)
		  << "    Time error estimate sum = " << scalContents[i]->est_t_sum 
		  << "    Time error estimate max = " << scalContents[i]->est_t_max 
262
263
		  << "    Time error low bound = " << scalContents[i]->timeErrLow  
		  << "    Time error high bound = " << scalContents[i]->timeTolerance << "\n";
264
      }
265
    }
266

267
    /// Returns \ref spaceIteration.
268
269
    inline int getSpaceIteration() 
    { 
270
      return spaceIteration; 
271
    }
272

273
    /// Sets \ref spaceIteration.
274
275
    inline void setSpaceIteration(int it) 
    { 
276
      spaceIteration = it; 
277
    }
278
  
279
    /// Returns \ref maxSpaceIteration.
280
281
    inline int getMaxSpaceIteration() 
    { 
282
      return maxSpaceIteration;
283
    }
284

285
    /// Sets \ref maxSpaceIteration.
286
287
    inline void setMaxSpaceIteration(int it) 
    { 
288
      maxSpaceIteration = it; 
289
    }
290
  
291
    /// Increments \ref spaceIteration by 1;
292
293
    inline void incSpaceIteration() 
    { 
294
      spaceIteration++; 
295
    }
296

297
    /// Sets \ref timestepIteration.
298
299
    inline void setTimestepIteration(int it) 
    { 
300
      timestepIteration = it; 
301
    }
302
  
303
    /// Returns \ref timestepIteration.
304
305
    inline int getTimestepIteration() 
    { 
306
      return timestepIteration; 
307
    }
308

309
    /// Increments \ref timestepIteration by 1;
310
311
    inline void incTimestepIteration() 
    { 
312
      timestepIteration++; 
313
    }
314

315
    /// Returns \ref maxTimestepIteration.
316
317
    inline int getMaxTimestepIteration() 
    { 
318
      return maxTimestepIteration; 
319
    }
320

321
    /// Sets \ref maxTimestepIteration.
322
323
    inline void setMaxTimestepIteration(int it) 
    { 
324
      maxTimestepIteration = it; 
325
    }
326
  
327
    /// Sets \ref timeIteration.
328
329
    inline void setTimeIteration(int it) 
    { 
330
      timeIteration = it; 
331
    }
332
  
333
    /// Returns \ref timeIteration.
334
335
    inline int getTimeIteration() 
    { 
336
      return timeIteration; 
337
    }
338

339
    /// Increments \ref timesIteration by 1;
340
341
    inline void incTimeIteration() 
    { 
342
      timeIteration++; 
343
    }
344

345
    /// Returns \ref maxTimeIteration.
346
347
    inline int getMaxTimeIteration() 
    { 
348
      return maxTimeIteration; 
349
    }
350

351
    /// Sets \ref maxTimeIteration.
352
353
    inline void setMaxTimeIteration(int it) 
    { 
354
      maxTimeIteration = it; 
355
    }
356
  
357
    /// Returns \ref timestepNumber.
358
359
    inline int getTimestepNumber() 
    { 
360
      return timestepNumber; 
361
    }
362

363
364
365
    /// Sets \ref timestepNumber.
    inline void setTimestepNumber(int num) 
    {
366
      timestepNumber = std::min(nTimesteps, num);
367
368
    }
    
Thomas Witkowski's avatar
* Bla  
Thomas Witkowski committed
369
    /// Returns \ref nTimesteps.
370
371
    inline int getNumberOfTimesteps() 
    {
Thomas Witkowski's avatar
* Bla  
Thomas Witkowski committed
372
373
374
      return nTimesteps;
    }

375
376
377
    /// Sets \ref nTimesteps.
    inline void setNumberOfTimesteps(int num) 
    {
378
      nTimesteps = std::max(0, num);
379
380
    }

381
    /// Increments \ref timestepNumber by 1;
382
383
    inline void incTimestepNumber() 
    { 
384
      timestepNumber++; 
385
    }
386

387
    /// Sets \ref est_sum.
388
389
    inline void setEstSum(double e, int index) 
    {
390
      scalContents[index]->est_sum = e;
391
    }
392

393
    /// Sets \ref est_max.
394
395
    inline void setEstMax(double e, int index) 
    {
396
      scalContents[index]->est_max = e;
397
    }
398

399
    /// Sets \ref est_max.
400
401
    inline void setTimeEstMax(double e, int index) 
    {
402
      scalContents[index]->est_t_max = e;
403
    }
404

405
    /// Sets \ref est_t_sum.
406
407
    inline void setTimeEstSum(double e, int index) 
    {
408
      scalContents[index]->est_t_sum = e;
409
    }
410

411
    /// Returns \ref est_sum.
412
413
    inline double getEstSum(int index) 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
414
415
416
417
418
      FUNCNAME("AdaptInfo::getEstSum()");

      TEST_EXIT_DBG(static_cast<unsigned int>(index) < scalContents.size())
	("Wrong index for adaptInfo!\n");

419
      return scalContents[index]->est_sum; 
420
    }
421

422
    /// Returns \ref est_t_sum.
423
424
    inline double getEstTSum(int index) 
    { 
425
      return scalContents[index]->est_t_sum; 
426
    }
427

428
    /// Returns \ref est_max.
429
430
    inline double getEstMax(int index) 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
431
432
433
434
435
      FUNCNAME("AdaptInfo::getEstSum()");

      TEST_EXIT_DBG(static_cast<unsigned int>(index) < scalContents.size())
	("Wrong index for adaptInfo!\n");

436
      return scalContents[index]->est_max; 
437
    }
438

439
    /// Returns \ref est_max.
440
441
    inline double getTimeEstMax(int index) 
    { 
442
      return scalContents[index]->est_t_max; 
443
    }
444

445
    /// Returns \ref est_t_sum.
446
447
    inline double getTimeEstSum(int index) 
    { 
448
      return scalContents[index]->est_t_sum; 
449
    }
450

451
    /// Returns \ref spaceTolerance.
452
453
    inline double getSpaceTolerance(int index) 
    { 
454
      return scalContents[index]->spaceTolerance; 
455
    }  
456

457
    /// Sets \ref spaceTolerance.
458
459
    inline void setSpaceTolerance(int index, double tol) 
    { 
460
      scalContents[index]->spaceTolerance = tol; 
461
    }  
462

463
    /// Returns \ref timeTolerance.
464
465
    inline double getTimeTolerance(int index) 
    { 
466
      return scalContents[index]->timeTolerance; 
467
    }  
468

469
    /// Sets \ref time
470
471
    inline double setTime(double t) 
    { 
472
      time = t; 
473
474
475
476
477
      if (time > endTime) 
	time = endTime;
      if (time < startTime) 
	time = startTime;

478
      return time;
479
    }
480

481
    /// Gets \ref time
482
483
    inline double getTime() 
    { 
484
      return time; 
485
    }  
486

487
    /// Gets \ref &time
488
489
    inline double *getTimePtr() 
    { 
490
      return &time; 
491
    }  
492

493
    /// Sets \ref timestep
494
495
    inline double setTimestep(double t) 
    { 
496
      timestep = t; 
497
      if (timestep > maxTimestep)
498
	timestep = maxTimestep;
499
      if (timestep < minTimestep)
500
	timestep = minTimestep;
501
      if (time + timestep > endTime)
502
503
	timestep = endTime - time;
      
504
      return timestep;
505
    }
506
507
508
509
510
511
512
513
514
515
516
517
518
    /// Gets \ref timestep
    inline double getTimestep() 
    { 
      return timestep; 
    }

    inline void setLastProcessedTimestep(double t){
	lastProcessedTimestep=t;
    } 

    inline double getLastProcessedTimestep(){
	return lastProcessedTimestep;
    } 
519

Thomas Witkowski's avatar
* Bla  
Thomas Witkowski committed
520
521
522
523
    /** \brief
     * Returns true, if the end time is reached and no more timestep
     * computations must be done.
     */
524
525
    inline bool reachedEndTime() 
    {
526
527
528
529
      if (nTimesteps > 0) 
	return !(timestepNumber < nTimesteps);

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

532

533
    /// Sets \ref minTimestep
534
535
    inline void setMinTimestep(double t) 
    { 
536
      minTimestep = t; 
537
    }
538

539
    /// Gets \ref minTimestep
540
541
    inline double getMinTimestep() 
    { 
542
      return minTimestep; 
543
    }  
544

545
    /// Sets \ref maxTimestep
546
547
    inline void setMaxTimestep(double t) 
    { 
548
      maxTimestep = t; 
549
    }
550

551
    /// Gets \ref maxTimestep
552
553
    inline double getMaxTimestep() 
    { 
554
      return maxTimestep; 
555
    }  
556
    
557
    /// Gets \ref &timestep
558
559
    inline double *getTimestepPtr() 
    { 
560
      return &timestep; 
561
    }  
562

563
    /// Sets \ref startTime = time
564
565
    inline void setStartTime(double time) 
    { 
566
      startTime = time; 
567
    }
568

569
    /// Sets \ref endTime = time
570
571
    inline void setEndTime(double time) 
    { 
572
      endTime = time; 
573
    }
574

575
    /// Returns \ref startTime
576
577
    inline double getStartTime() 
    { 
578
      return startTime; 
579
    }
580

581
    /// Returns \ref endTime
582
583
    inline double getEndTime() 
    { 
584
      return endTime; 
585
    }
586

587
    /// Returns \ref timeErrLow.
588
589
    inline double getTimeErrLow(int index) 
    { 
590
      return scalContents[index]->timeErrLow; 
591
    }  
592

593
    /// Returns whether coarsening is allowed or not.
594
595
    inline bool isCoarseningAllowed(int index) 
    {
596
      return (scalContents[index]->coarsenAllowed == 1);
597
    }
598

599
    /// Returns whether coarsening is allowed or not.
600
601
    inline bool isRefinementAllowed(int index) 
    {
602
      return (scalContents[index]->refinementAllowed == 1);
603
    }
604

605
    ///
606
607
    inline void allowRefinement(bool allow, int index) 
    {
608
      scalContents[index]->refinementAllowed = allow;
609
    }
610

611
    ///
612
613
    inline void allowCoarsening(bool allow, int index) 
    {
614
      scalContents[index]->coarsenAllowed = allow;
615
    }
616

617
    /// Returns \ref refineBisections
618
619
    inline const int getRefineBisections(int index) const 
    {
620
      return scalContents[index]->refineBisections;
621
    }
622

623
    /// Returns \ref coarseBisections
624
625
    inline const int getCoarseBisections(int index) const 
    {
626
      return scalContents[index]->coarseBisections;
627
    }    
628

629
630
    inline int getSize() 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
631
      return scalContents.size(); 
632
    }
633

634
635
636
637
638
    inline bool getRosenbrockMode()
    {
      return rosenbrockMode;
    }

639
640
    inline void setSolverIterations(int it) 
    {
641
      solverIterations = it;
642
    }
643
  
644
645
    inline int getSolverIterations() 
    {
646
      return solverIterations;
647
    }
648
  
649
650
    inline void setMaxSolverIterations(int it) 
    {
651
      maxSolverIterations = it;
652
    }
653
  
654
655
    inline int getMaxSolverIterations() 
    {
656
      return maxSolverIterations;
657
    }
658
  
659
660
    inline void setSolverTolerance(double tol) 
    {
661
      solverTolerance = tol;
662
    }
663
  
664
665
    inline double getSolverTolerance() 
    {
666
      return solverTolerance;
667
    }
668
  
669
670
    inline void setSolverResidual(double res) 
    {
671
      solverResidual = res;
672
    }
673
  
674
675
    inline double getSolverResidual() 
    {
676
      return solverResidual;
677
    }
678

679
    /// Returns true, if the adaptive procedure was deserialized from a file.
680
681
    const bool isDeserialized() const 
    {
682
      return deserialized;
683
684
    }

685
686
    inline void setIsDeserialized(bool b) 
    {
687
688
689
690
691
692
      deserialized = b;
    }

    inline void setRosenbrockMode(bool b)
    {
      rosenbrockMode = b;
693
694
    }

695
    /// Creates new scalContents with the given size.
696
697
    void setScalContents(int newSize);

698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
    /** \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
713
    void serialize(std::ostream& out);
714

Thomas Witkowski's avatar
Thomas Witkowski committed
715
    void deserialize(std::istream& in);
716
717

  protected:
718
    /// Name.
Thomas Witkowski's avatar
Thomas Witkowski committed
719
    std::string name;
720

721
    /// Current space iteration
722
723
724
725
726
727
728
729
    int spaceIteration;

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

730
    /// Current timestep iteration
731
732
    int timestepIteration;

733
    /// Maximal number of iterations for choosing a timestep
734
735
    int maxTimestepIteration;

736
    /// Current time iteration
737
738
    int timeIteration;

739
    /// Maximal number of time iterations
740
741
    int maxTimeIteration;

742
    /// Actual time, end of time interval for current time step
743
744
    double time;

745
    /// Initial time
746
747
    double startTime;

748
    /// Final time
749
750
    double endTime;

751
    ///Time step size to be used
752
753
    double timestep;

754