RosenbrockAdaptInstationary.cc 5.77 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
32
      tau(1.0),
      tauGamma(1.0)
33
  {
34
    FUNCNAME("RosenbrockAdaptInstationary::RosenbrockAdaptInstationary()");
35
36
    initConstructor(problemStat);
  }
37

38
39
40
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),
      tauGamma(1.0)
  {
    FUNCNAME("RosenbrockAdaptInstationary::RosenbrockAdaptInstationary()");
    initConstructor(&problemStat);
  }
56
57
58

  void RosenbrockAdaptInstationary::initConstructor(RosenbrockStationary *problemStat)
  {
59
    std::string str;
60
    Parameters::get(name + "->rosenbrock method", str);
61
    RosenbrockMethodCreator *creator = 
62
      dynamic_cast<RosenbrockMethodCreator*>(CreatorMap<RosenbrockMethod>::getCreator(str));
63
    rosenbrockMethod = creator->create();
64
    
65
    TEST_EXIT_DBG(rosenbrockMethod)("This should not happen!\n");
66
    
67
    Parameters::get(name + "->fix first timesteps", fixFirstTimesteps);
68
    problemStat->setRosenbrockMethod(rosenbrockMethod);
69
    
70
    adaptInfo->setRosenbrockMode(true);
71
    
72
    problemStat->setTauGamma(&tauGamma);
73
    problemStat->setTau(&tau);
74
  }
75

76
77
  void RosenbrockAdaptInstationary::oneTimestep()
  {
78
    FUNCNAME("RosenbrockAdaptInstationary::oneTimestep()");
79
    
80
81
82
    // estimate before first adaption
    if (adaptInfo->getTime() <= adaptInfo->getStartTime())
      problemIteration->oneIteration(adaptInfo, ESTIMATE);
83
    
84
    bool rejected = false;
85
86
    double timeTol = adaptInfo->getTimeTolerance(0);
    
87
88
89
90
    do {
      // increment time
      adaptInfo->setTime(adaptInfo->getTime() + adaptInfo->getTimestep());
      problemTime->setTime(adaptInfo);
91
92
      tau = adaptInfo->getTimestep();
      tauGamma = 1.0 / (tau * rosenbrockMethod->getGamma());
93
94
95
      
      INFO(info, 6)("time = %e   timestep = %e\n",
		    adaptInfo->getTime(), adaptInfo->getTimestep());
96
      
97
98
99
100
101
102
103
104
105
      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;
106
      double order = rosenbrockMethod->getOrder();
107

108
      if (errorEst < timeTol || fixFirstTimesteps > 0 || firstTimestep) {
109
	double fac = 1.0;
110
111
112
	
	if (fixFirstTimesteps > 0) {
	  newTimestep = adaptInfo->getTimestep();
113
	} else {
114
115
116
117
118
119
120
121
122
	  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;
123
124
125
126
	}

	tauAcc = adaptInfo->getTimestep();
	estAcc = errorEst;
127
	
128
129
130
131
132
	lastTimestepRejected = false;
	succRejection = false;
      } else {
	if (lastTimestepRejected) {
	  succRejection = true;  
133
	  
134
135
	  double reducedP = 
	    log(errorEst / estRej) / log(adaptInfo->getTimestep() / tauRej);
136
	  
137
138
	  if (reducedP < order && reducedP > 0.0)
	    order = reducedP;
139
140
	} 
	
141
142
143
	newTimestep = pow((timeTol / errorEst), 1.0 / order) * adaptInfo->getTimestep();
	newTimestep *= 0.95;

144
145
	tauRej = adaptInfo->getTimestep();
	estRej = errorEst;
146
	
147
148
	lastTimestepRejected = true;
      }
149
150
      
      
151
152
      if (errorEst < timeTol || fixFirstTimesteps 
	  || !(adaptInfo->getTimestep()>adaptInfo->getMinTimestep()) ) {
153
	INFO(info, 6)("Accepted timestep at time = %e   with timestep = %e\n",
154
155
156
157
158
159
160
		      adaptInfo->getTime()  - adaptInfo->getTimestep(),
		      adaptInfo->getTimestep());

	adaptInfo->setTimestep(newTimestep);

	rejected = false;

161
	for (int i = 0; i < adaptInfo->getSize(); i++)
162
	  INFO(info, 6)("time estimate for component %d = %e\n", 
163
			i, adaptInfo->getTimeEstSum(i));
164
165
166
167

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

168
      } else {
169
170
171
172
	for (int i = 0; i < adaptInfo->getSize(); i++)
	  INFO(info, 6)("time estimate for component %d = %e\n", 
			i, adaptInfo->getTimeEstSum(i));

173
	INFO(info, 6)("Rejected timestep at time = %e   with timestep = %e\n",
174
175
176
177
178
179
180
181
182
		      adaptInfo->getTime()  - adaptInfo->getTimestep(),
		      adaptInfo->getTimestep());
	
	adaptInfo->setTime(adaptInfo->getTime() - adaptInfo->getTimestep());
	adaptInfo->setTimestep(newTimestep);

	rejected = true;
      }

183
184
185
186
187
188
      if (firstTimestep)
	firstTimestep = false;

      if (fixFirstTimesteps > 0)
	fixFirstTimesteps--;

189
190
191
192
193
194
195
196
197
    } while (rejected);

    rosenbrockStat->acceptTimestep();

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

}