AdaptInfo.hpp 17.8 KB
Newer Older
1
2
3
#pragma once

// std c++ headers
4
#include <algorithm>
5
#include <cmath>
6
#include <forward_list>
7
#include <limits>
8
#include <map>
9
#include <string>
10
#include <utility>
11
12
13
#include <vector>

// AMDiS includes
14
#include <amdis/Output.hpp>
15
#include <amdis/common/ConceptsBase.hpp>
16
#include <amdis/common/Math.hpp>
17
#include <amdis/utility/TreePath.hpp>
18
19
20
21
22
23
24
25

namespace AMDiS
{

  /**
   * \ingroup Adaption
   *
   * \brief
26
   * Holds adapt parameters and infos about the problem.
27
28
29
   */
  class AdaptInfo
  {
30
31
32
  public:
    using Key = std::string;

33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
  protected:
    /** \brief
     * Stores adapt infos for a scalar problem or for one component of a
     * vector valued problem.
     */
    class ScalContent
    {
    public:
      /// Constructor.
      explicit ScalContent(std::string prefix);

      /// Sum of all error estimates
      double est_sum = 0.0;

      /// Sum of all time error estimates
      double est_t_sum = 0.0;

      /// maximal local error estimate.
      double est_max = 0.0;

      /// Maximum of all time error estimates
      double est_t_max = 0.0;

      /// factors to combine max and integral time estimate
      double fac_max = 0.0, fac_sum = 1.0;

      /// Tolerance for the (absolute or relative) error
      double spaceTolerance = 0.0;

      /// Time tolerance.
      double timeTolerance = 0.0;

      /// Relative time tolerance
      double timeRelativeTolerance = 0.0;

      /// Lower bound for the time error.
      double timeErrLow = 0.0;

      /// true if coarsening is allowed, false otherwise.
      int coarsenAllowed = 0;

      /// true if refinement is allowed, false otherwise.
      int refinementAllowed = 1;
    };

  public:
    /// Constructor.
80
    explicit AdaptInfo(std::string name);
81
82

    /// Destructor.
83
    virtual ~AdaptInfo() {}
84
85
86
87
88
89
90

    /// Resets all variables to zero (or something equivalent)
    void reset();

    /// Returns whether space tolerance is reached.
    virtual bool spaceToleranceReached() const
    {
91
92
      for (auto const& scalContent : scalContents)      {
        if (!(scalContent.second.est_sum < scalContent.second.spaceTolerance))
93
94
95
96
97
98
          return false;
      }

      return true;
    }

99
    /// Returns whether space tolerance of component associated with key is reached.
100
    virtual bool spaceToleranceReached(Key key) const
101
    {
102
      if (!(getScalContent(key).est_sum < getScalContent(key).spaceTolerance))
103
104
105
106
107
        return false;
      else
        return true;
    }

108
    template <class TP, REQUIRES( Concepts::PreTreePath<TP> )>
109
    bool spaceToleranceReached(const TP& tp) const
110
111
112
113
    {
      return spaceToleranceReached(to_string(tp));
    }

114
115
116
    /// Returns whether time tolerance is reached.
    virtual bool timeToleranceReached() const
    {
117
118
      for (auto const& scalContent : scalContents)
        if (!(getTimeEstCombined(scalContent.first) < scalContent.second.timeTolerance))
119
120
121
122
123
          return false;

      return true;
    }

124
    /// Returns whether time tolerance of component associated with key is reached.
125
    virtual bool timeToleranceReached(Key key) const
126
    {
127
      if (!(getTimeEstCombined(key) < getScalContent(key).timeTolerance))
128
129
130
131
132
        return false;
      else
        return true;
    }

133
    template <class TP, REQUIRES( Concepts::PreTreePath<TP> )>
134
    bool timeToleranceReached(const TP& tp) const
135
136
137
138
    {
      return timeToleranceReached(to_string(tp));
    }

139
140
141
    /// Returns whether time error is under its lower bound.
    virtual bool timeErrorLow() const
    {
142
143
      for (auto const& scalContent : scalContents)
        if (!(getTimeEstCombined(scalContent.first) < scalContent.second.timeErrLow))
144
145
146
147
          return false;

      return true;
    }
148

149
150
    /// Returns the time estimation as a combination
    /// of maximal and integral time error
151
    double getTimeEstCombined(Key key) const
152
153
    {
      return
154
155
        getScalContent(key).est_t_max * getScalContent(key).fac_max +
        getScalContent(key).est_t_sum * getScalContent(key).fac_sum;
156
157
    }

158
    template <class TP, REQUIRES( Concepts::PreTreePath<TP> )>
159
    double getTimeEstCombined(const TP& tp) const
160
161
162
    {
      return getTimeEstCombined(to_string(tp));
    }
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287

