RosenbrockAdaptInstationary.cc 2.78 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
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
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
#include "time/RosenbrockAdaptInstationary.h"
#include "AdaptInfo.h"
#include "ProblemTimeInterface.h"

namespace AMDiS {

  void RosenbrockAdaptInstationary::oneTimestep()
  {
    FUNCNAME("RosenbrockAdaptInstationary::oneTimestep()");

    // estimate before first adaption
    if (adaptInfo->getTime() <= adaptInfo->getStartTime())
      problemIteration->oneIteration(adaptInfo, ESTIMATE);

    bool rejected = false;

    do {
      // increment time
      adaptInfo->setTime(adaptInfo->getTime() + adaptInfo->getTimestep());
      problemTime->setTime(adaptInfo);
      
      INFO(info, 6)("time = %e   timestep = %e\n",
		    adaptInfo->getTime(), adaptInfo->getTimestep());

      adaptInfo->setSpaceIteration(0);
      
      // do the iteration
      problemIteration->beginIteration(adaptInfo);
      problemIteration->oneIteration(adaptInfo, FULL_ITERATION);
      problemIteration->endIteration(adaptInfo);

      double errorEst = adaptInfo->getTimeEstSum(0);
      double newTimestep = 0.0;

      if (errorEst < timeTol) {
	double fac = 1.0;

	if (firstTimestep || succRejection) {
	  fac = pow((timeTol / errorEst), 1.0 / 3.0);
	} else {
	  fac = adaptInfo->getTimestep() / tauAcc * 
	    pow((timeTol * estAcc / (errorEst * errorEst)), 1.0 / 3.0);
	}
	fac = std::min(fac, 3.0);
	newTimestep = fac * adaptInfo->getTimestep();

	tauAcc = adaptInfo->getTimestep();
	estAcc = errorEst;

	lastTimestepRejected = false;
	succRejection = false;
      } else {
	double p = 3.0;

	if (lastTimestepRejected) {
	  succRejection = true;  

	  double reducedP = log(errorEst / estRej) / log(adaptInfo->getTimestep() / tauRej);
	  
	  if (reducedP < p && reducedP > 0.0)
	    p = reducedP;
	} 

	newTimestep = pow((timeTol / errorEst), 1.0 / p) * adaptInfo->getTimestep();
	
	tauRej = adaptInfo->getTimestep();
	estRej = errorEst;

	lastTimestepRejected = true;
      }
	
      newTimestep *= 0.95;

      if (errorEst < timeTol || firstTimestep) {
	INFO(info, 6)("Accepted timestep at time = %f   with timestep = %f\n",
		      adaptInfo->getTime()  - adaptInfo->getTimestep(),
		      adaptInfo->getTimestep());

	adaptInfo->setTimestep(newTimestep);

	if (firstTimestep)
	  firstTimestep = false;

	rejected = false;

	for (int i = 0; i < 3; i++)
	  INFO(info, 6)("time estimator for component %d = %f\n", 
			i, adaptInfo->getTimeEstSum(i));
      } else {
	INFO(info, 6)("Rejected timestep at time = %f   with timestep = %f\n",
		      adaptInfo->getTime()  - adaptInfo->getTimestep(),
		      adaptInfo->getTimestep());
	
	adaptInfo->setTime(adaptInfo->getTime() - adaptInfo->getTimestep());
	adaptInfo->setTimestep(newTimestep);

	rejected = true;
      }

    } while (rejected);

    rosenbrockStat->acceptTimestep();

    adaptInfo->setLastProcessedTimestep(adaptInfo->getTimestep());
    adaptInfo->incTimestepNumber();
  }

}