AdaptInfo.h 19.8 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 47 48 49 50 51 52 53 54
// ============================================================================
// ==                                                                        ==
// == 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 "MemoryManager.h"
#include "MatrixVector.h"
#include "Parameters.h"
#include "Serializable.h"

namespace AMDiS {

  // ===========================================================================
  // ===== class AdaptInfo =====================================================
  // ===========================================================================

  /**
   * \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:
      /** \brief
       * Constructor.
       */
Thomas Witkowski's avatar
Thomas Witkowski committed
55
      ScalContent(std::string prefix) 
56 57 58 59 60 61 62 63 64 65
	: 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),
66
  	  coarseBisections(1)
67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86
      {
	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;
87
      }
88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152

      /** \brief
       * Sum of all error estimates
       */
      double est_sum;

      /** \brief
       * Sum of all time error estimates
       */
      double est_t_sum;

      /** \brief
       * maximal local error estimate.
       */
      double est_max;

      /** \brief
       * Maximum of all time error estimates
       */
      double est_t_max;

      /** \brief
       * Tolerance for the (absolute or relative) error
       */
      double spaceTolerance;

      /** \brief
       * Time tolerance.
       */
      double timeTolerance;

      /** \brief
       * Lower bound for the time error.
       */
      double timeErrLow;

      /** \brief
       * true if coarsening is allowed, false otherwise.
       */
      int coarsenAllowed;

      /** \brief
       * true if refinement is allowed, false otherwise.
       */
      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
       */
      int  refineBisections;

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

  public:
    /** \brief
     * Constructor.
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
153
    AdaptInfo(std::string name_, int size = 1) 
154 155 156 157 158 159 160 161 162 163 164 165 166 167
      : 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
168
	nTimesteps(0),
169 170 171 172
	solverIterations(0),
	maxSolverIterations(0),
	solverTolerance(1e-8),
	solverResidual(0.0),
173 174
        scalContents(size),
	isDeserialized_(false)
175 176 177 178 179 180 181 182 183 184 185 186
    {
      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
187 188
      GET_PARAMETER(0, name_ + "->number of timesteps", "%d", &nTimesteps);

189 190 191 192 193 194
      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
195
	  scalContents[i] = new ScalContent(name + std::string(number)); 
196 197
	}
      }
198
    }
199 200 201 202 203

    /** \brief
     * Destructor.
     */
    virtual ~AdaptInfo() {
204
      for (int i = 0;  i < scalContents.getSize(); i++) {
205 206
	delete scalContents[i];
      }
207
    }
208 209 210 211 212 213 214 215 216 217 218

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

      GET_PARAMETER(0, name + "->timestep", "%f", &timestep);
221
    }
222 223 224 225 226

    /** \brief
     * Returns whether space tolerance is reached.
     */
    virtual bool spaceToleranceReached() {
227 228 229
      int size = scalContents.getSize();
      for (int i = 0; i < size; i++) {
	if (!(scalContents[i]->est_sum < scalContents[i]->spaceTolerance)) {
230 231 232 233
	  return false;
	}
      }
      return true;
234
    }
235 236 237 238 239

    /** \brief
     * Returns whether space tolerance of component i is reached.
     */
    virtual bool spaceToleranceReached(int i) {
240
      if (!(scalContents[i]->est_sum < scalContents[i]->spaceTolerance)) {
241 242 243 244
	return false;
      } else {
	return true;
      }
245
    }
246 247 248 249 250

    /** \brief
     * Returns whether time tolerance is reached.
     */
    virtual bool timeToleranceReached() {
251 252 253
      int size = scalContents.getSize();
      for (int i = 0; i < size; i++) {
	if (!(scalContents[i]->est_t_sum < scalContents[i]->timeTolerance)) {
254 255 256 257
	  return false;
	}
      }
      return true;
258
    }
259 260 261 262 263

    /** \brief
     * Returns whether time tolerance of component i is reached.
     */
    virtual bool timeToleranceReached(int i) {
264
      if (!(scalContents[i]->est_t_sum < scalContents[i]->timeTolerance)) {
265 266 267 268
	return false;
      } else {
	return true;
      }
269
    }
270 271 272 273 274 275 276 277 278 279 280 281

    /** \brief
     * Returns whether time error is under its lower bound.
     */
    virtual bool timeErrorLow() {
      int size = scalContents.getSize();
      for (int i = 0; i < size; i++) {
	if (!(scalContents[i]->est_t_sum < scalContents[i]->timeErrLow)) {
	  return false;
	}
      }
      return true;
282 283 284 285 286 287 288 289 290 291 292
    }

    /** \brief
     * Print debug information about time error and its bound.
     */
    void printTimeErrorLowInfo() {
      for (int i = 0; i < scalContents.getSize(); i++) {
	std::cout << "    Time error estimate = " << scalContents[i]->est_t_sum 
		  << "    Time error bound = " << scalContents[i]->timeErrLow << "\n";
      }
    }
