RosenbrockAdaptInstationary.cc 4.72 KB
Newer Older
1
2
3
4
5
6
#include "time/RosenbrockAdaptInstationary.h"
#include "AdaptInfo.h"
#include "ProblemTimeInterface.h"

namespace AMDiS {

7
8
9
10
11
12
13
14
15
16
17
  RosenbrockAdaptInstationary::RosenbrockAdaptInstationary(std::string name, 
							   RosenbrockStationary *problemStat,
							   AdaptInfo *info,
							   ProblemTimeInterface *problemInstat,
							   AdaptInfo *initialInfo,
							   time_t initialTimestamp)
    : AdaptInstationary(name, problemStat, info, problemInstat, initialInfo, initialTimestamp),
      rosenbrockStat(problemStat),
      firstTimestep(true),
      lastTimestepRejected(false),
      succRejection(false),
18
      fixFirstTimesteps(0),
19
20
      tau(1.0),
      tauGamma(1.0)
21
  {
22
23
    FUNCNAME("RosenbrockAdaptInstationary::RosenbrockAdaptInstationary()");

24
25
26
27
28
29
    std::string str;
    GET_PARAMETER(0, name + "->rosenbrock method", &str);
    RosenbrockMethodCreator *creator = 
      dynamic_cast<RosenbrockMethodCreator*>(CreatorMap<RosenbrockMethod>::getCreator(str));
    rosenbrockMethod = creator->create();

30
31
32
    TEST_EXIT_DBG(rosenbrockMethod)("This should not happen!\n");

    GET_PARAMETER(0, name + "->fix first timesteps", "%d", &fixFirstTimesteps);
33
    problemStat->setRosenbrockMethod(rosenbrockMethod);
34
35

    adaptInfo->setRosenbrockMode(true);
36
37

    problemStat->setTauGamma(&tauGamma);
38
    problemStat->setTau(&tau);
39
40
41
  }


42
43
44
45
46
47
48
49
50
  void RosenbrockAdaptInstationary::oneTimestep()
  {
    FUNCNAME("RosenbrockAdaptInstationary::oneTimestep()");

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

    bool rejected = false;
51
52
    double timeTol = adaptInfo->getTimeTolerance(0);
    
53
54
55
56
    do {
      // increment time
      adaptInfo->setTime(adaptInfo->getTime() + adaptInfo->getTimestep());
      problemTime->setTime(adaptInfo);
57
58
      tau = adaptInfo->getTimestep();
      tauGamma = 1.0 / (tau * rosenbrockMethod->getGamma());
59
60
61
62
63
64
65
66
67
68
69
70
71
      
      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;
72
      double order = rosenbrockMethod->getOrder();
73

74
      if (errorEst < timeTol || fixFirstTimesteps > 0 || firstTimestep) {
75
	double fac = 1.0;
76
77
78
	
	if (fixFirstTimesteps > 0) {
	  newTimestep = adaptInfo->getTimestep();
79
	} else {
80
81
82
83
84
85
86
87
88
	  if (firstTimestep || succRejection) {
	    fac = pow((timeTol / errorEst), 1.0 / order);
	  } else {
	    fac = adaptInfo->getTimestep() / tauAcc * 
	      pow((timeTol * estAcc / (errorEst * errorEst)), 1.0 / order);
	  }
	  fac = std::min(fac, 3.0);
	  newTimestep = fac * adaptInfo->getTimestep();
	  newTimestep *= 0.95;
89
90
91
92
	}

	tauAcc = adaptInfo->getTimestep();
	estAcc = errorEst;
93
	
94
95
96
97
98
	lastTimestepRejected = false;
	succRejection = false;
      } else {
	if (lastTimestepRejected) {
	  succRejection = true;  
99
	  
100
101
	  double reducedP = 
	    log(errorEst / estRej) / log(adaptInfo->getTimestep() / tauRej);
102
	  
103
104
	  if (reducedP < order && reducedP > 0.0)
	    order = reducedP;
105
106
	} 
	
107
108
109
	newTimestep = pow((timeTol / errorEst), 1.0 / order) * adaptInfo->getTimestep();
	newTimestep *= 0.95;

110
111
	tauRej = adaptInfo->getTimestep();
	estRej = errorEst;
112
	
113
114
	lastTimestepRejected = true;
      }
115
              
116

117
118
      if (errorEst < timeTol || fixFirstTimesteps 
	  || !(adaptInfo->getTimestep()>adaptInfo->getMinTimestep()) ) {
119
	INFO(info, 6)("Accepted timestep at time = %e   with timestep = %e\n",
120
121
122
123
124
125
126
		      adaptInfo->getTime()  - adaptInfo->getTimestep(),
		      adaptInfo->getTimestep());

	adaptInfo->setTimestep(newTimestep);

	rejected = false;

127
	for (int i = 0; i < adaptInfo->getSize(); i++)
128
	  INFO(info, 6)("time estimate for component %d = %e\n", 
129
			i, adaptInfo->getTimeEstSum(i));
130
131
132
133

        if(errorEst > timeTol)
	  WARNING("Accepted timestep but tolerance not reached \n" );

134
      } else {
135
136
137
138
	for (int i = 0; i < adaptInfo->getSize(); i++)
	  INFO(info, 6)("time estimate for component %d = %e\n", 
			i, adaptInfo->getTimeEstSum(i));

139
	INFO(info, 6)("Rejected timestep at time = %e   with timestep = %e\n",
140
141
142
143
144
145
146
147
148
		      adaptInfo->getTime()  - adaptInfo->getTimestep(),
		      adaptInfo->getTimestep());
	
	adaptInfo->setTime(adaptInfo->getTime() - adaptInfo->getTimestep());
	adaptInfo->setTimestep(newTimestep);

	rejected = true;
      }

149
150
151
152
153
154
      if (firstTimestep)
	firstTimestep = false;

      if (fixFirstTimesteps > 0)
	fixFirstTimesteps--;

155
156
157
158
159
160
161
162
163
    } while (rejected);

    rosenbrockStat->acceptTimestep();

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

}