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 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
      {
	double timeTheta2 = 0.3;

66 67 68 69
	// TODO: obsolete parameters timeTheta2, relTimeErr, relSpaceErr

	GET_PARAMETER(0, prefix + "->tolerance", "%f", &spaceTolerance);
	GET_PARAMETER(0, prefix + "->time tolerance", "%f", &timeTolerance);
70 71 72 73 74
	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);
75 76
	GET_PARAMETER(0, prefix + "->sum factor", "%f", &fac_sum);
	GET_PARAMETER(0, prefix + "->max factor", "%f", &fac_max);
77

78
	timeErrLow = timeTolerance * timeTheta2;
79
      }
80

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

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

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

90
      /// Maximum of all time error estimates
91
      double est_t_max;
92 93 94
      
      /// factors to combine max and integral time estimate
      double fac_max, fac_sum;
95

96
      /// Tolerance for the (absolute or relative) error
97 98
      double spaceTolerance;

99
      /// Time tolerance.
100 101
      double timeTolerance;

102
      /// Lower bound for the time error.
103 104
      double timeErrLow;

105
      /// true if coarsening is allowed, false otherwise.
106 107
      int coarsenAllowed;

108
      /// true if refinement is allowed, false otherwise.
109 110 111 112 113 114 115
      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
       */
116
      int refineBisections;
117 118 119 120 121 122

      /** \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
       */                          
123
      int coarseBisections;    
124 125 126
    };

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

166 167 168 169 170 171
      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
172
	  scalContents[i] = new ScalContent(name + std::string(number)); 
173 174
	}
      }
175
    }
176

177
    /// Destructor.
Thomas Witkowski's avatar
Thomas Witkowski committed
178 179 180
    virtual ~AdaptInfo() 
    {
      for (unsigned int i = 0;  i < scalContents.size(); i++)
181
	delete scalContents[i];
182
    }
183 184 185 186 187 188 189 190 191 192 193

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

      GET_PARAMETER(0, name + "->timestep", "%f", &timestep);
196
      lastProcessedTimestep=timestep;
197
    }
198

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

210
      return true;
211
    }
212

213
    /// Returns whether space tolerance of component i is reached.
214 215
    virtual bool spaceToleranceReached(int i) 
    {
216
      if (!(scalContents[i]->est_sum < scalContents[i]->spaceTolerance))
217
	return false;
218
      else
219
	return true;
220
    }
221

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

229
      return true;
230
    }
231

232
    /// Returns whether time tolerance of component i is reached.
233 234
    virtual bool timeToleranceReached(int i) 
    {
235
      if (!(getTimeEstCombined(i) < scalContents[i]->timeTolerance))
236
	return false;
237
      else
238
	return true;
239
    }
240

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

248
      return true;
249
    }
250 251 252 253 254 255 256 257
    /// Returns the time estimation as a combination 
    /// of maximal and integral time error 
    double getTimeEstCombined(unsigned i) const 
    { 
      return scalContents[i]->est_t_max*scalContents[i]->fac_max
	+scalContents[i]->est_t_sum*scalContents[i]->fac_sum; 
    }

258

259
    /// Print debug information about time error and its bound.
260 261
    void printTimeErrorLowInfo() 
    {
262 263 264 265
      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 
266 267
		  << "    Time error low bound = " << scalContents[i]->timeErrLow  
		  << "    Time error high bound = " << scalContents[i]->timeTolerance << "\n";
268
      }
269
    }
270

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

423
      return scalContents[index]->est_sum; 
424
    }
425

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

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

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

440
      return scalContents[index]->est_max; 
441
    }
442

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

449
    /// Returns \ref est_t_sum.
450 451
    inline double getTimeEstSum(int index) 
    { 
452
      return scalContents[index]->est_t_sum; 
453
    }
454

455
    /// Returns \ref spaceTolerance.
456 457
    inline double getSpaceTolerance(int index) 
    { 
458
      return scalContents[index]->spaceTolerance; 
459
    }  
460

461
    /// Sets \ref spaceTolerance.
462 463
    inline void setSpaceTolerance(int index, double tol) 
    { 
464
      scalContents[index]->spaceTolerance = tol; 
465
    }  
466