293 294 295 296 297 298

    /** \brief
     * Returns \ref spaceIteration.
     */
    inline int getSpaceIteration() { 
      return spaceIteration; 
299
    }
300 301 302 303 304 305

    /** \brief
     * Sets \ref spaceIteration.
     */
    inline void setSpaceIteration(int it) { 
      spaceIteration = it; 
306
    }
307 308 309 310 311 312
  
    /** \brief
     * Returns \ref maxSpaceIteration.
     */
    inline int getMaxSpaceIteration() { 
      return maxSpaceIteration;
313
    }
314 315 316 317 318 319

    /** \brief
     * Sets \ref maxSpaceIteration.
     */
    inline void setMaxSpaceIteration(int it) { 
      maxSpaceIteration = it; 
320
    }
321 322 323 324 325 326
  
    /** \brief
     * Increments \ref spaceIteration by 1;
     */
    inline void incSpaceIteration() { 
      spaceIteration++; 
327
    }
328 329 330 331 332 333

    /** \brief
     * Sets \ref timestepIteration.
     */
    inline void setTimestepIteration(int it) { 
      timestepIteration = it; 
334
    }
335 336 337 338 339 340
  
    /** \brief
     * Returns \ref timestepIteration.
     */
    inline int getTimestepIteration() { 
      return timestepIteration; 
341
    }
342 343 344 345 346 347

    /** \brief
     * Increments \ref timestepIteration by 1;
     */
    inline void incTimestepIteration() { 
      timestepIteration++; 
348
    }
349 350 351 352 353 354

    /** \brief
     * Returns \ref maxTimestepIteration.
     */
    inline int getMaxTimestepIteration() { 
      return maxTimestepIteration; 
355
    }
356 357 358 359 360 361

    /** \brief
     * Sets \ref maxTimestepIteration.
     */
    inline void setMaxTimestepIteration(int it) { 
      maxTimestepIteration = it; 
362
    }
363 364 365 366 367 368
  
    /** \brief
     * Sets \ref timeIteration.
     */
    inline void setTimeIteration(int it) { 
      timeIteration = it; 
369
    }
370 371 372 373 374 375
  
    /** \brief
     * Returns \ref timeIteration.
     */
    inline int getTimeIteration() { 
      return timeIteration; 
376
    }
377 378 379 380 381 382

    /** \brief
     * Increments \ref timesIteration by 1;
     */
    inline void incTimeIteration() { 
      timeIteration++; 
383
    }
384 385 386 387 388 389

    /** \brief
     * Returns \ref maxTimeIteration.
     */
    inline int getMaxTimeIteration() { 
      return maxTimeIteration; 
390
    }
391 392 393 394 395 396

    /** \brief
     * Sets \ref maxTimeIteration.
     */
    inline void setMaxTimeIteration(int it) { 
      maxTimeIteration = it; 
397
    }
398 399 400 401 402 403
  
    /** \brief
     * Returns \ref timestepNumber.
     */
    inline int getTimestepNumber() { 
      return timestepNumber; 
404
    }
405

Thomas Witkowski's avatar
* Bla  
Thomas Witkowski committed
406 407 408 409 410
    /// Returns \ref nTimesteps.
    inline int getNumberOfTimesteps() {
      return nTimesteps;
    }

411 412 413 414 415
    /** \brief
     * Increments \ref timestepNumber by 1;
     */
    inline void incTimestepNumber() { 
      timestepNumber++; 
416
    }
417 418 419 420 421 422

    /** \brief
     * Sets \ref est_sum.
     */
    inline void setEstSum(double e, int index) {
      scalContents[index]->est_sum = e;
423
    }
424 425 426 427 428 429

    /** \brief
     * Sets \ref est_max.
     */
    inline void setEstMax(double e, int index) {
      scalContents[index]->est_max = e;
430
    }
431 432 433 434 435 436

    /** \brief
     * Sets \ref est_max.
     */
    inline void setTimeEstMax(double e, int index) {
      scalContents[index]->est_t_max = e;
437
    }
438 439 440 441 442 443

    /** \brief
     * Sets \ref est_t_sum.
     */
    inline void setTimeEstSum(double e, int index) {
      scalContents[index]->est_t_sum = e;
444
    }
445 446 447 448 449 450

    /** \brief
     * Returns \ref est_sum.
     */
    inline double getEstSum(int index) { 
      return scalContents[index]->est_sum; 
451
    }
452

453 454 455 456 457
    /** \brief
     * Returns \ref est_t_sum.
     */
    inline double getEstTSum(int index) { 
      return scalContents[index]->est_t_sum; 
458
    }
459

460 461 462 463 464
    /** \brief
     * Returns \ref est_max.
     */
    inline double getEstMax(int index) { 
      return scalContents[index]->est_max; 
465
    }
466 467 468 469 470 471

    /** \brief
     * Returns \ref est_max.
     */
    inline double getTimeEstMax(int index) { 
      return scalContents[index]->est_t_max; 
472
    }
473 474 475 476 477 478

    /** \brief
     * Returns \ref est_t_sum.
     */
    inline double getTimeEstSum(int index) { 
      return scalContents[index]->est_t_sum; 
479
    }
480 481 482 483 484 485

    /** \brief
     * Returns \ref spaceTolerance.
     */
    inline double getSpaceTolerance(int index) { 
      return scalContents[index]->spaceTolerance; 
486
    }  
487 488 489 490 491 492

    /** \brief
     * Sets \ref spaceTolerance.
     */
    inline void setSpaceTolerance(int index, double tol) { 
      scalContents[index]->spaceTolerance = tol; 
493
    }  
494 495 496 497 498 499

    /** \brief
     * Returns \ref timeTolerance.
     */
    inline double getTimeTolerance(int index) { 
      return scalContents[index]->timeTolerance; 
500
    }  
501 502 503 504 505 506

    /** \brief
     * Sets \ref time
     */
    inline double setTime(double t) { 
      time = t; 
507 508 509 510 511
      if (time > endTime) 
	time = endTime;
      if (time < startTime) 
	time = startTime;

512
      return time;
513
    }
514 515 516 517 518 519

    /** \brief
     * Gets \ref time
     */
    inline double getTime() { 
      return time; 
520
    }  
521 522 523 524 525 526

    /** \brief
     * Gets \ref &time
     */
    inline double *getTimePtr() { 
      return &time; 
527
    }  
528 529 530 531 532 533

    /** \brief
     * Sets \ref timestep
     */
    inline double setTimestep(double t) { 
      timestep = t; 
534
      if (timestep > maxTimestep) {
535
	timestep = maxTimestep;
536 537
      }
      if (timestep < minTimestep) {
538
	timestep = minTimestep;
539 540
      }
      if (time + timestep > endTime) {
541
	timestep = endTime - time;
542
      }
543
      
544
      return timestep;
545
    }
546

Thomas Witkowski's avatar
* Bla  
Thomas Witkowski committed
547 548 549 550 551
    /** \brief
     * Returns true, if the end time is reached and no more timestep
     * computations must be done.
     */
    inline bool reachedEndTime() {
552 553 554 555
      if (nTimesteps > 0) 
	return !(timestepNumber < nTimesteps);

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

558 559 560 561 562
    /** \brief
     * Gets \ref timestep
     */
    inline double getTimestep() { 
      return timestep; 
563
    }
564 565 566 567 568 569

    /** \brief
     * Sets \ref minTimestep
     */
    inline void setMinTimestep(double t) { 
      minTimestep = t; 
570
    }
571 572 573 574 575 576

    /** \brief
     * Gets \ref minTimestep
     */
    inline double getMinTimestep() { 
      return minTimestep; 
577
    }  
578 579 580 581 582 583

    /** \brief
     * Sets \ref maxTimestep
     */
    inline void setMaxTimestep(double t) { 
      maxTimestep = t; 
584
    }
585 586 587 588 589 590

    /** \brief
     * Gets \ref maxTimestep
     */
    inline double getMaxTimestep() { 
      return maxTimestep; 
591
    }  
592 593 594 595 596 597

    /** \brief
     * Gets \ref &timestep
     */
    inline double *getTimestepPtr() { 
      return &timestep; 
598
    }  
599 600 601 602 603 604

    /** \brief
     * Sets \ref startTime = time
     */
    inline void setStartTime(double time) { 
      startTime = time; 
605
    }
606 607 608 609 610 611

    /** \brief
     * Sets \ref endTime = time
     */
    inline void setEndTime(double time) { 
      endTime = time; 
612
    }
613 614 615 616 617 618

    /** \brief
     * Returns \ref startTime
     */
    inline double getStartTime() { 
      return startTime; 
619
    }
620 621 622 623 624 625

    /** \brief
     * Returns \ref endTime
     */
    inline double getEndTime() { 
      return endTime; 
626
    }
627 628 629 630 631 632

    /** \brief
     * Returns \ref timeErrLow.
     */
    inline double getTimeErrLow(int index) { 
      return scalContents[index]->timeErrLow; 
633
    }  
634 635 636 637 638 639

    /** \brief
     * Returns whether coarsening is allowed or not.
     */
    inline bool isCoarseningAllowed(int index) {
      return (scalContents[index]->coarsenAllowed == 1);
640
    }
641 642 643 644 645 646

    /** \brief
     * Returns whether coarsening is allowed or not.
     */
    inline bool isRefinementAllowed(int index) {
      return (scalContents[index]->refinementAllowed == 1);
647
    }
648 649 650 651 652 653

    /** \brief
     * 
     */
    inline void allowRefinement(bool allow, int index) {
      scalContents[index]->refinementAllowed = allow;
654
    }
655 656 657 658 659 660

    /** \brief
     * 
     */
    inline void allowCoarsening(bool allow, int index) {
      scalContents[index]->coarsenAllowed = allow;
661
    }
662 663 664 665 666 667

    /** \brief
     * Returns \ref refineBisections
     */
    inline const int getRefineBisections(int index) const {
      return scalContents[index]->refineBisections;
668
    }
669 670 671 672 673 674

    /** \brief
     * Returns \ref coarseBisections
     */
    inline const int getCoarseBisections(int index) const {
      return scalContents[index]->coarseBisections;
675
    }
676 677 678

    inline int getSize() { 
      return scalContents.getSize(); 
679
    }
680 681 682

    inline void setSolverIterations(int it) {
      solverIterations = it;
683
    }
684 685 686
  
    inline int getSolverIterations() {
      return solverIterations;
687
    }
688 689 690
  
    inline void setMaxSolverIterations(int it) {
      maxSolverIterations = it;
691
    }
692 693 694
  
    inline int getMaxSolverIterations() {
      return maxSolverIterations;
695
    }
696 697 698
  
    inline void setSolverTolerance(double tol) {
      solverTolerance = tol;
699
    }
700 701 702
  
    inline double getSolverTolerance() {
      return solverTolerance;
703
    }
704 705 706
  
    inline void setSolverResidual(double res) {
      solverResidual = res;
707
    }
708 709 710
  
    inline double getSolverResidual() {
      return solverResidual;
711
    }
712

713 714 715 716 717 718 719 720 721 722 723
    /** \brief
     * Returns true, if the adaptive procedure was deserialized from a file.
     */
    const bool isDeserialized() const {
      return isDeserialized_;
    }

    inline void setIsDeserialized(bool b) {
      isDeserialized_ = b;
    }

724 725 726 727 728
    /** \brief
     * Creates new scalContents with the given size.
     */
    void setScalContents(int newSize);

729 730 731 732 733 734 735 736 737 738 739 740 741 742 743
    /** \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;
    }

744
    // ===== Serialiazable implementation =====
Thomas Witkowski's avatar
Thomas Witkowski committed
745
    void serialize(std::ostream& out);
746

Thomas Witkowski's avatar
Thomas Witkowski committed
747
    void deserialize(std::istream& in);
748 749 750 751 752

  protected:
    /** \brief
     * Name.
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
753
    std::string name;
754 755 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 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819

    /** \brief
     * Current space iteration
     */
    int spaceIteration;

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

    /** \brief
     * Current timestep iteration
     */
    int timestepIteration;

    /** \brief
     * Maximal number of iterations for choosing a timestep
     */
    int maxTimestepIteration;

    /** \brief
     * Current time iteration
     */
    int timeIteration;

    /** \brief
     * Maximal number of time iterations
     */
    int maxTimeIteration;

    /** \brief
     * Actual time, end of time interval for current time step
     */
    double time;

    /** \brief
     * Initial time
     */
    double startTime;

    /** \brief
     * Final time
     */
    double endTime;

    /** \brief
     * Current time step size
     */
    double timestep;

    /** \brief
     * Minimal step size
     */
    double minTimestep;

    /** \brief
     * Maximal step size
     */
    double maxTimestep;

    /** \brief
     * Number of current time step
     */
    int timestepNumber;
Thomas Witkowski's avatar
* Bla  
Thomas Witkowski committed
820 821 822 823 824 825 826

    /** \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;
827 828 829 830 831 832 833 834 835 836 837 838
  
    /** \brief
     * number of iterations needed of linear or nonlinear solver
     */
    int solverIterations;


    /** \brief
     * maximal number of iterations needed of linear or nonlinear solver
     */
    int maxSolverIterations;

839 840 841
    /** \brief
     *
     */
842 843
    double solverTolerance;

844 845 846
    /** \brief
     *
     */
847 848 849 850 851 852
    double solverResidual;

    /** \brief
     * Scalar adapt infos.
     */
    Vector<ScalContent*> scalContents;
853 854 855 856 857

    /** \brief
     * Is true, if the adaptive procedure was deserialized from a file.
     */
    bool isDeserialized_;
858 859 860 861 862
  };

}

#endif //  AMDIS_ADAPTINFO_H