AdaptInfo.h 18.4 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.
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
      Parameters::get( name_ + "->start time", startTime);
149
      time = startTime;
150
151
152
153
154
155
156
157
158
159
160
161
162
      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);

      char number[5];
      for (int i = 0; i < size; i++) {
	sprintf(number, "[%d]", i);
	scalContents[i] = new ScalContent(name + std::string(number)); 	
163
      }
164
    }
165

166
    /// Destructor.
Thomas Witkowski's avatar
Thomas Witkowski committed
167
168
169
    virtual ~AdaptInfo() 
    {
      for (unsigned int i = 0;  i < scalContents.size(); i++)
170
	delete scalContents[i];
171
    }
172
173
174
175
176
177
178
179
180
181
182

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

184
      Parameters::get(name + "->timestep", timestep);
185
      lastProcessedTimestep=timestep;
186
    }
187

188
    /// Returns whether space tolerance is reached.
189
190
    virtual bool spaceToleranceReached() 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
191
      for (unsigned int i = 0; i < scalContents.size(); i++) {
192
	if (!(scalContents[i]->est_sum < scalContents[i]->spaceTolerance))
193
	  return false;
194
      }
195

196
      return true;
197
    }
198

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

208
    /// Returns whether time tolerance is reached.
209
210
    virtual bool timeToleranceReached() 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
211
      for (unsigned int i = 0; i < scalContents.size(); i++)
212
	if (!(getTimeEstCombined(i) < scalContents[i]->timeTolerance))
213
	  return false;
214

215
      return true;
216
    }
217

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

227
    /// Returns whether time error is under its lower bound.
228
229
    virtual bool timeErrorLow() 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
230
      for (unsigned int i = 0; i < scalContents.size(); i++)
231
	if (!(getTimeEstCombined(i) < scalContents[i]->timeErrLow))
232
	  return false;
233

234
      return true;
235
    }
236
237
238
239
    /// 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
240
241
242
      return 
	scalContents[i]->est_t_max * scalContents[i]->fac_max +
	scalContents[i]->est_t_sum * scalContents[i]->fac_sum; 
243
244
    }

245

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

263
    /// Returns \ref spaceIteration.
264
265
    inline int getSpaceIteration() 
    { 
266
      return spaceIteration; 
267
    }
268

269
    /// Sets \ref spaceIteration.
270
271
    inline void setSpaceIteration(int it) 
    { 
272
      spaceIteration = it; 
273
    }
274
  
275
    /// Returns \ref maxSpaceIteration.
276
277
    inline int getMaxSpaceIteration() 
    { 
278
      return maxSpaceIteration;
279
    }
280

281
    /// Sets \ref maxSpaceIteration.
282
283
    inline void setMaxSpaceIteration(int it) 
    { 
284
      maxSpaceIteration = it; 
285
    }
286
  
287
    /// Increments \ref spaceIteration by 1;
288
289
    inline void incSpaceIteration() 
    { 
290
      spaceIteration++; 
291
    }
292

293
    /// Sets \ref timestepIteration.
294
295
    inline void setTimestepIteration(int it) 
    { 
296
      timestepIteration = it; 
297
    }
298
  
299
    /// Returns \ref timestepIteration.
300
301
    inline int getTimestepIteration() 
    { 
302
      return timestepIteration; 
303
    }
304

305
    /// Increments \ref timestepIteration by 1;
306
307
    inline void incTimestepIteration() 
    { 
308
      timestepIteration++; 
309
    }
310

311
    /// Returns \ref maxTimestepIteration.
312
313
    inline int getMaxTimestepIteration() 
    { 
314
      return maxTimestepIteration; 
315
    }
316

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

335
    /// Increments \ref timesIteration by 1;
336
337
    inline void incTimeIteration() 
    { 
338
      timeIteration++; 
339
    }
340

341
    /// Returns \ref maxTimeIteration.
342
343
    inline int getMaxTimeIteration() 
    { 
344
      return maxTimeIteration; 
345
    }
346

347
    /// Sets \ref maxTimeIteration.