467
    /// Returns \ref timeTolerance.
468 469
    inline double getTimeTolerance(int index) 
    { 
470
      return scalContents[index]->timeTolerance; 
471
    }  
472

473
    /// Sets \ref time
474 475
    inline double setTime(double t) 
    { 
476
      time = t; 
477 478 479 480 481
      if (time > endTime) 
	time = endTime;
      if (time < startTime) 
	time = startTime;

482
      return time;
483
    }
484

485
    /// Gets \ref time
486 487
    inline double getTime() 
    { 
488
      return time; 
489
    }  
490

491
    /// Gets \ref &time
492 493
    inline double *getTimePtr() 
    { 
494
      return &time; 
495
    }  
496

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

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

    inline double getLastProcessedTimestep(){
	return lastProcessedTimestep;
    } 
523

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

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

536

537
    /// Sets \ref minTimestep
538 539
    inline void setMinTimestep(double t) 
    { 
540
      minTimestep = t; 
541
    }
542

543
    /// Gets \ref minTimestep
544 545
    inline double getMinTimestep() 
    { 
546
      return minTimestep; 
547
    }  
548

549
    /// Sets \ref maxTimestep
550 551
    inline void setMaxTimestep(double t) 
    { 
552
      maxTimestep = t; 
553
    }
554

555
    /// Gets \ref maxTimestep
556 557
    inline double getMaxTimestep() 
    { 
558
      return maxTimestep; 
559
    }  
560
    
561
    /// Gets \ref &timestep
562 563
    inline double *getTimestepPtr() 
    { 
564
      return &timestep; 
565
    }  
566

567
    /// Sets \ref startTime = time
568 569
    inline void setStartTime(double time) 
    { 
570
      startTime = time; 
571
    }
572

573
    /// Sets \ref endTime = time
574 575
    inline void setEndTime(double time) 
    { 
576
      endTime = time; 
577
    }
578

579
    /// Returns \ref startTime
580 581
    inline double getStartTime() 
    { 
582
      return startTime; 
583
    }
584

585
    /// Returns \ref endTime
586 587
    inline double getEndTime() 
    { 
588
      return endTime; 
589
    }
590

591
    /// Returns \ref timeErrLow.
592 593
    inline double getTimeErrLow(int index) 
    { 
594
      return scalContents[index]->timeErrLow; 
595
    }  
596

597
    /// Returns whether coarsening is allowed or not.
598 599
    inline bool isCoarseningAllowed(int index) 
    {
600
      return (scalContents[index]->coarsenAllowed == 1);
601
    }
602

603
    /// Returns whether coarsening is allowed or not.
604 605
    inline bool isRefinementAllowed(int index) 
    {
606
      return (scalContents[index]->refinementAllowed == 1);
607
    }
608

609
    ///
610 611
    inline void allowRefinement(bool allow, int index) 
    {
612
      scalContents[index]->refinementAllowed = allow;
613
    }
614

615
    ///
616 617
    inline void allowCoarsening(bool allow, int index) 
    {
618
      scalContents[index]->coarsenAllowed = allow;
619
    }
620

621
    /// Returns \ref refineBisections
622 623
    inline const int getRefineBisections(int index) const 
    {
624
      return scalContents[index]->refineBisections;
625
    }
626

627
    /// Returns \ref coarseBisections
628 629
    inline const int getCoarseBisections(int index) const 
    {
630
      return scalContents[index]->coarseBisections;
631
    }    
632

633 634
    inline int getSize() 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
635
      return scalContents.size(); 
636
    }
637

638 639 640 641 642
    inline bool getRosenbrockMode()
    {
      return rosenbrockMode;
    }

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

683
    /// Returns true, if the adaptive procedure was deserialized from a file.
684 685
    const bool isDeserialized() const 
    {
686
      return deserialized;
687 688
    }

689 690
    inline void setIsDeserialized(bool b) 
    {
691 692 693 694 695 696
      deserialized = b;
    }

    inline void setRosenbrockMode(bool b)
    {
      rosenbrockMode = b;
697 698
    }

699
    /// Creates new scalContents with the given size.
700 701
    void setScalContents(int newSize);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
719
    void deserialize(std::istream& in);
720 721

  protected: