RosenbrockAdaptInstationary.cc 7.1 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
//
// 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.


13
14
15
16
17
18
#include "time/RosenbrockAdaptInstationary.h"
#include "AdaptInfo.h"
#include "ProblemTimeInterface.h"

namespace AMDiS {

19
20
21
22
23
24
25
26
27
28
29
  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),
30
      fixFirstTimesteps(0),
31
      tau(1.0),
Thomas Witkowski's avatar
Thomas Witkowski committed
32
33
34
35
36
      tauGamma(1.0),
      minusTauGamma(-1.0),
      invTauGamma(1.0),
      minusInvTauGamma(-1.0),
      dbgTimestepStudy(false)
37
  {
38
    FUNCNAME("RosenbrockAdaptInstationary::RosenbrockAdaptInstationary()");
39
40
    initConstructor(problemStat);
  }
41

42

43
44
45
46
47
48
49
50
51
52
53
54
55
  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),
      fixFirstTimesteps(0),
      tau(1.0),
Thomas Witkowski's avatar
Thomas Witkowski committed
56
57
58
      tauGamma(1.0),
      minusTauGamma(-1.0),
      invTauGamma(1.0),
59
60
      minusInvTauGamma(-1.0),
      dbgTimestepStudy(false)
61
62
63
64
  {
    FUNCNAME("RosenbrockAdaptInstationary::RosenbrockAdaptInstationary()");
    initConstructor(&problemStat);
  }
65

66

67
68
  void RosenbrockAdaptInstationary::initConstructor(RosenbrockStationary *problemStat)
  {
69
70
71
72
73
    FUNCNAME("RosenbrockAdaptInstationary::initConstructor()");

    std::string str("");
    std::string initFileStr(name + "->rosenbrock method");
    Parameters::get(initFileStr, str);
74
    RosenbrockMethodCreator *creator = 
75
      dynamic_cast<RosenbrockMethodCreator*>(CreatorMap<RosenbrockMethod>::getCreator(str, initFileStr));
76
    rosenbrockMethod = creator->create();
77
    
78
    TEST_EXIT_DBG(rosenbrockMethod)("This should not happen!\n");
79
    
80
    Parameters::get(name + "->fix first timesteps", fixFirstTimesteps);
81
    problemStat->setRosenbrockMethod(rosenbrockMethod);
82
    
83
    adaptInfo->setRosenbrockMode(true);
84
    
Thomas Witkowski's avatar
Thomas Witkowski committed
85
86
    problemStat->setTauGamma(&tauGamma, &minusTauGamma, 
			     &invTauGamma, &minusInvTauGamma);
87
    problemStat->setTau(&tau);
Thomas Witkowski's avatar
Thomas Witkowski committed
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
    
    Parameters::get(name + "->rosenbrock->timestep study", dbgTimestepStudy);
    Parameters::get(name + "->rosenbrock->timestep study steps", dbgTimesteps);
  }


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

    firstTimestep = true;
    lastTimestepRejected = false;
    succRejection = false;
    fixFirstTimesteps = 0;
    tau = 1.0;
    tauGamma = 1.0;    
    minusTauGamma = -1.0;
    invTauGamma = 1.0;
    minusInvTauGamma = -1.0;
107
  }
108

109

