AdaptInfo.h 18.9 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

/** \file AdaptInfo.h */

#ifndef AMDIS_ADAPTINFO_H
#define AMDIS_ADAPTINFO_H

#include "MatrixVector.h"
27
#include "Initfile.h"
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
#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.
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
66
67
68
69
70
71
	Parameters::get(prefix + "->tolerance", spaceTolerance);
	Parameters::get(prefix + "->time tolerance", timeTolerance);
	Parameters::get(prefix + "->coarsen allowed", coarsenAllowed);
	Parameters::get(prefix + "->refinement allowed", refinementAllowed);
	Parameters::get(prefix + "->refine bisections", refineBisections);
	Parameters::get(prefix + "->coarsen bisections", coarseBisections);
	Parameters::get(prefix + "->sum factor", fac_sum);
	Parameters::get(prefix + "->max factor", 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
      init();
149
150
      char number[5];
      for (int i = 0; i < size; i++) {
151
152
	sprintf(number, "[%d]", i);
	scalContents[i] = new ScalContent(name + std::string(number));  
153
      }
154
    }
155

156
    /// Destructor.
Thomas Witkowski's avatar
Thomas Witkowski committed
157
158
159
    virtual ~AdaptInfo() 
    {
      for (unsigned int i = 0;  i < scalContents.size(); i++)
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
  delete scalContents[i];
    }

    /// Sets initial values to time/timestep variables
    inline void init()
    {
      Parameters::get(name + "->start time", startTime);
      time = startTime;
      Parameters::get(name + "->timestep", timestep);
      Parameters::get(name + "->end time", endTime);
      Parameters::get(name + "->max iteration", maxSpaceIteration);
      Parameters::get(name + "->max timestep iteration", maxTimestepIteration);
      Parameters::get(name + "->max time iteration", maxTimeIteration);
      Parameters::get(name + "->min timestep", minTimestep);
      Parameters::get(name + "->max timestep", maxTimestep);
      Parameters::get(name + "->number of timesteps", nTimesteps);
176
    }
177

178
    /// Resets all variables to zero (or something equivalent)
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
      Parameters::get(name + "->timestep", 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
      for (unsigned int i = 0; i < scalContents.size(); i++) {
198
	if (!(scalContents[i]->est_sum < scalContents[i]->spaceTolerance))
199
	  return false;
200
      }
201

202
      return true;
203
    }
204

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

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

221
      return true;
222
    }
223

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

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

240
      return true;
241
    }
242
243
244
245
    /// 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
246
247
248
      return 
	scalContents[i]->est_t_max * scalContents[i]->fac_max +
	scalContents[i]->est_t_sum * scalContents[i]->fac_sum; 
249
250
    }

251

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

421
      return scalContents[index]->est_sum; 
422
    }
423

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

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

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

438
      return scalContents[index]->est_max; 
439
    }
440

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

447
    /// Returns \ref est_t_sum.