    /// Print debug information about time error and its bound.
    void printTimeErrorLowInfo() const;

    /// Returns \ref spaceIteration.
    int getSpaceIteration() const
    {
      return spaceIteration;
    }

    /// Sets \ref spaceIteration.
    void setSpaceIteration(int it)
    {
      spaceIteration = it;
    }

    /// Returns \ref maxSpaceIteration.
    int getMaxSpaceIteration() const
    {
      return maxSpaceIteration;
    }

    /// Sets \ref maxSpaceIteration.
    void setMaxSpaceIteration(int it)
    {
      maxSpaceIteration = it;
    }

    /// Increments \ref spaceIteration by 1;
    void incSpaceIteration()
    {
      spaceIteration++;
    }

    /// Sets \ref timestepIteration.
    void setTimestepIteration(int it)
    {
      timestepIteration = it;
    }

    /// Returns \ref timestepIteration.
    int getTimestepIteration() const
    {
      return timestepIteration;
    }

    /// Increments \ref timestepIteration by 1;
    void incTimestepIteration()
    {
      timestepIteration++;
    }

    /// Returns \ref maxTimestepIteration.
    int getMaxTimestepIteration() const
    {
      return maxTimestepIteration;
    }

    /// Sets \ref maxTimestepIteration.
    void setMaxTimestepIteration(int it)
    {
      maxTimestepIteration = it;
    }

    /// Sets \ref timeIteration.
    void setTimeIteration(int it)
    {
      timeIteration = it;
    }

    /// Returns \ref timeIteration.
    int getTimeIteration() const
    {
      return timeIteration;
    }

    /// Increments \ref timesIteration by 1;
    void incTimeIteration()
    {
      timeIteration++;
    }

    /// Returns \ref maxTimeIteration.
    int getMaxTimeIteration() const
    {
      return maxTimeIteration;
    }

    /// Sets \ref maxTimeIteration.
    void setMaxTimeIteration(int it)
    {
      maxTimeIteration = it;
    }

    /// Returns \ref timestepNumber.
    int getTimestepNumber() const
    {
      return timestepNumber;
    }

    /// Sets \ref timestepNumber.
    void setTimestepNumber(int num)
    {
      timestepNumber = std::min(nTimesteps, num);
    }

    /// Returns \ref nTimesteps.
    int getNumberOfTimesteps() const
    {
      return nTimesteps;
    }

    /// Sets \ref nTimesteps.
    void setNumberOfTimesteps(int num)
    {
      nTimesteps = std::max(0, num);
    }

    /// Increments \ref timestepNumber by 1;
    void incTimestepNumber()
    {
      timestepNumber++;
    }