110
111
  void RosenbrockAdaptInstationary::oneTimestep()
  {
112
    FUNCNAME("RosenbrockAdaptInstationary::oneTimestep()");
113
    
114
115
116
    // estimate before first adaption
    if (adaptInfo->getTime() <= adaptInfo->getStartTime())
      problemIteration->oneIteration(adaptInfo, ESTIMATE);
117
    
118
    bool rejected = false;
119
    double timeTol = adaptInfo->getTimeTolerance(0);
Thomas Witkowski's avatar
Thomas Witkowski committed
120
121
122

    int studyTimestep = -1;

123
    do {
Thomas Witkowski's avatar
Thomas Witkowski committed
124
125
126
127
128
129
130
131
132
133
      if (dbgTimestepStudy) {
	if (studyTimestep >= 0)
	  adaptInfo->setTime(adaptInfo->getTime() - dbgTimesteps[studyTimestep]);

	studyTimestep++;
	TEST_EXIT(studyTimestep < dbgTimesteps.size())("Should not happen!\n");
	
	adaptInfo->setTimestep(dbgTimesteps[studyTimestep]);
      }

134
135
136
      // increment time
      adaptInfo->setTime(adaptInfo->getTime() + adaptInfo->getTimestep());
      problemTime->setTime(adaptInfo);
137
      tau = adaptInfo->getTimestep();
Thomas Witkowski's avatar
Thomas Witkowski committed
138
139
140
141
142

      tauGamma = tau * rosenbrockMethod->getGamma();
      minusTauGamma = -tauGamma;
      invTauGamma = 1.0 / (tau * rosenbrockMethod->getGamma());
      minusInvTauGamma = -1.0;
143
144
145
      
      INFO(info, 6)("time = %e   timestep = %e\n",
		    adaptInfo->getTime(), adaptInfo->getTimestep());
146
      
147
148
149
150
151
152
153
154
155
      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;
156
      double order = rosenbrockMethod->getOrder();
157

158
      if (errorEst < timeTol || fixFirstTimesteps > 0 || firstTimestep) {
159
	double fac = 1.0;
160
161
162
	
	if (fixFirstTimesteps > 0) {
	  newTimestep = adaptInfo->getTimestep();
163
	} else {
164
165
166
167
168
169
170
171
172
	  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;
173
174
175
176
	}

	tauAcc = adaptInfo->getTimestep();
	estAcc = errorEst;
177
	
178
179
180
181
182
	lastTimestepRejected = false;
	succRejection = false;
      } else {
	if (lastTimestepRejected) {
	  succRejection = true;  
183
	  
184
185
	  double reducedP = 
	    log(errorEst / estRej) / log(adaptInfo->getTimestep() / tauRej);
186
	  
187
188
	  if (reducedP < order && reducedP > 0.0)
	    order = reducedP;
189
190
	} 
	
191
192
193
	newTimestep = pow((timeTol / errorEst), 1.0 / order) * adaptInfo->getTimestep();
	newTimestep *= 0.95;

194
195
	tauRej = adaptInfo->getTimestep();
	estRej = errorEst;
196
	
197
198
	lastTimestepRejected = true;
      }
199
200
      
      
201
202
      if (errorEst < timeTol || fixFirstTimesteps 
	  || !(adaptInfo->getTimestep()>adaptInfo->getMinTimestep()) ) {
203
	INFO(info, 6)("Accepted timestep at time = %e   with timestep = %e\n",
204
205
206
207
208
209
210
		      adaptInfo->getTime()  - adaptInfo->getTimestep(),
		      adaptInfo->getTimestep());

	adaptInfo->setTimestep(newTimestep);

	rejected = false;

211
	for (int i = 0; i < adaptInfo->getSize(); i++)
212
	  INFO(info, 6)("time estimate for component %d = %e\n", 
213
			i, adaptInfo->getTimeEstSum(i));
214

215
        if (errorEst > timeTol)
216
	  MSG("Accepted timestep but tolerance not reached \n");
217

218
      } else {
219
220
221
222
	for (int i = 0; i < adaptInfo->getSize(); i++)
	  INFO(info, 6)("time estimate for component %d = %e\n", 
			i, adaptInfo->getTimeEstSum(i));

223
	INFO(info, 6)("Rejected timestep at time = %e   with timestep = %e\n",
224
225
226
227
228
229
230
231
232
		      adaptInfo->getTime()  - adaptInfo->getTimestep(),
		      adaptInfo->getTimestep());
	
	adaptInfo->setTime(adaptInfo->getTime() - adaptInfo->getTimestep());
	adaptInfo->setTimestep(newTimestep);

	rejected = true;
      }

233
234
235
      if (firstTimestep)
	firstTimestep = false;

Thomas Witkowski's avatar
Thomas Witkowski committed
236
237
238
      
    } while (rejected || 
	     (dbgTimestepStudy && (studyTimestep + 1 < dbgTimesteps.size())));
239
240
241
242
243
244
245
246

    rosenbrockStat->acceptTimestep();

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

}