448
449
    inline double getTimeEstSum(int index) 
    { 
450
      return scalContents[index]->est_t_sum; 
451
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
452
453
454
455
456
457
458
459
460
461
462
    
    /// Returns \ref est_t the estimated overall time error
    inline double getTimeEst()
    {
      return est_t;
    }
    
    inline void setTimeEst(double value)
    {
      est_t = value;
    }
463

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

470
    /// Sets \ref spaceTolerance.
471
472
    inline void setSpaceTolerance(int index, double tol) 
    { 
473
      scalContents[index]->spaceTolerance = tol; 
474
    }  
475

476
    /// Returns \ref timeTolerance.
477
478
    inline double getTimeTolerance(int index) 
    { 
479
      return scalContents[index]->timeTolerance; 
480
    }  
481

482
    /// Sets \ref time
483
484
    inline double setTime(double t) 
    { 
485
      time = t; 
486
487
488
489
490
      if (time > endTime) 
	time = endTime;
      if (time < startTime) 
	time = startTime;

491
      return time;
492
    }
493

494
    /// Gets \ref time
495
496
    inline double getTime() 
    { 
497
      return time; 
498
    }  
499

500
    /// Gets \ref &time
501
502
    inline double *getTimePtr() 
    { 
503
      return &time; 
504
    }  
505

506
    /// Sets \ref timestep
507
508
    inline double setTimestep(double t) 
    { 
509
      timestep = t; 
510
      if (timestep > maxTimestep)
511
	timestep = maxTimestep;
512
      if (timestep < minTimestep)
513
	timestep = minTimestep;
514
      if (time + timestep > endTime)
515
516
	timestep = endTime - time;
      
517
      return timestep;
518
    }
519
520
521
522
523
524
    /// Gets \ref timestep
    inline double getTimestep() 
    { 
      return timestep; 
    }

525
526
527
    inline void setLastProcessedTimestep(double t)
    {
      lastProcessedTimestep = t;
528
529
    } 

530
531
    inline double getLastProcessedTimestep()
    {
532
533
	return lastProcessedTimestep;
    } 
534

535
536
    /// Returns true, if the end time is reached and no more timestep
    /// computations must be done.
537
538
    inline bool reachedEndTime() 
    {
539
540
541
      if (nTimesteps > 0) 
	return !(timestepNumber < nTimesteps);

542
      return !(fabs(time - endTime) > DBL_TOL);
Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
543
544
    }

545

546
    /// Sets \ref minTimestep
547
548
    inline void setMinTimestep(double t) 
    { 
549
      minTimestep = t; 
550
    }
551

552
    /// Gets \ref minTimestep
553
554
    inline double getMinTimestep() 
    { 
555
      return minTimestep; 
556
    }  
557

558
    /// Sets \ref maxTimestep
559
560
    inline void setMaxTimestep(double t) 
    { 
561
      maxTimestep = t; 
562
    }
563

564
    /// Gets \ref maxTimestep
565
566
    inline double getMaxTimestep() 
    { 
567
      return maxTimestep; 
568
    }  
569
    
570
    /// Gets \ref &timestep
571
572
    inline double *getTimestepPtr() 
    { 
573
      return &timestep; 
574
    }  
575

576
    /// Sets \ref startTime = time
577
578
    inline void setStartTime(double time) 
    { 
579
      startTime = time; 
580
    }
581

582
    /// Sets \ref endTime = time
583
584
    inline void setEndTime(double time) 
    { 
585
      endTime = time; 
586
    }
587

588
    /// Returns \ref startTime
589
590
    inline double getStartTime() 
    { 
591
      return startTime; 
592
    }
593

594
    /// Returns \ref endTime
595
596
    inline double getEndTime() 
    { 
597
      return endTime; 
598
    }
599

600
    /// Returns \ref timeErrLow.
601
602
    inline double getTimeErrLow(int index) 
    { 
603
      return scalContents[index]->timeErrLow; 
604
    }  
605

606
    /// Returns whether coarsening is allowed or not.
607
608
    inline bool isCoarseningAllowed(int index) 
    {
609
      return (scalContents[index]->coarsenAllowed == 1);
610
    }
611

612
    /// Returns whether coarsening is allowed or not.
613
614
    inline bool isRefinementAllowed(int index) 
    {
615
      return (scalContents[index]->refinementAllowed == 1);
616
    }
617

618
    ///
619
620
    inline void allowRefinement(bool allow, int index) 
    {
621
      scalContents[index]->refinementAllowed = allow;
622
    }
623

624
    ///
625
626
    inline void allowCoarsening(bool allow, int index) 
    {
627
      scalContents[index]->coarsenAllowed = allow;
628
    }
629

630
    /// Returns \ref refineBisections
631
    inline int getRefineBisections(int index) const 
632
    {
633
      return scalContents[index]->refineBisections;
634
    }
635

636
    /// Returns \ref coarseBisections
637
    inline int getCoarseBisections(int index) const 
638
    {
639
      return scalContents[index]->coarseBisections;
640
    }    
641

642
643
    inline int getSize() 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
644
      return scalContents.size(); 
645
    }