    /// Sets \ref est_sum.
288
    void setEstSum(double e, Key key)
289
    {
290
      getScalContent(key).est_sum = e;
291
292
    }

293
    template <class TP, REQUIRES( Concepts::PreTreePath<TP> )>
294
    void setEstSum(double e, const TP& tp)
295
296
297
298
    {
      setEstSum(e, to_string(tp));
    }

299
    /// Sets \ref est_max.
300
    void setEstMax(double e, Key key)
301
    {
302
      getScalContent(key).est_max = e;
303
304
    }

305
    template <class TP, REQUIRES( Concepts::PreTreePath<TP> )>
306
    void setEstMax(double e, const TP& tp)
307
308
309
310
    {
      setEstMax(e, to_string(tp));
    }

311
    /// Sets \ref est_max.
312
    void setTimeEstMax(double e, Key key)
313
    {
314
      getScalContent(key).est_t_max = e;
315
316
    }

317
    template <class TP, REQUIRES( Concepts::PreTreePath<TP> )>
318
    void setTimeEstMax(double e, const TP& tp)
319
320
321
322
    {
      setTimeEstMax(e, to_string(tp));
    }

323
    /// Sets \ref est_t_sum.
324
    void setTimeEstSum(double e, Key key)
325
    {
326
      getScalContent(key).est_t_sum = e;
327
328
    }

329
    template <class TP, REQUIRES( Concepts::PreTreePath<TP> )>
330
    void setTimeEstSum(double e, const TP& tp)
331
332
333
334
    {
      setTimeEstSum(e, to_string(tp));
    }

335
    /// Returns \ref est_sum.
336
    double getEstSum(Key key) const
337
    {
338
      return getScalContent(key).est_sum;
339
340
    }

341
    template <class TP, REQUIRES( Concepts::PreTreePath<TP> )>
342
    double getEstSum(const TP& tp)
343
344
345
346
    {
      return getEstSum(to_string(tp));
    }

347
    /// Returns \ref est_t_sum.
348
    double getEstTSum(Key key) const
349
    {
350
      return getScalContent(key).est_t_sum;
351
352
    }

353
    template <class TP, REQUIRES( Concepts::PreTreePath<TP> )>
354
    double getEstTSum(const TP& tp)
355
356
357
358
    {
      return getEstTSum(to_string(tp));
    }

359
    /// Returns \ref est_max.
360
    double getEstMax(Key key) const
361
    {
362
      return getScalContent(key).est_max;
363
364
    }

365
    template <class TP, REQUIRES( Concepts::PreTreePath<TP> )>
366
    double getEstMax(const TP& tp)
367
368
369
370
    {
      return getEstMax(to_string(tp));
    }

371
    /// Returns \ref est_max.
372
    double getTimeEstMax(Key key) const
373
    {
374
      return getScalContent(key).est_t_max;
375
376
    }

377
    template <class TP, REQUIRES( Concepts::PreTreePath<TP> )>
378
    double getTimeEstmax(const TP& tp)
379
380
381
382
    {
      return getTimeEstMax(to_string(tp));
    }

383
    /// Returns \ref est_t_sum.
384
    double getTimeEstSum(Key key) const
385
    {
386
      return getScalContent(key).est_t_sum;
387
388
    }

389
    template <class TP, REQUIRES( Concepts::PreTreePath<TP> )>
390
    double getTimeEstSum(const TP& tp)
391
392
393
394
    {
      return getTimeEstSum(to_string(tp));
    }

395
396
397
398
399
400
401
402
403
404
405
406
    /// Returns \ref est_t the estimated overall time error
    double getTimeEst() const
    {
      return est_t;
    }

    void setTimeEst(double value)
    {
      est_t = value;
    }

    /// Returns \ref spaceTolerance.
407
    double getSpaceTolerance(Key key) const
408
    {
409
      return getScalContent(key).spaceTolerance;
410
411
    }

412
    template <class TP, REQUIRES( Concepts::PreTreePath<TP> )>
413
    double getSpaceTolerance(const TP& tp)
414
415
416
417
    {
      return getSpaceTolerance(to_string(tp));
    }

418
    /// Sets \ref spaceTolerance.
419
    void setSpaceTolerance(Key key, double tol)
420
    {
421
      getScalContent(key).spaceTolerance = tol;
422
423
    }

424
    template <class TP, REQUIRES( Concepts::PreTreePath<TP> )>
425
    void setSpaceTolerance(const TP& tp, double tol)
426
427
428
429
    {
      return setSpaceTolerance(to_string(tp), tol);
    }

430
    /// Returns \ref timeTolerance.
431
    double getTimeTolerance(Key key) const
432
    {
433
      return getScalContent(key).timeTolerance;
434
435
    }

436
    template <class TP, REQUIRES( Concepts::PreTreePath<TP> )>
437
    double getTimeTolerance(const TP& tp)
438
439
440
441
    {
      return getTimeTolerance(to_string(tp));
    }

442
    /// Returns \ref timeRelativeTolerance.
443
    double getTimeRelativeTolerance(Key key) const
444
    {
445
      return getScalContent(key).timeRelativeTolerance;
446
447
    }

448
    template <class TP, REQUIRES( Concepts::PreTreePath<TP> )>
449
    double getTimeRelativeTolerance(const TP& tp)
450
451
452
453
    {
      return getTimeRelativeTolerance(to_string(tp));
    }

454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
    /// Sets \ref time
    double setTime(double t)
    {
      time = t;
      if (time > endTime)
        time = endTime;
      if (time < startTime)
        time = startTime;

      return time;
    }

