AdaptInstationary.cc 14.3 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: 
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
 * 
 ******************************************************************************/
20
21


22
#include "AdaptInstationary.h"
23
#include "Initfile.h"
24
#include "est/Estimator.h"
25
26
27
28
#include "ProblemIterationInterface.h"
#include "ProblemTimeInterface.h"
#include "Serializer.h"

29
30
#if HAVE_PARALLEL_DOMAIN_AMDIS
#include "parallel/MeshDistributor.h"
31
#ifdef HAVE_PARALLEL_PETSC
32
33
#include <petsc.h>
#endif
Naumann, Andreas's avatar
Naumann, Andreas committed
34
#endif
35

36
37
using namespace std;

38
39
namespace AMDiS {

Thomas Witkowski's avatar
Thomas Witkowski committed
40
  AdaptInstationary::AdaptInstationary(std::string name,
41
				       ProblemIterationInterface *problemStat,  
42
43
44
				       AdaptInfo *info,
				       ProblemTimeInterface *problemInstat,
				       AdaptInfo *initialInfo,
45
				       time_t initialTimestampSet)
46
    : AdaptBase(name, problemStat, info, problemInstat, initialInfo),
47
48
      breakWhenStable(0),
      dbgMode(false)
49
50
51
  {
    FUNCNAME("AdaptInstationary::AdaptInstationary()");

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

55
56
    initConstructor(problemStat, info, initialInfo, initialTimestampSet);
 }
57

Thomas Witkowski's avatar
Thomas Witkowski committed
58

59
60
61
62
63
  AdaptInstationary::AdaptInstationary(std::string name,
				       ProblemIterationInterface &problemStat,  
				       AdaptInfo &info,
				       ProblemTimeInterface &problemInstat,
				       AdaptInfo &initialInfo,
64
				       time_t initialTimestampSet)
65
66
67
68
69
    : AdaptBase(name, &problemStat, &info, &problemInstat, &initialInfo),
      breakWhenStable(0),
      dbgMode(false)
  {

70
    initConstructor(&problemStat, &info, &initialInfo, initialTimestampSet);
71
72
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
73

74
75
76
  void AdaptInstationary::initConstructor(ProblemIterationInterface *problemStat,  
					  AdaptInfo *info,
					  AdaptInfo *initialInfo,
77
					  time_t initialTimestampSet)
78
  {
79
    FUNCNAME("AdaptInstationary::initConstructor()");
Thomas Witkowski's avatar
Thomas Witkowski committed
80
    initialize(name);
81

Thomas Witkowski's avatar
Thomas Witkowski committed
82
    fixedTimestep = (info->getMinTimestep() == info->getMaxTimestep());
83
84
 
    if (initialTimestampSet == 0)
85
      initialTimestamp = time(0);
86
    else
87
88
      initialTimestamp = initialTimestampSet;

89
    // Check if the problem should be deserialized because of the -rs parameter.
Thomas Witkowski's avatar
Thomas Witkowski committed
90
    std::string serializationFilename = "";
91
    Parameters::get("argv->rs", serializationFilename);
92
93
94
95

    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
96
      MSG("Deserialization from file: %s\n", queueSerializationFilename.c_str());
97

98
      std::ifstream in(queueSerializationFilename.c_str() , ios::in);
99
100
101
102
103
104

      // Read the revision number of the AMDiS version which was used to create 
      // the serialization file.
      int revNumber = -1;
      SerUtil::deserialize(in, revNumber);
      
105
106
107
      deserialize(in);
      in.close();

108
109
      info->setIsDeserialized(true);
      initialInfo->setIsDeserialized(true);
110
111
112
113
    } else {
      int readSerialization = 0;
      int readSerializationWithAdaptInfo = 0;

114
115
116
117
      Parameters::get((*problemStat).getName() + "->input->read serialization",
		      readSerialization);
      Parameters::get((*problemStat).getName() + "->input->serialization with adaptinfo",
		      readSerializationWithAdaptInfo);
118
119

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

122
123
	Parameters::get((*problemStat).getName() + "->input->serialization filename", 
			serializationFilename);
124
125
126
	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
127
	std::ifstream in(serializationFilename.c_str());
128
129
130
131
132
133

	// Read the revision number of the AMDiS version which was used to create 
	// the serialization file.
	int revNumber = -1;
	SerUtil::deserialize(in, revNumber);

134
135
136
137
138
139
	deserialize(in);
	in.close();
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
140

141
142
143
144
145
  void AdaptInstationary::explicitTimeStrategy()
  {
    FUNCNAME("AdaptInstationary::explicitTimeStrategy()");

    // estimate before first adaption
Thomas Witkowski's avatar
Thomas Witkowski committed
146
    if (adaptInfo->getTime() <= adaptInfo->getStartTime())
Thomas Witkowski's avatar
Thomas Witkowski committed
147
      problemIteration->oneIteration(adaptInfo, ESTIMATE);
Thomas Witkowski's avatar
Thomas Witkowski committed
148
149


150
    // increment time
Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
151
    adaptInfo->setTime(adaptInfo->getTime() + adaptInfo->getTimestep());
152

Thomas Witkowski's avatar
Thomas Witkowski committed
153
    problemTime->setTime(adaptInfo);
154

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

Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
158
    adaptInfo->setSpaceIteration(0);
159
160
  
    // do the iteration
Thomas Witkowski's avatar
Thomas Witkowski committed
161
162
163
    problemIteration->beginIteration(adaptInfo);
    problemIteration->oneIteration(adaptInfo, FULL_ITERATION);
    problemIteration->endIteration(adaptInfo);
164
    adaptInfo->setLastProcessedTimestep(adaptInfo->getTimestep()); 
165
166
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
167

168
169
170
171
172
  void AdaptInstationary::implicitTimeStrategy()
  {
    FUNCNAME("AdaptInstationary::implicitTimeStrategy()");

    do {
Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
173
      adaptInfo->setTime(adaptInfo->getTime() + adaptInfo->getTimestep());
Thomas Witkowski's avatar
Thomas Witkowski committed
174
      problemTime->setTime(adaptInfo);
175

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

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

Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
181
      adaptInfo->incTimestepIteration();
182

Thomas Witkowski's avatar
Thomas Witkowski committed
183
      if (!fixedTimestep && 
Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
184
	  !adaptInfo->timeToleranceReached() &&
185
	  adaptInfo->getTimestepIteration() <= adaptInfo->getMaxTimestepIteration() &&
Thomas Witkowski's avatar
Thomas Witkowski committed
186
	  !(adaptInfo->getTimestep() <= adaptInfo->getMinTimestep())) {
187
188
189
190
  
	adaptInfo->setTime(adaptInfo->getTime() - adaptInfo->getTimestep());
	adaptInfo->setTimestep(adaptInfo->getTimestep() * timeDelta1);
	continue;
191
      }
192

193

Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
194
      adaptInfo->setSpaceIteration(0);
195
196


Praetorius, Simon's avatar
Praetorius, Simon committed
197
      // === Do space iterations only if the maximum is higher than 0. === 
198

Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
199
      if (adaptInfo->getMaxSpaceIteration() > 0) {
200
    
201
	// === Space iterations. ===
202
	do {
Thomas Witkowski's avatar
Thomas Witkowski committed
203
	  problemIteration->beginIteration(adaptInfo);
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
	      
	  if (dbgMode) {
	    std::cout << "=== ADAPT INFO DEBUG MODE           ===\n";
	    std::cout << "=== in implicitTimeStrategy() ===\n";
	    std::cout << "=== space/time iteration "<< adaptInfo->getSpaceIteration()
	      <<" : "<< adaptInfo->getTimestepIteration() <<" ===\n";
	    adaptInfo->printTimeErrorLowInfo();
	  }
	  
	  Flag adapted = problemIteration->oneIteration(adaptInfo, FULL_ITERATION);
#if HAVE_PARALLEL_DOMAIN_AMDIS
          int recvAllValues = 0;
          int isAdapted = static_cast<bool>(adapted);
          MPI::COMM_WORLD.Allreduce(&isAdapted, &recvAllValues, 1, MPI_INT, MPI_SUM);
          if (recvAllValues) {
#else
          if (adapted) {
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
222
	    if (!fixedTimestep && 
Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
223
		!adaptInfo->timeToleranceReached() &&
Thomas Witkowski's avatar
Thomas Witkowski committed
224
225
		!(adaptInfo->getTimestep() <= adaptInfo->getMinTimestep())) {
	      adaptInfo->setTime(adaptInfo->getTime() - adaptInfo->getTimestep());
Thomas Witkowski's avatar
Thomas Witkowski committed
226
227
	      adaptInfo->setTimestep(adaptInfo->getTimestep() * timeDelta2);
	      problemIteration->endIteration(adaptInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
228
229
230
	      adaptInfo->incSpaceIteration();
	      break;
	    }	
231
232
	  }

Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
233
	  adaptInfo->incSpaceIteration();
Thomas Witkowski's avatar
Thomas Witkowski committed
234
	  problemIteration->endIteration(adaptInfo);
235
	  
Thomas Witkowski's avatar
Thomas Witkowski committed
236
237
	} while (!adaptInfo->spaceToleranceReached() && 
		 adaptInfo->getSpaceIteration() <= adaptInfo->getMaxSpaceIteration());
238
239

      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
240
	problemIteration->endIteration(adaptInfo);
241
242
      }

243

Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
244
    } while(!adaptInfo->timeToleranceReached() &&
Thomas Witkowski's avatar
Thomas Witkowski committed
245
	    !(adaptInfo->getTimestep() <= adaptInfo->getMinTimestep()) && 
Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
246
	    adaptInfo->getTimestepIteration() <= adaptInfo->getMaxTimestepIteration());  
247

248
    adaptInfo->setLastProcessedTimestep(adaptInfo->getTimestep()); 
249
250
251
252

    // 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
253
254
    if (!fixedTimestep && adaptInfo->timeErrorLow()) {
      adaptInfo->setTimestep(adaptInfo->getTimestep() * timeDelta2);
255
256
257
258
259
260
261
      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
262
	if (fixedTimestep)
263
	  std::cout << "   fixedTimestep = true\n";	
Thomas Witkowski's avatar
Thomas Witkowski committed
264
	if (!adaptInfo->timeErrorLow())
Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
265
	  adaptInfo->printTimeErrorLowInfo();
266
      }
267
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
268

269
    // Second, check for decrease of timestep
Thomas Witkowski's avatar
Thomas Witkowski committed
270
    if (!fixedTimestep &&
271
	!adaptInfo->timeToleranceReached() &&
Thomas Witkowski's avatar
Thomas Witkowski committed
272
273
	!(adaptInfo->getTimestep() <= adaptInfo->getMinTimestep()))
	adaptInfo->setTimestep(adaptInfo->getTimestep() * timeDelta1);    
274
275
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
276
277
278

  void AdaptInstationary::simpleAdaptiveTimeStrategy()
  {
279
    FUNCNAME("AdaptInstationary::simpleAdaptiveTimeStrategy()");
Thomas Witkowski's avatar
Thomas Witkowski committed
280
281
282
283
284
285
286
287
288
289
290
291
292
293

    // 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());        
294
295
296
297
298
299

    if (dbgMode) {
	std::cout << "=== ADAPT INFO DEBUG MODE           ===\n";
	std::cout << "=== in simpleAdaptiveTimeStrategy() ===\n";
	adaptInfo->printTimeErrorLowInfo();
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
300
301
302
303
304
305
306
307
308
309
310
311
312
    
    // 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);    
  }


313
314
  void AdaptInstationary::oneTimestep()
  {
315
    FUNCNAME("AdaptInstationary::oneTimestep()");
316

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

319
320
321
322
323
324
325
    switch (strategy) {
    case 0:
      explicitTimeStrategy();
      break;
    case 1:
      implicitTimeStrategy();
      break;
Thomas Witkowski's avatar
Thomas Witkowski committed
326
327
328
    case 2:
      simpleAdaptiveTimeStrategy();
      break;
329
    default:
Thomas Witkowski's avatar
Thomas Witkowski committed
330
      ERROR_EXIT("Unknown strategy = %d!\n", strategy);
331
    }
332

Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
333
    adaptInfo->incTimestepNumber();
334
335
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
336

337
338
339
340
341
  int AdaptInstationary::adapt()
  {
    FUNCNAME("AdaptInstationary::adapt()");
    int errorCode = 0;

Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
342
    TEST_EXIT(adaptInfo->getTimestep() >= adaptInfo->getMinTimestep())
343
      ("timestep < min timestep\n");
Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
344
    TEST_EXIT(adaptInfo->getTimestep() <= adaptInfo->getMaxTimestep())
345
346
      ("timestep > max timestep\n");

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

349
#if HAVE_PARALLEL_DOMAIN_AMDIS
350
    Parallel::MeshDistributor::globalMeshDistributor->initParallelization(); 
351
352
#endif

Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
353
354
    if (adaptInfo->getTimestepNumber() == 0) {
      adaptInfo->setTime(adaptInfo->getStartTime());
Thomas Witkowski's avatar
Thomas Witkowski committed
355
356
      initialAdaptInfo->setStartTime(adaptInfo->getStartTime());
      initialAdaptInfo->setTime(adaptInfo->getStartTime());
357

Thomas Witkowski's avatar
Thomas Witkowski committed
358
      problemTime->setTime(adaptInfo);
359
360

      // initial adaption
Thomas Witkowski's avatar
Thomas Witkowski committed
361
362
      problemTime->solveInitialProblem(initialAdaptInfo);
      problemTime->transferInitialSolution(adaptInfo);
363
364
    }

Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
365
    while (!adaptInfo->reachedEndTime()) {
366
      iterationTimestamp = time(0);
367

Thomas Witkowski's avatar
Thomas Witkowski committed
368
      problemTime->initTimestep(adaptInfo);
369
      oneTimestep();
Thomas Witkowski's avatar
Thomas Witkowski committed
370
      problemTime->closeTimestep(adaptInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
371

372
      if (breakWhenStable && (adaptInfo->getSolverIterations() == 0)) 
373
	break;
374

375
376
377
378
379
380
381
382
383
384
385
386
      // 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;
      }
    }

    return errorCode;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
387

Thomas Witkowski's avatar
Thomas Witkowski committed
388
  void AdaptInstationary::initialize(std::string aName)
389
390
  {
    strategy = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
391
392
393
394
    timeDelta1 = 0.7071;
    timeDelta2 = 1.4142;
    queueRuntime = -1;
    queueSerializationFilename = "__serialized_problem.ser";
395

396
397
398
399
400
    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);
401
    Parameters::get(aName + "->time adaptivity debug mode", dbgMode);
402
403
404
    Parameters::get(aName + "->queue->runtime", queueRuntime);
    Parameters::get(aName + "->queue->serialization filename", 
		    queueSerializationFilename);
405
406
407
408
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    // debug output only on rank 0
    if (MPI::COMM_WORLD.Get_rank() != 0) dbgMode=false;
#endif
409
410
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
411

Thomas Witkowski's avatar
Thomas Witkowski committed
412
  void AdaptInstationary::serialize(std::ostream &out)
413
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
414
    problemIteration->serialize(out);
Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
415
    adaptInfo->serialize(out);
Thomas Witkowski's avatar
Thomas Witkowski committed
416
417
    if (problemTime)
      problemTime->serialize(out);    
418
419
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
420

Thomas Witkowski's avatar
Thomas Witkowski committed
421
  void AdaptInstationary::deserialize(std::istream &in)
422
423
  {
    FUNCNAME("AdaptInstationary::deserialize()");
424
425
426
427

    if (in.fail())
      ERROR_EXIT("File not found for deserialization!\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
428
    problemIteration->deserialize(in);
Thomas Witkowski's avatar
* Bla    
Thomas Witkowski committed
429
    adaptInfo->deserialize(in);
Thomas Witkowski's avatar
Thomas Witkowski committed
430
431
    if (problemTime)
      problemTime->deserialize(in);    
432
433
434
435
436
437
  }


  bool AdaptInstationary::checkQueueRuntime()
  {
    // If there is no time limited runtime queue, there is also nothing to check.
Thomas Witkowski's avatar
Thomas Witkowski committed
438
    if (queueRuntime == -1)
439
      return false;
Thomas Witkowski's avatar
Thomas Witkowski committed
440
    
441
    // Get the current time.
442
    time_t currentTimestamp = time(0);
443
444

    // Update list with the last iteration runtimes.
Thomas Witkowski's avatar
Thomas Witkowski committed
445
    lastIterationsDuration.push(currentTimestamp - iterationTimestamp);
446
    // The list should not contain more than 5 elements. If so, delete the oldest one.
Thomas Witkowski's avatar
Thomas Witkowski committed
447
448
    if (lastIterationsDuration.size() > 5)
      lastIterationsDuration.pop();    
449
450

    // Calculate the avarage of the last iterations.
Thomas Witkowski's avatar
Thomas Witkowski committed
451
    std::queue<int> tmpQueue = lastIterationsDuration;
452
453
    int avrgLastIterations = 0;
    while (!tmpQueue.empty()) {
454
455
	avrgLastIterations += tmpQueue.front();
	tmpQueue.pop();
456
    } 
Thomas Witkowski's avatar
Thomas Witkowski committed
457
    avrgLastIterations /= lastIterationsDuration.size();
458

459
    // Check if there is enough time for a further iteration.
Thomas Witkowski's avatar
Thomas Witkowski committed
460
461
    if (initialTimestamp + queueRuntime - currentTimestamp < avrgLastIterations * 2) {
      std::ofstream out(queueSerializationFilename.c_str());
462
463
464
465
466
467
468
469
470
471
      serialize(out);
      out.close();

      return true;
    }

    return false;
  }

}