646

647
648
649
650
651
    inline bool getRosenbrockMode()
    {
      return rosenbrockMode;
    }

652
653
    inline void setSolverIterations(int it) 
    {
654
      solverIterations = it;
655
    }
656
  
657
658
    inline int getSolverIterations() 
    {
659
      return solverIterations;
660
    }
661
  
662
663
    inline void setMaxSolverIterations(int it) 
    {
664
      maxSolverIterations = it;
665
    }
666
  
667
668
    inline int getMaxSolverIterations() 
    {
669
      return maxSolverIterations;
670
    }
671
  
672
673
    inline void setSolverTolerance(double tol) 
    {
674
      solverTolerance = tol;
675
    }
676
  
677
678
    inline double getSolverTolerance() 
    {
679
      return solverTolerance;
680
    }
681
  
682
683
    inline void setSolverResidual(double res) 
    {
684
      solverResidual = res;
685
    }
686
  
687
688
    inline double getSolverResidual() 
    {
689
      return solverResidual;
690
    }
691

692
    /// Returns true, if the adaptive procedure was deserialized from a file.
693
    bool isDeserialized() const 
694
    {
695
      return deserialized;
696
697
    }

698
699
    inline void setIsDeserialized(bool b) 
    {
700
701
702
703
704
705
      deserialized = b;
    }

    inline void setRosenbrockMode(bool b)
    {
      rosenbrockMode = b;
706
707
    }

708
    /// Creates new scalContents with the given size.
709
710
    void setScalContents(int newSize);

711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
    /** \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
726
    void serialize(std::ostream& out);
727

Thomas Witkowski's avatar
Thomas Witkowski committed
728
    void deserialize(std::istream& in);
729
730

  protected:
731
    /// Name.
Thomas Witkowski's avatar
Thomas Witkowski committed
732
    std::string name;
733

734
    /// Current space iteration
735
736
737
738
739
740
741
742
    int spaceIteration;

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

743
    /// Current timestep iteration
744
745
    int timestepIteration;

746
    /// Maximal number of iterations for choosing a timestep
747
748
    int maxTimestepIteration;

749
    /// Current time iteration
750
751
    int timeIteration;

752
    /// Maximal number of time iterations
753
754
    int maxTimeIteration;

755
    /// Actual time, end of time interval for current time step
756
757
    double time;

758
    /// Initial time
759
760
    double startTime;

761
    /// Final time
762
763
    double endTime;

764
    ///Time step size to be used
765
766
    double timestep;

767
768
769
    /// Last processed time step size of finished iteration
    double lastProcessedTimestep;

770
    /// Minimal step size
771
772
    double minTimestep;

773
    /// Maximal step size
774
    double maxTimestep;
775
    
776
    /// Number of current time step
777
    int timestepNumber;
Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
778
779
780
781
782
783
784

    /** \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;
785
  
786
    /// number of iterations needed of linear or nonlinear solver
787
788
    int solverIterations;

789
    /// maximal number of iterations needed of linear or nonlinear solver
790
791
    int maxSolverIterations;

792
    ///
793
794
    double solverTolerance;

795
    ///
796
797
    double solverResidual;

798
    /// Scalar adapt infos.
Thomas Witkowski's avatar
Thomas Witkowski committed
799
    std::vector<ScalContent*> scalContents;
800

801
    /// Is true, if the adaptive procedure was deserialized from a file.
802
803
804
805
    bool deserialized;

    /// Is true, if the time adaption is controlled by a Rosenbrock Method.
    bool rosenbrockMode;
Praetorius, Simon's avatar
Praetorius, Simon committed
806
807
808
    
    /// overall time error estimate
    double est_t;
809
810
811
812
813
  };

}

#endif //  AMDIS_ADAPTINFO_H