    /// Gets \ref time
    double getTime() const
    {
      return time;
    }

    /// Gets \ref &time
    double* getTimePtr()
    {
      return &time;
    }

    /// Sets \ref timestep
    double setTimestep(double t)
    {
      timestep = t;
      if (timestep > maxTimestep)
        timestep = maxTimestep;
      if (timestep < minTimestep)
        timestep = minTimestep;
      if (time + timestep > endTime)
        timestep = endTime - time;

      return timestep;
    }
    /// Gets \ref timestep
    double getTimestep() const
    {
      return timestep;
    }

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

    double getLastProcessedTimestep() const
    {
      return lastProcessedTimestep;
    }

    /// Returns true, if the end time is reached and no more timestep
    /// computations must be done.
    bool reachedEndTime() const
    {
      if (nTimesteps > 0)
        return !(timestepNumber < nTimesteps);

514
      return !(std::abs(time - endTime) > std::numeric_limits<double>::epsilon());
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
    }


    /// Sets \ref minTimestep
    void setMinTimestep(double t)
    {
      minTimestep = t;
    }

    /// Gets \ref minTimestep
    double getMinTimestep() const
    {
      return minTimestep;
    }

    /// Sets \ref maxTimestep
    void setMaxTimestep(double t)
    {
      maxTimestep = t;
    }

    /// Gets \ref maxTimestep
    double getMaxTimestep() const
    {
      return maxTimestep;
    }

    /// Gets \ref &timestep
    double* getTimestepPtr()
    {
      return &timestep;
    }

    /// Sets \ref startTime = time
    void setStartTime(double time)
    {
      startTime = time;
    }

    /// Sets \ref endTime = time
    void setEndTime(double time)
    {
      endTime = time;
    }

    /// Returns \ref startTime
    double getStartTime() const
    {
      return startTime;
    }

    /// Returns \ref endTime
    double getEndTime() const
    {
      return endTime;
    }

    /// Returns \ref timeErrLow.
573
    double getTimeErrLow(Key key) const
574
    {
575
      return getScalContent(key).timeErrLow;
576
577
    }

578
    template <class TP, REQUIRES( Concepts::PreTreePath<TP> )>
579
    double getTimeErrLow(const TP& tp)
580
581
582
583
    {
      return getTimeErrLow(to_string(tp));
    }

584
    /// Returns whether coarsening is allowed or not.
585
    bool isCoarseningAllowed(Key key) const
586
    {
587
      return (getScalContent(key).coarsenAllowed == 1);
588
589
    }

590
    template <class TP, REQUIRES( Concepts::PreTreePath<TP> )>
591
    bool isCoarseningAllowed(const TP& tp)
592
593
594
595
    {
      return isCoarseningAllowed(to_string(tp));
    }

596
    /// Returns whether coarsening is allowed or not.
597
    bool isRefinementAllowed(Key key) const
598
    {
599
      return (getScalContent(key).refinementAllowed == 1);
600
601
    }

602
    template <class TP, REQUIRES( Concepts::PreTreePath<TP> )>
603
    bool isRefinementAllowed(const TP& tp)
604
605
606
607
    {
      return isRefinementAllowed(to_string(tp));
    }

608
    ///
609
    void allowRefinement(bool allow, Key key)
610
    {
611
      getScalContent(key).refinementAllowed = allow;
612
613
    }

614
    template <class TP, REQUIRES( Concepts::PreTreePath<TP> )>
615
    void allowRefinement(bool allow, const TP& tp)
616
617
618
619
    {
      return allowRefinement(allow, to_string(tp));
    }

620
    ///
621
    void allowCoarsening(bool allow, Key key)
622
    {
623
      getScalContent(key).coarsenAllowed = allow;
624
625
    }

626
    template <class TP, REQUIRES( Concepts::PreTreePath<TP> )>
627
    void allowCoarsening(bool allow, const TP& tp)
628
629
630
631
    {
      return allowCoarsening(allow, to_string(tp));
    }

632
633
    int getSize() const
    {
634
      return int(scalContents.size());
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
    }