348
349
    inline void setMaxTimeIteration(int it) 
    { 
350
      maxTimeIteration = it; 
351
    }
352
  
353
    /// Returns \ref timestepNumber.
354
355
    inline int getTimestepNumber() 
    { 
356
      return timestepNumber; 
357
    }
358

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

371
372
373
    /// Sets \ref nTimesteps.
    inline void setNumberOfTimesteps(int num) 
    {
374
      nTimesteps = std::max(0, num);
375
376
    }

377
    /// Increments \ref timestepNumber by 1;
378
379
    inline void incTimestepNumber() 
    { 
380
      timestepNumber++; 
381
    }
382

383
    /// Sets \ref est_sum.
384
385
    inline void setEstSum(double e, int index) 
    {
386
      scalContents[index]->est_sum = e;
387
    }
388

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

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

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

407
    /// Returns \ref est_sum.
408
409
    inline double getEstSum(int index) 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
410
411
412
413
414
      FUNCNAME("AdaptInfo::getEstSum()");

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

415
      return scalContents[index]->est_sum; 
416
    }
417

418
    /// Returns \ref est_t_sum.
419
420
    inline double getEstTSum(int index) 
    { 
421
      return scalContents[index]->est_t_sum; 
422
    }
423

424
    /// Returns \ref est_max.
425
426
    inline double getEstMax(int index) 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
427
428
429
430
431
      FUNCNAME("AdaptInfo::getEstSum()");

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

432
      return scalContents[index]->est_max; 
433
    }
434

435
    /// Returns \ref est_max.
436
437
    inline double getTimeEstMax(int index) 
    { 
438
      return scalContents[index]->est_t_max; 
439
    }
440

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

447
    /// Returns \ref spaceTolerance.
448
449
    inline double getSpaceTolerance(int index) 
    { 
450
      return scalContents[index]->spaceTolerance; 
451
    }  
452

453
    /// Sets \ref spaceTolerance.
454
455
    inline void setSpaceTolerance(int index, double tol) 
    { 
456
      scalContents[index]->spaceTolerance = tol; 
457
    }  
458

459
    /// Returns \ref timeTolerance.
460
461
    inline double getTimeTolerance(int index) 
    { 
462
      return scalContents[index]->timeTolerance; 
463
    }  
464

465
    /// Sets \ref time
466
467
    inline double setTime(double t) 
    { 
468
      time = t; 
469
470
471
472
473
      if (time > endTime) 
	time = endTime;
      if (time < startTime) 
	time = startTime;

474
      return time;
475
    }
476

477
    /// Gets \ref time
478
479
    inline double getTime() 
    { 
480
      return time; 
481
    }  
482

483
    /// Gets \ref &time
484
485
    inline double *getTimePtr() 
    { 
486
      return &time; 
487
    }  
488

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

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

    inline double getLastProcessedTimestep(){
	return lastProcessedTimestep;
    } 
515

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

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

528

529
    /// Sets \ref minTimestep
530
531
    inline void setMinTimestep(double t) 
    { 
532
      minTimestep = t; 
533
    }
534

535
    /// Gets \ref minTimestep
536
537
    inline double getMinTimestep() 
    { 
538
      return minTimestep; 
539
    }  
540

541
    /// Sets \ref maxTimestep
542
543
    inline void setMaxTimestep(double t) 
    { 
544
      maxTimestep = t; 
545
    }
546

547
    /// Gets \ref maxTimestep
548
549
    inline double getMaxTimestep() 
    { 
550
      return maxTimestep; 
551
    }  
552
    
553
    /// Gets \ref &timestep
554
555
    inline double *getTimestepPtr() 
    { 
556
      return &timestep; 
557
    }  
558

559
    /// Sets \ref startTime = time
560
561
    inline void setStartTime(double time) 
    { 
562
      startTime = time; 
563
    }
564

565
    /// Sets \ref endTime = time
566
567
    inline void setEndTime(double time) 
    { 
568
      endTime = time; 
569
    }
570

571
    /// Returns \ref startTime
572
573
    inline double getStartTime() 
    { 
574
      return startTime; 
575
    }
576

