AdaptInfo.h 17.9 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 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
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  crystal growth group                                                  ==
// ==                                                                        ==
// ==  Stiftung caesar                                                       ==
// ==  Ludwig-Erhard-Allee 2                                                 ==
// ==  53175 Bonn                                                            ==
// ==  germany                                                               ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  http://www.caesar.de/cg/AMDiS                                         ==
// ==                                                                        ==
// ============================================================================

/** \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:
47
      /// Constructor.
Thomas Witkowski's avatar
Thomas Witkowski committed
48
      ScalContent(std::string prefix) 
49 50 51 52 53 54 55 56 57 58
	: est_sum(0.0),
	  est_t_sum(0.0),
	  est_max(0.0),
	  est_t_max(0.0),
	  spaceTolerance(1.0),
	  timeTolerance(1.0),
	  timeErrLow(1.0),
	  coarsenAllowed(0),
	  refinementAllowed(1),
	  refineBisections(1),
59
  	  coarseBisections(1)
60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79
      {
	double totalTol = 1.0;
	double relSpaceErr = 1.0;
	double relTimeErr = 0.5;
	double timeTheta1 = 1.0;
	double timeTheta2 = 0.3;

	GET_PARAMETER(0, prefix + "->tolerance", "%f", &totalTol);
	GET_PARAMETER(0, prefix + "->rel space error", "%f", &relSpaceErr);
	GET_PARAMETER(0, prefix + "->rel time error", "%f", &relTimeErr);
	GET_PARAMETER(0, prefix + "->time theta 1", "%f", &timeTheta1);
	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);

	spaceTolerance = totalTol * relSpaceErr;
	timeTolerance = totalTol * relTimeErr * timeTheta1;
	timeErrLow = totalTol * relTimeErr * timeTheta2;
80
      }
81

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

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

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

91
      /// Maximum of all time error estimates
92 93
      double est_t_max;

94
      /// Tolerance for the (absolute or relative) error
95 96
      double spaceTolerance;

97
      /// Time tolerance.
98 99
      double timeTolerance;

100
      /// Lower bound for the time error.
101 102
      double timeErrLow;

103
      /// true if coarsening is allowed, false otherwise.
104 105
      int coarsenAllowed;

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

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

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

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

173
    /// Destructor.
174
    virtual ~AdaptInfo() {
175
      for (int i = 0;  i < scalContents.getSize(); 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
    }
192

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

204
      return true;
205
    }
206

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

216
    /// Returns whether time tolerance is reached.
217 218
    virtual bool timeToleranceReached() 
    {
219
      int size = scalContents.getSize();
220 221
      for (int i = 0; i < size; i++)
	if (!(scalContents[i]->est_t_sum < 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 (!(scalContents[i]->est_t_sum < 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() 
    {
239
      int size = scalContents.getSize();
240 241
      for (int i = 0; i < size; i++)
	if (!(scalContents[i]->est_t_sum < scalContents[i]->timeErrLow))
242
	  return false;
243

244
      return true;
245 246
    }

247
    /// Print debug information about time error and its bound.
248 249
    void printTimeErrorLowInfo() 
    {
250
      for (int i = 0; i < scalContents.getSize(); i++)
251 252 253
	std::cout << "    Time error estimate = " << scalContents[i]->est_t_sum 
		  << "    Time error bound = " << scalContents[i]->timeErrLow << "\n";
    }
254

255
    /// Returns \ref spaceIteration.
256 257
    inline int getSpaceIteration() 
    { 
258
      return spaceIteration; 
259
    }
260

261
    /// Sets \ref spaceIteration.
262 263
    inline void setSpaceIteration(int it) 
    { 
264
      spaceIteration = it; 
265
    }
266
  
267
    /// Returns \ref maxSpaceIteration.
268 269
    inline int getMaxSpaceIteration() 
    { 
270
      return maxSpaceIteration;
271
    }
272

273
    /// Sets \ref maxSpaceIteration.
274 275
    inline void setMaxSpaceIteration(int it) 
    { 
276
      maxSpaceIteration = it; 
277
    }
278
  
279
    /// Increments \ref spaceIteration by 1;
280 281
    inline void incSpaceIteration() 
    { 
282
      spaceIteration++; 
283
    }
284

285
    /// Sets \ref timestepIteration.
286 287
    inline void setTimestepIteration(int it) 
    { 
288
      timestepIteration = it; 
289
    }
290
  
291
    /// Returns \ref timestepIteration.
292 293
    inline int getTimestepIteration() 
    { 
294
      return timestepIteration; 
295
    }
296

297
    /// Increments \ref timestepIteration by 1;
298 299
    inline void incTimestepIteration() 
    { 
300
      timestepIteration++; 
301
    }
302

303
    /// Returns \ref maxTimestepIteration.
304 305
    inline int getMaxTimestepIteration() 
    { 
306
      return maxTimestepIteration; 
307
    }
308

309
    /// Sets \ref maxTimestepIteration.
310 311
    inline void setMaxTimestepIteration(int it) 
    { 
312
      maxTimestepIteration = it; 
313
    }
314
  
315
    /// Sets \ref timeIteration.
316 317
    inline void setTimeIteration(int it) 
    { 
318
      timeIteration = it; 
319
    }
320
  
321
    /// Returns \ref timeIteration.
322 323
    inline int getTimeIteration() 
    { 
324
      return timeIteration; 
325
    }
326

327
    /// Increments \ref timesIteration by 1;
328 329
    inline void incTimeIteration() 
    { 
330
      timeIteration++; 
331
    }
332

333
    /// Returns \ref maxTimeIteration.
334 335
    inline int getMaxTimeIteration() 
    { 
336
      return maxTimeIteration; 
337
    }
338

339
    /// Sets \ref maxTimeIteration.
340 341
    inline void setMaxTimeIteration(int it) 
    { 
342
      maxTimeIteration = it; 
343
    }
344
  
345
    /// Returns \ref timestepNumber.
346 347
    inline int getTimestepNumber() 
    { 
348
      return timestepNumber; 
349
    }
350

Thomas Witkowski's avatar
* Bla  
Thomas Witkowski committed
351
    /// Returns \ref nTimesteps.
352 353
    inline int getNumberOfTimesteps() 
    {
Thomas Witkowski's avatar
* Bla  
Thomas Witkowski committed
354 355 356
      return nTimesteps;
    }

357
    /// Increments \ref timestepNumber by 1;
358 359
    inline void incTimestepNumber() 
    { 
360
      timestepNumber++; 
361
    }
362

363
    /// Sets \ref est_sum.
364 365
    inline void setEstSum(double e, int index) 
    {
366
      scalContents[index]->est_sum = e;
367
    }
368

369
    /// Sets \ref est_max.
370 371
    inline void setEstMax(double e, int index) 
    {
372
      scalContents[index]->est_max = e;
373
    }
374

375
    /// Sets \ref est_max.
376 377
    inline void setTimeEstMax(double e, int index) 
    {
378
      scalContents[index]->est_t_max = e;
379
    }
380

381
    /// Sets \ref est_t_sum.
382 383
    inline void setTimeEstSum(double e, int index) 
    {
384
      scalContents[index]->est_t_sum = e;
385
    }
386

387
    /// Returns \ref est_sum.
388 389
    inline double getEstSum(int index) 
    { 
390
      return scalContents[index]->est_sum; 
391
    }
392

393
    /// Returns \ref est_t_sum.
394 395
    inline double getEstTSum(int index) 
    { 
396
      return scalContents[index]->est_t_sum; 
397
    }
398

399
    /// Returns \ref est_max.
400 401
    inline double getEstMax(int index) 
    { 
402
      return scalContents[index]->est_max; 
403
    }
404

405
    /// Returns \ref est_max.
406 407
    inline double getTimeEstMax(int index) 
    { 
408
      return scalContents[index]->est_t_max; 
409
    }
410

411
    /// Returns \ref est_t_sum.
412 413
    inline double getTimeEstSum(int index) 
    { 
414
      return scalContents[index]->est_t_sum; 
415
    }
416

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

423
    /// Sets \ref spaceTolerance.
424 425
    inline void setSpaceTolerance(int index, double tol) 
    { 
426
      scalContents[index]->spaceTolerance = tol; 
427
    }  
428

429
    /// Returns \ref timeTolerance.
430 431
    inline double getTimeTolerance(int index) 
    { 
432
      return scalContents[index]->timeTolerance; 
433
    }  
434

435
    /// Sets \ref time
436 437
    inline double setTime(double t) 
    { 
438
      time = t; 
439 440 441 442 443
      if (time > endTime) 
	time = endTime;
      if (time < startTime) 
	time = startTime;

444
      return time;
445
    }
446

447
    /// Gets \ref time
448 449
    inline double getTime() 
    { 
450
      return time; 
451
    }  
452

453
    /// Gets \ref &time
454 455
    inline double *getTimePtr() 
    { 
456
      return &time; 
457
    }  
458

459
    /// Sets \ref timestep
460 461
    inline double setTimestep(double t) 
    { 
462
      timestep = t; 
463
      if (timestep > maxTimestep)
464
	timestep = maxTimestep;
465
      if (timestep < minTimestep)
466
	timestep = minTimestep;
467
      if (time + timestep > endTime)
468 469
	timestep = endTime - time;
      
470
      return timestep;
471
    }
472

Thomas Witkowski's avatar
* Bla  
Thomas Witkowski committed
473 474 475 476
    /** \brief
     * Returns true, if the end time is reached and no more timestep
     * computations must be done.
     */
