AdaptInstationary.cc 13.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
#include "AdaptInstationary.h"
14
#include "Initfile.h"
15
16
17
18
19
#include "Estimator.h"
#include "ProblemIterationInterface.h"
#include "ProblemTimeInterface.h"
#include "Serializer.h"

20
21
22
23
24
#if HAVE_PARALLEL_DOMAIN_AMDIS
#include "parallel/MeshDistributor.h"
#include <petsc.h>
#endif

25
26
namespace AMDiS {

Thomas Witkowski's avatar
Thomas Witkowski committed
27
  AdaptInstationary::AdaptInstationary(std::string name,
28
				       ProblemIterationInterface *problemStat,  
29
30
31
				       AdaptInfo *info,
				       ProblemTimeInterface *problemInstat,
				       AdaptInfo *initialInfo,
32
				       time_t initialTimestampSet)
33
    : AdaptBase(name, problemStat, info, problemInstat, initialInfo),
34
35
      breakWhenStable(0),
      dbgMode(false)
36
37
38
  {
    FUNCNAME("AdaptInstationary::AdaptInstationary()");

39
40
    MSG("You make use of the obsolete constructor AdaptInstationary::AdaptInstationary(...)!\n");
    MSG("Please use the constructor that uses references instead of pointers!\n");
41

42
43
    initConstructor(problemStat, info, initialInfo, initialTimestampSet);
 }
44

Thomas Witkowski's avatar
Thomas Witkowski committed
45

46
47
48
49
50
  AdaptInstationary::AdaptInstationary(std::string name,
				       ProblemIterationInterface &problemStat,  
				       AdaptInfo &info,
				       ProblemTimeInterface &problemInstat,
				       AdaptInfo &initialInfo,
51
				       time_t initialTimestampSet)
52
53
54
55
56
57
    : AdaptBase(name, &problemStat, &info, &problemInstat, &initialInfo),
      breakWhenStable(0),
      dbgMode(false)
  {
    FUNCNAME("AdaptInstationary::AdaptInstationary()");

58
    initConstructor(&problemStat, &info, &initialInfo, initialTimestampSet);
59
60
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
61

62
63
64
  void AdaptInstationary::initConstructor(ProblemIterationInterface *problemStat,  
					  AdaptInfo *info,
					  AdaptInfo *initialInfo,
65
					  time_t initialTimestampSet)
66
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
67
    initialize(name);
68

Thomas Witkowski's avatar
Thomas Witkowski committed
69
    fixedTimestep = (info->getMinTimestep() == info->getMaxTimestep());
70
71
 
    if (initialTimestampSet == 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
72
      initialTimestamp = time(NULL);
73
    else
74
75
      initialTimestamp = initialTimestampSet;

76
    // Check if the problem should be deserialized because of the -rs parameter.
Thomas Witkowski's avatar
Thomas Witkowski committed
77
    std::string serializationFilename = "";
78
    Parameters::get("argv->rs", serializationFilename);
79
80
81
82

    if (serializationFilename.compare("")) {
      // The value of the -rs argument is ignored, because we want to use the 
      // serialization file mentioned in the used init file.
Thomas Witkowski's avatar
Thomas Witkowski committed
83
      MSG("Deserialization from file: %s\n", queueSerializationFilename.c_str());
84

85
      std::ifstream in(queueSerializationFilename.c_str() , ios::in);
86
87
88
      deserialize(in);
      in.close();

89
90
      info->setIsDeserialized(true);
      initialInfo->setIsDeserialized(true);
91
92
93
94
    } else {
      int readSerialization = 0;
      int readSerializationWithAdaptInfo = 0;

95
96
97
98
      Parameters::get((*problemStat).getName() + "->input->read serialization",
		      readSerialization);
      Parameters::get((*problemStat).getName() + "->input->serialization with adaptinfo",
		      readSerializationWithAdaptInfo);
99
100

      if (readSerialization && readSerializationWithAdaptInfo) {
Thomas Witkowski's avatar
Thomas Witkowski committed
101
	std::string serializationFilename = "";
102

103
104
	Parameters::get((*problemStat).getName() + "->input->serialization filename", 
			serializationFilename);
105
106
107
	TEST_EXIT(serializationFilename != "")("no serialization file\n");

	MSG("Deserialization with AdaptInfo from file: %s\n", serializationFilename.c_str());
Thomas Witkowski's avatar
Thomas Witkowski committed
108
	std::ifstream in(serializationFilename.c_str());
109
110
111
112
113
114
	deserialize(in);
	in.close();
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
115

116
117
118
119
120
  void AdaptInstationary::explicitTimeStrategy()
  {
    FUNCNAME("AdaptInstationary::explicitTimeStrategy()");

    // estimate before first adaption
Thomas Witkowski's avatar
Thomas Witkowski committed
121
    if (adaptInfo->getTime() <= adaptInfo->getStartTime())
Thomas Witkowski's avatar
Thomas Witkowski committed
122
      problemIteration->oneIteration(adaptInfo, ESTIMATE);
Thomas Witkowski's avatar
Thomas Witkowski committed
123
124


125
    // increment time
Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
126
    adaptInfo->setTime(adaptInfo->getTime() + adaptInfo->getTimestep());
127

Thomas Witkowski's avatar
Thomas Witkowski committed
128
    problemTime->setTime(adaptInfo);
129

130
    INFO(info, 6)("time = %e, timestep = %e\n",
Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
131
		  adaptInfo->getTime(), adaptInfo->getTimestep());
132

Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
133
    adaptInfo->setSpaceIteration(0);
134
135
  
    // do the iteration
Thomas Witkowski's avatar
Thomas Witkowski committed
136
137
138
    problemIteration->beginIteration(adaptInfo);
    problemIteration->oneIteration(adaptInfo, FULL_ITERATION);
    problemIteration->endIteration(adaptInfo);
139
    adaptInfo->setLastProcessedTimestep(adaptInfo->getTimestep()); 
140
141
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
142

143
144
145
146
147
  void AdaptInstationary::implicitTimeStrategy()
  {
    FUNCNAME("AdaptInstationary::implicitTimeStrategy()");

    do {
Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
148
      adaptInfo->setTime(adaptInfo->getTime() + adaptInfo->getTimestep());
Thomas Witkowski's avatar
Thomas Witkowski committed
149
      problemTime->setTime(adaptInfo);
150

Thomas Witkowski's avatar
Thomas Witkowski committed
151
152
      INFO(info,6)("time = %e, try timestep = %e\n",
		   adaptInfo->getTime(), adaptInfo->getTimestep());
153

Thomas Witkowski's avatar
Thomas Witkowski committed
154
      problemIteration->oneIteration(adaptInfo, NO_ADAPTION);
155

Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
156
      adaptInfo->incTimestepIteration();
157

Thomas Witkowski's avatar
Thomas Witkowski committed
158
      if (!fixedTimestep && 
Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
159
	  !adaptInfo->timeToleranceReached() &&
160
	  adaptInfo->getTimestepIteration() <= adaptInfo->getMaxTimestepIteration() &&
Thomas Witkowski's avatar
Thomas Witkowski committed
161
	  !(adaptInfo->getTimestep() <= adaptInfo->getMinTimestep())) {
162
163
164
165
  
	adaptInfo->setTime(adaptInfo->getTime() - adaptInfo->getTimestep());
	adaptInfo->setTimestep(adaptInfo->getTimestep() * timeDelta1);
	continue;
166
      }
167

168

Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
169
      adaptInfo->setSpaceIteration(0);
170
171


172
      // === Do only space iterations only if the maximum is higher than 0. === 
173

Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
174
      if (adaptInfo->getMaxSpaceIteration() > 0) {
175
    
176
	// === Space iterations. ===
177
	do {
Thomas Witkowski's avatar
Thomas Witkowski committed
178
	  problemIteration->beginIteration(adaptInfo);
179
	  
Thomas Witkowski's avatar
Thomas Witkowski committed
180
181
	  if (problemIteration->oneIteration(adaptInfo, FULL_ITERATION)) {
	    if (!fixedTimestep && 
Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
182
		!adaptInfo->timeToleranceReached() &&
Thomas Witkowski's avatar
Thomas Witkowski committed
183
184
		!(adaptInfo->getTimestep() <= adaptInfo->getMinTimestep())) {
	      adaptInfo->setTime(adaptInfo->getTime() - adaptInfo->getTimestep());
Thomas Witkowski's avatar
Thomas Witkowski committed
185
186
	      adaptInfo->setTimestep(adaptInfo->getTimestep() * timeDelta2);
	      problemIteration->endIteration(adaptInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
187
188
189
	      adaptInfo->incSpaceIteration();
	      break;
	    }	
190
191
	  }

Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
192
	  adaptInfo->incSpaceIteration();
Thomas Witkowski's avatar
Thomas Witkowski committed
193
	  problemIteration->endIteration(adaptInfo);
194
	  
Thomas Witkowski's avatar
Thomas Witkowski committed
195
196
	} while (!adaptInfo->spaceToleranceReached() && 
		 adaptInfo->getSpaceIteration() <= adaptInfo->getMaxSpaceIteration());
197
198

      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
199
	problemIteration->endIteration(adaptInfo);
200
201
      }

202

Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
203
    } while(!adaptInfo->timeToleranceReached() &&
Thomas Witkowski's avatar
Thomas Witkowski committed
204
	    !(adaptInfo->getTimestep() <= adaptInfo->getMinTimestep()) && 
Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
205
	    adaptInfo->getTimestepIteration() <= adaptInfo->getMaxTimestepIteration());  
206

207
    adaptInfo->setLastProcessedTimestep(adaptInfo->getTimestep()); 
208
209
210
211

    // After successful iteration/timestep the timestep will be changed according 
    // adaption rules for next timestep. 
    // First, check for increase of timestep
Thomas Witkowski's avatar
Thomas Witkowski committed
212
213
    if (!fixedTimestep && adaptInfo->timeErrorLow()) {
      adaptInfo->setTimestep(adaptInfo->getTimestep() * timeDelta2);
214
215
216
217
218
219
220
      if (dbgMode) {
	// print information about timestep increase
      }
    } else {
      if (dbgMode) {
	std::cout << "=== ADAPT INFO DEBUG MODE ===\n";
	std::cout << " Do not increase timestep: \n";
Thomas Witkowski's avatar
Thomas Witkowski committed
221
	if (fixedTimestep)
222
	  std::cout << "   fixedTimestep = true\n";	
Thomas Witkowski's avatar
Thomas Witkowski committed
223
	if (!adaptInfo->timeErrorLow())
Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
224
	  adaptInfo->printTimeErrorLowInfo();
225
      }
226
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
227

228
    // Second, check for decrease of timestep
Thomas Witkowski's avatar
Thomas Witkowski committed
229
    if (!fixedTimestep &&
230
	!adaptInfo->timeToleranceReached() &&
Thomas Witkowski's avatar
Thomas Witkowski committed
231
232
	!(adaptInfo->getTimestep() <= adaptInfo->getMinTimestep()))
	adaptInfo->setTimestep(adaptInfo->getTimestep() * timeDelta1);    
233
234
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
235
236
237

  void AdaptInstationary::simpleAdaptiveTimeStrategy()
  {
238
    FUNCNAME("AdaptInstationary::simpleAdaptiveTimeStrategy()");
Thomas Witkowski's avatar
Thomas Witkowski committed
239
240
241
242
243
244
245
246
247
248
249
250
251
252

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

    adaptInfo->setTime(adaptInfo->getTime() + adaptInfo->getTimestep());
    problemTime->setTime(adaptInfo);
    
    INFO(info,6)("time = %e, timestep = %e\n",
		 adaptInfo->getTime(), adaptInfo->getTimestep());
    
    problemIteration->oneIteration(adaptInfo, FULL_ITERATION);

    adaptInfo->setLastProcessedTimestep(adaptInfo->getTimestep());        
253
254
255
256
257
258

    if (dbgMode) {
	std::cout << "=== ADAPT INFO DEBUG MODE           ===\n";
	std::cout << "=== in simpleAdaptiveTimeStrategy() ===\n";
	adaptInfo->printTimeErrorLowInfo();
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
259
260
261
262
263
264
265
266
267
268
269
270
271
    
    // First, check for increase of timestep
    if (!fixedTimestep && adaptInfo->timeErrorLow())
      adaptInfo->setTimestep(adaptInfo->getTimestep() * timeDelta2);
    
    // Second, check for decrease of timestep
    if (!fixedTimestep &&
	!adaptInfo->timeToleranceReached() &&
	!(adaptInfo->getTimestep() <= adaptInfo->getMinTimestep()))
      adaptInfo->setTimestep(adaptInfo->getTimestep() * timeDelta1);    
  }


272
273
  void AdaptInstationary::oneTimestep()
  {
274
    FUNCNAME("AdaptInstationary::oneTimestep()");
275

276
277
    MSG("ONE TIMESTEP!\n");

Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
278
    adaptInfo->setTimestepIteration(0);
279

280
281
282
283
284
285
286
    switch (strategy) {
    case 0:
      explicitTimeStrategy();
      break;
    case 1:
      implicitTimeStrategy();
      break;
Thomas Witkowski's avatar
Thomas Witkowski committed
287
288
289
    case 2:
      simpleAdaptiveTimeStrategy();
      break;
290
    default:
Thomas Witkowski's avatar
Thomas Witkowski committed
291
      ERROR_EXIT("Unknown strategy = %d!\n", strategy);
292
    }
293

Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
294
    adaptInfo->incTimestepNumber();
295
296
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
297

298
299
300
301
302
  int AdaptInstationary::adapt()
  {
    FUNCNAME("AdaptInstationary::adapt()");
    int errorCode = 0;

Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
303
    TEST_EXIT(adaptInfo->getTimestep() >= adaptInfo->getMinTimestep())
304
      ("timestep < min timestep\n");
Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
305
    TEST_EXIT(adaptInfo->getTimestep() <= adaptInfo->getMaxTimestep())
306
307
      ("timestep > max timestep\n");

Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
308
    TEST_EXIT(adaptInfo->getTimestep() > 0)("timestep <= 0!\n");
309

310
311
312
313
#if HAVE_PARALLEL_DOMAIN_AMDIS
    MeshDistributor::globalMeshDistributor->initParallelization(); 
#endif

Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
314
315
    if (adaptInfo->getTimestepNumber() == 0) {
      adaptInfo->setTime(adaptInfo->getStartTime());
Thomas Witkowski's avatar
Thomas Witkowski committed
316
317
      initialAdaptInfo->setStartTime(adaptInfo->getStartTime());
      initialAdaptInfo->setTime(adaptInfo->getStartTime());
318

Thomas Witkowski's avatar
Thomas Witkowski committed
319
      problemTime->setTime(adaptInfo);
320
321

      // initial adaption
Thomas Witkowski's avatar
Thomas Witkowski committed
322
323
      problemTime->solveInitialProblem(initialAdaptInfo);
      problemTime->transferInitialSolution(adaptInfo);
324
325
    }

Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
326
    while (!adaptInfo->reachedEndTime()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
327
      iterationTimestamp = time(NULL);
328

Thomas Witkowski's avatar
Thomas Witkowski committed
329
      problemTime->initTimestep(adaptInfo);
330
      oneTimestep();
Thomas Witkowski's avatar
Thomas Witkowski committed
331
      problemTime->closeTimestep(adaptInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
332

333
      if (breakWhenStable && (adaptInfo->getSolverIterations() == 0)) 
334
	break;
335

336
337
338
339
340
341
342
343
344
      // Check if there is a runtime limitation. If there is a runtime limitation
      // and there is no more time for a next adaption loop, than return the error
      // code for rescheduling the problem and break the adaption loop.
      if (checkQueueRuntime()) {
	errorCode = RescheduleErrorCode;
	break;
      }
    }

345
346
347
348
349
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    MeshDistributor::globalMeshDistributor->exitParallelization();
    PetscFinalize();
#endif

350
351
352
    return errorCode;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
353

Thomas Witkowski's avatar
Thomas Witkowski committed
354
  void AdaptInstationary::initialize(std::string aName)
355
356
357
358
  {
    FUNCNAME("AdaptInstationary::initialize()");

    strategy = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
359
360
361
362
    timeDelta1 = 0.7071;
    timeDelta2 = 1.4142;
    queueRuntime = -1;
    queueSerializationFilename = "__serialized_problem.ser";
363

364
365
366
367
368
    Parameters::get(aName + "->strategy", strategy);
    Parameters::get(aName + "->time delta 1", timeDelta1);
    Parameters::get(aName + "->time delta 2", timeDelta2);
    Parameters::get(aName + "->info", info);
    Parameters::get(aName + "->break when stable", breakWhenStable);
369
    Parameters::get(aName + "->time adaptivity debug mode", dbgMode);
370
371
372
    Parameters::get(aName + "->queue->runtime", queueRuntime);
    Parameters::get(aName + "->queue->serialization filename", 
		    queueSerializationFilename);
373
374
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
375

Thomas Witkowski's avatar
Thomas Witkowski committed
376
  void AdaptInstationary::serialize(std::ostream &out)
377
378
379
  {
    FUNCNAME("AdaptInstationary::serialize()");

380
    SerUtil::serialize(out, amdisRevisionNumber);
Thomas Witkowski's avatar
Thomas Witkowski committed
381
    problemIteration->serialize(out);
Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
382
    adaptInfo->serialize(out);
Thomas Witkowski's avatar
Thomas Witkowski committed
383
384
    if (problemTime)
      problemTime->serialize(out);    
385
386
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
387

Thomas Witkowski's avatar
Thomas Witkowski committed
388
  void AdaptInstationary::deserialize(std::istream &in)
389
390
  {
    FUNCNAME("AdaptInstationary::deserialize()");
391
392
    if(in.fail())
      ERROR_EXIT("File not found for deserialization  \n");
Thomas Witkowski's avatar
Thomas Witkowski committed
393
    problemIteration->deserialize(in);
Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
394
    adaptInfo->deserialize(in);
Thomas Witkowski's avatar
Thomas Witkowski committed
395
396
    if (problemTime)
      problemTime->deserialize(in);    
397
398
399
400
401
402
  }


  bool AdaptInstationary::checkQueueRuntime()
  {
    // If there is no time limited runtime queue, there is also nothing to check.
Thomas Witkowski's avatar
Thomas Witkowski committed
403
    if (queueRuntime == -1) {
404
405
406
407
408
409
      return false;
    }
    // Get the current time.
    time_t currentTimestamp = time(NULL);

    // Update list with the last iteration runtimes.
Thomas Witkowski's avatar
Thomas Witkowski committed
410
    lastIterationsDuration.push(currentTimestamp - iterationTimestamp);
411
    // The list should not contain more than 5 elements. If so, delete the oldest one.
Thomas Witkowski's avatar
Thomas Witkowski committed
412
413
    if (lastIterationsDuration.size() > 5)
      lastIterationsDuration.pop();    
414
415

    // Calculate the avarage of the last iterations.
Thomas Witkowski's avatar
Thomas Witkowski committed
416
    std::queue<int> tmpQueue = lastIterationsDuration;
417
418
    int avrgLastIterations = 0;
    while (!tmpQueue.empty()) {
419
420
	avrgLastIterations += tmpQueue.front();
	tmpQueue.pop();
421
    } 
Thomas Witkowski's avatar
Thomas Witkowski committed
422
    avrgLastIterations /= lastIterationsDuration.size();
423

424
    // Check if there is enough time for a further iteration.
Thomas Witkowski's avatar
Thomas Witkowski committed
425
426
    if (initialTimestamp + queueRuntime - currentTimestamp < avrgLastIterations * 2) {
      std::ofstream out(queueSerializationFilename.c_str());
427
428
429
430
431
432
433
434
435
436
      serialize(out);
      out.close();

      return true;
    }

    return false;
  }

}