577
    /// Returns \ref endTime
578
579
    inline double getEndTime() 
    { 
580
      return endTime; 
581
    }
582

583
    /// Returns \ref timeErrLow.
584
585
    inline double getTimeErrLow(int index) 
    { 
586
      return scalContents[index]->timeErrLow; 
587
    }  
588

589
    /// Returns whether coarsening is allowed or not.
590
591
    inline bool isCoarseningAllowed(int index) 
    {
592
      return (scalContents[index]->coarsenAllowed == 1);
593
    }
594

595
    /// Returns whether coarsening is allowed or not.
596
597
    inline bool isRefinementAllowed(int index) 
    {
598
      return (scalContents[index]->refinementAllowed == 1);
599
    }
600

601
    ///
602
603
    inline void allowRefinement(bool allow, int index) 
    {
604
      scalContents[index]->refinementAllowed = allow;
605
    }
606

607
    ///
608
609
    inline void allowCoarsening(bool allow, int index) 
    {
610
      scalContents[index]->coarsenAllowed = allow;
611
    }
612

613
    /// Returns \ref refineBisections
614
615
    inline const int getRefineBisections(int index) const 
    {
616
      return scalContents[index]->refineBisections;
617
    }
618

619
    /// Returns \ref coarseBisections
620
621
    inline const int getCoarseBisections(int index) const 
    {
622
      return scalContents[index]->coarseBisections;
623
    }    
624

625
626
    inline int getSize() 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
627
      return scalContents.size(); 
628
    }
629

630
631
632
633
634
    inline bool getRosenbrockMode()
    {
      return rosenbrockMode;
    }

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

675
    /// Returns true, if the adaptive procedure was deserialized from a file.
676
677
    const bool isDeserialized() const 
    {
678
      return deserialized;
679
680
    }

681
682
    inline void setIsDeserialized(bool b) 
    {
683
684
685
686
687
688
      deserialized = b;
    }

    inline void setRosenbrockMode(bool b)
    {
      rosenbrockMode = b;
689
690
    }

691
    /// Creates new scalContents with the given size.
692
693
    void setScalContents(int newSize);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
711
    void deserialize(std::istream& in);
712
713

  protected:
714
    /// Name.
Thomas Witkowski's avatar
Thomas Witkowski committed
715
    std::string name;
716

717
    /// Current space iteration
718
719
720
721
722
723
724
725
    int spaceIteration;

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

726
    /// Current timestep iteration
727
728
    int timestepIteration;

729
    /// Maximal number of iterations for choosing a timestep
730
731
    int maxTimestepIteration;

732
    /// Current time iteration
733
734
    int timeIteration;

735
    /// Maximal number of time iterations
736
737
    int maxTimeIteration;

738
    /// Actual time, end of time interval for current time step
739
740
    double time;

741
    /// Initial time
742
743
    double startTime;

744
    /// Final time
745
746
    double endTime;

747
    ///Time step size to be used
748
749
    double timestep;

750
751
752
    /// Last processed time step size of finished iteration
    double lastProcessedTimestep;

753
    /// Minimal step size
754
755
    double minTimestep;

756
    /// Maximal step size
757
    double maxTimestep;
758
    
759
    /// Number of current time step
760
    int timestepNumber;
Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
761
762
763
764
765
766
767

    /** \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;
768
  
769
    /// number of iterations needed of linear or nonlinear solver
770
771
    int solverIterations;

772
    /// maximal number of iterations needed of linear or nonlinear solver
773
774
    int maxSolverIterations;

775
    ///
776
777
    double solverTolerance;

778
    ///
779
780
    double solverResidual;

781
    /// Scalar adapt infos.
Thomas Witkowski's avatar
Thomas Witkowski committed
782
    std::vector<ScalContent*> scalContents;
783

784
    /// Is true, if the adaptive procedure was deserialized from a file.
785
786
787
788
    bool deserialized;

    /// Is true, if the time adaption is controlled by a Rosenbrock Method.
    bool rosenbrockMode;
789
790
791
792
793
  };

}

#endif //  AMDIS_ADAPTINFO_H