477 478
    inline bool reachedEndTime() 
    {
479 480 481 482
      if (nTimesteps > 0) 
	return !(timestepNumber < nTimesteps);

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

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

491
    /// Sets \ref minTimestep
492 493
    inline void setMinTimestep(double t) 
    { 
494
      minTimestep = t; 
495
    }
496

497
    /// Gets \ref minTimestep
498 499
    inline double getMinTimestep() 
    { 
500
      return minTimestep; 
501
    }  
502

503
    /// Sets \ref maxTimestep
504 505
    inline void setMaxTimestep(double t) 
    { 
506
      maxTimestep = t; 
507
    }
508

509
    /// Gets \ref maxTimestep
510 511
    inline double getMaxTimestep() 
    { 
512
      return maxTimestep; 
513
    }  
514

515
    /// Gets \ref &timestep
516 517
    inline double *getTimestepPtr() 
    { 
518
      return &timestep; 
519
    }  
520

521
    /// Sets \ref startTime = time
522 523
    inline void setStartTime(double time) 
    { 
524
      startTime = time; 
525
    }
526

527
    /// Sets \ref endTime = time
528 529
    inline void setEndTime(double time) 
    { 
530
      endTime = time; 
531
    }
532

533
    /// Returns \ref startTime
534 535
    inline double getStartTime() 
    { 
536
      return startTime; 
537
    }
538

539
    /// Returns \ref endTime
540 541
    inline double getEndTime() 
    { 
542
      return endTime; 
543
    }
544

545
    /// Returns \ref timeErrLow.
546 547
    inline double getTimeErrLow(int index) 
    { 
548
      return scalContents[index]->timeErrLow; 
549
    }  
550

551
    /// Returns whether coarsening is allowed or not.
552 553
    inline bool isCoarseningAllowed(int index) 
    {
554
      return (scalContents[index]->coarsenAllowed == 1);
555
    }
556

557
    /// Returns whether coarsening is allowed or not.
558 559
    inline bool isRefinementAllowed(int index) 
    {
560
      return (scalContents[index]->refinementAllowed == 1);
561
    }
562

563
    ///
564 565
    inline void allowRefinement(bool allow, int index) 
    {
566
      scalContents[index]->refinementAllowed = allow;
567
    }
568

569
    ///
570 571
    inline void allowCoarsening(bool allow, int index) 
    {
572
      scalContents[index]->coarsenAllowed = allow;
573
    }
574

575
    /// Returns \ref refineBisections
576 577
    inline const int getRefineBisections(int index) const 
    {
578
      return scalContents[index]->refineBisections;
579
    }
580

581
    /// Returns \ref coarseBisections
582 583
    inline const int getCoarseBisections(int index) const 
    {
584
      return scalContents[index]->coarseBisections;
585
    }
586

587 588
    inline int getSize() 
    { 
589
      return scalContents.getSize(); 
590
    }
591

592 593
    inline void setSolverIterations(int it) 
    {
594
      solverIterations = it;
595
    }
596
  
597 598
    inline int getSolverIterations() 
    {
599
      return solverIterations;
600
    }
601
  
602 603
    inline void setMaxSolverIterations(int it) 
    {
604
      maxSolverIterations = it;
605
    }
606
  
607 608
    inline int getMaxSolverIterations() 
    {
609
      return maxSolverIterations;
610
    }
611
  
612 613
    inline void setSolverTolerance(double tol) 
    {
614
      solverTolerance = tol;
615
    }
616
  
617 618
    inline double getSolverTolerance() 
    {
619
      return solverTolerance;
620
    }
621
  
622 623
    inline void setSolverResidual(double res) 
    {
624
      solverResidual = res;
625
    }
626
  
627 628
    inline double getSolverResidual() 
    {
629
      return solverResidual;
630
    }
631

632
    /// Returns true, if the adaptive procedure was deserialized from a file.
633 634
    const bool isDeserialized() const 
    {
635 636 637
      return isDeserialized_;
    }

638 639
    inline void setIsDeserialized(bool b) 
    {
640 641 642
      isDeserialized_ = b;
    }

643
    /// Creates new scalContents with the given size.
644 645
    void setScalContents(int newSize);

646 647 648 649 650 651 652 653 654 655 656 657 658 659 660
    /** \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
661
    void serialize(std::ostream& out);
662

Thomas Witkowski's avatar
Thomas Witkowski committed
663
    void deserialize(std::istream& in);
664 665

  protected:
666
    /// Name.
Thomas Witkowski's avatar
Thomas Witkowski committed
667
    std::string name;
668

669
    /// Current space iteration
670 671 672 673 674 675 676 677
    int spaceIteration;

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

678
    /// Current timestep iteration
679 680
    int timestepIteration;

681
    /// Maximal number of iterations for choosing a timestep
682 683
    int maxTimestepIteration;

684
    /// Current time iteration
685 686
    int timeIteration;

687
    /// Maximal number of time iterations
688 689
    int maxTimeIteration;

690
    /// Actual time, end of time interval for current time step
691 692
    double time;

693
    /// Initial time
694 695
    double startTime;

696
    /// Final time
697 698
    double endTime;

699
    /// Current time step size
700 701
    double timestep;