    // TODO: remove from AdaptInfo
    bool getRosenbrockMode() const
    {
      return rosenbrockMode;
    }

    void setSolverIterations(int it)
    {
      solverIterations = it;
    }

    int getSolverIterations() const
    {
      return solverIterations;
    }

    void setMaxSolverIterations(int it)
    {
      maxSolverIterations = it;
    }

    int getMaxSolverIterations() const
    {
      return maxSolverIterations;
    }

    void setSolverTolerance(double tol)
    {
      solverTolerance = tol;
    }

    double getSolverTolerance() const
    {
      return solverTolerance;
    }

    void setSolverResidual(double res)
    {
      solverResidual = res;
    }

    double getSolverResidual() const
    {
      return solverResidual;
    }

    void setGlobalTimeTolerance(double tol)
    {
      globalTimeTolerance = tol;
    }

    double getGlobalTimeTolerance() const
    {
      return globalTimeTolerance;
    }

    // TODO: remove from AdaptInfo
    void setRosenbrockMode(bool b)
    {
      rosenbrockMode = b;
    }

    /** \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;
    }

714
  private:
715
    ScalContent& getScalContent(Key key) const
716
717
718
719
720
    {
      auto result = scalContents.emplace(std::piecewise_construct, std::forward_as_tuple(key), std::forward_as_tuple(name + "[" + key + "]") );
      return result.first->second;
    }

721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
  protected:
    /// Name.
    std::string name;

    /// Current space iteration
    int spaceIteration = -1;

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

    /// Current timestep iteration
    int timestepIteration = 0;

    /// Maximal number of iterations for choosing a timestep
    int maxTimestepIteration = 30;

    /// Current time iteration
    int timeIteration = 0;

    /// Maximal number of time iterations
    int maxTimeIteration = 30;

    /// Actual time, end of time interval for current time step
    double time = 0.0;

    /// Initial time
    double startTime = 0.0;

    /// Final time
    double endTime = 1.0;

755
    /// Time step size to be used
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
    double timestep = 0.0;

    /// Last processed time step size of finished iteration
    double lastProcessedTimestep = 0.0;

    /// Minimal step size
    double minTimestep = 0.0;

    /// Maximal step size
    double maxTimestep = 1.0;

    /// Number of current time step
    int timestepNumber = 0;

    /** \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 = 0;

    /// number of iterations needed of linear or nonlinear solver
    int solverIterations = 0;

    /// maximal number of iterations needed of linear or nonlinear solver
    int maxSolverIterations = 0;

    ///
    double solverTolerance = 1.e-8;

    ///
    double solverResidual = 0.0;

    /// tolerance for the overall time error
    double globalTimeTolerance = 1.0;

792
    /// Scalar adapt infos
793
    mutable std::map<Key, ScalContent> scalContents;
794
795
796
797
798
799
800
801
802
803
804
805

    /// Is true, if the adaptive procedure was deserialized from a file. TODO: remove deserialization
    bool deserialized = false;

    /// Is true, if the time adaption is controlled by a Rosenbrock Method.
    bool rosenbrockMode = false;

    /// overall time error estimate
    double est_t = 0.0;
  };

} // end namespace AMDiS