ProblemStat.cc 51.7 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


Thomas Witkowski's avatar
Thomas Witkowski committed
22
#include <sstream>
Thomas Witkowski's avatar
Thomas Witkowski committed
23
#include <boost/lexical_cast.hpp>
24
25

#include "ProblemStat.h"
26
27
28
29
30
31
32
33
#include "Serializer.h"
#include "AbstractFunction.h"
#include "Operator.h"
#include "SystemVector.h"
#include "DOFMatrix.h"
#include "FiniteElemSpace.h"
#include "Marker.h"
#include "AdaptInfo.h"
34
#include "io/FileWriter.h"
35
36
#include "CoarseningManager.h"
#include "RefinementManager.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
37
#include "DualTraverse.h"
38
#include "Mesh.h"
39
#include "solver/LinearSolver.h"
40
41
42
43
#include "DirichletBC.h"
#include "RobinBC.h"
#include "PeriodicBC.h"
#include "Lagrange.h"
44
#include "Bubble.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
45
#include "Flag.h"
46
#include "est/Estimator.h"
47
48
#include "io/VtkWriter.h"
#include "io/ValueReader.h"
49
#include "ProblemStatDbg.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
50
#include "Debug.h"
51
52
53

namespace AMDiS {

54
  using namespace std;
55
56
  using boost::lexical_cast;

57
58
59
60
61
62
63
  ProblemStatSeq::ProblemStatSeq(string nameStr,
				 ProblemIterationInterface *problemIteration)
    : StandardProblemIteration(this),
      name(nameStr),
      nComponents(-1),
      nMeshes(0),
      traverseInfo(0),
64
65
66
67
      solver(nullptr),
      solution(nullptr),
      rhs(nullptr),
      systemMatrix(nullptr),
68
      useGetBound(true),
69
70
      refinementManager(nullptr),
      coarseningManager(nullptr),
71
72
73
74
      info(10),
      deserialized(false),
      computeExactError(false),
      boundaryConditionSet(false),
75
76
77
      writeAsmInfo(false),
      solutionTime(0.0),
      buildTime(0.0)
78
79
  {
    Parameters::get(name + "->components", nComponents);
80
81
    TEST_EXIT(nComponents > 0)("No value set for parameter \"%s->components\"!\n",
			       name.c_str());    
82
83
    estimator.resize(nComponents, nullptr);
    marker.resize(nComponents, nullptr);
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
    
    assembleMatrixOnlyOnce.resize(nComponents);
    assembledMatrix.resize(nComponents);
    for (int i = 0; i < nComponents; i++) {
      assembleMatrixOnlyOnce[i].resize(nComponents);
      assembledMatrix[i].resize(nComponents);
      for (int j = 0; j < nComponents; j++) {
	assembleMatrixOnlyOnce[i][j] = false;
	assembledMatrix[i][j] = false;
      }
    }
    
    exactSolutionFcts.resize(nComponents);


    // === Initialize name of components. ===

    componentNames.resize(nComponents, "");
    for (int i = 0; i < nComponents; i++)
      componentNames[i] = "solution[" + lexical_cast<string>(i) + "]";

    Parameters::get(name + "->name", componentNames);
106
107
    componentNames.resize(nComponents);

108
109
110
    for (int i = 0; i < nComponents; i++)
      Parameters::get(name + "->name[" + lexical_cast<string>(i) + "]",
		      componentNames[i]);
Thomas Witkowski's avatar
Thomas Witkowski committed
111
112

    Parameters::get("debug->write asm info", writeAsmInfo);
113
114
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
115
116
117

  ProblemStatSeq::~ProblemStatSeq()
  {
118
119
    if (rhs)
      delete rhs;
120
    rhs = nullptr;
121
122
    if (solution)
      delete solution;
123
    solution = nullptr;
Thomas Witkowski's avatar
Thomas Witkowski committed
124

125
126
127
128
129
    if (systemMatrix) {    
      for (int i = 0; i < nComponents; i++)
	for (int j = 0; j < nComponents; j++)
	  if ((*systemMatrix)[i][j]) {
	    delete (*systemMatrix)[i][j];
130
	    (*systemMatrix)[i][j] = nullptr;
131
132
133
	  }

      delete systemMatrix;
134
      systemMatrix = nullptr;
135
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
136
137

    for (unsigned int i = 0; i < meshes.size(); i++)
138
139
      if (meshes[i]) {
// 	delete meshes[i];
140
// 	meshes[i] = nullptr;
141
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
142
143
    
    for (unsigned int i = 0; i < estimator.size(); i++)
144
      if (estimator[i]) {
145
	delete estimator[i];
146
	estimator[i] = nullptr;
147
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
148
149

    for (unsigned int i = 0; i < marker.size(); i++)
150
      if (marker[i]) {
151
	delete marker[i];
152
	marker[i] = nullptr;
153
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
154
155
156
  }


157
158
159
  void ProblemStatSeq::initialize(Flag initFlag,
				  ProblemStatSeq *adoptProblem,
				  Flag adoptFlag)
160
  {
161
    FUNCNAME("ProblemStat::initialize()");
162

163
    // === create meshes ===
Thomas Witkowski's avatar
Thomas Witkowski committed
164
    if (meshes.size() != 0) { 
165
166
167
      WARNING("meshes already created\n");
    } else {
      if (initFlag.isSet(CREATE_MESH) || 
Thomas Witkowski's avatar
Thomas Witkowski committed
168
	  (!adoptFlag.isSet(INIT_MESH) &&
Thomas Witkowski's avatar
Thomas Witkowski committed
169
	   (initFlag.isSet(INIT_SYSTEM) || initFlag.isSet(INIT_FE_SPACE))))
170
	createMesh();      
Thomas Witkowski's avatar
Thomas Witkowski committed
171

172
      if (adoptProblem && 
Thomas Witkowski's avatar
Thomas Witkowski committed
173
	  (adoptFlag.isSet(INIT_MESH) ||
174
175
	   adoptFlag.isSet(INIT_SYSTEM) ||
	   adoptFlag.isSet(INIT_FE_SPACE))) {
Thomas Witkowski's avatar
Thomas Witkowski committed
176
	meshes = adoptProblem->getMeshes();
Thomas Witkowski's avatar
Thomas Witkowski committed
177
	componentMeshes = adoptProblem->componentMeshes;
178
179
180

	// If the adopt problem has fewer components than this problem, but only one
	// mesh for all component, than scal up the componentMeshes array.
Thomas Witkowski's avatar
Thomas Witkowski committed
181
	if (adoptProblem->getNumComponents() < nComponents) {
Thomas Witkowski's avatar
Thomas Witkowski committed
182
	  TEST_EXIT(meshes.size() == 1)("Daran muss ich noch arbeiten!\n");
183
	  
Thomas Witkowski's avatar
Thomas Witkowski committed
184
	  componentMeshes.resize(nComponents);
185
	  for (int i = adoptProblem->getNumComponents(); i < nComponents; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
186
	    componentMeshes[i] = componentMeshes[0];
187
	}
188
189
190
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
191
    if (meshes.size() == 0) 
192
193
      WARNING("no mesh created\n");

194
    // === create refinement/corasening-manager ===
195
    if (refinementManager != nullptr && coarseningManager != nullptr) { 
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
      WARNING("refinement-/coarseningmanager already created\n");
    } else {
      if (initFlag.isSet(CREATE_MESH) || 
          (!adoptFlag.isSet(INIT_MESH) &&
           (initFlag.isSet(INIT_SYSTEM) || initFlag.isSet(INIT_FE_SPACE))))
        createRefCoarseManager();

      if (adoptProblem && 
          (adoptFlag.isSet(INIT_MESH) ||
           adoptFlag.isSet(INIT_SYSTEM) ||
           adoptFlag.isSet(INIT_FE_SPACE))) {
        refinementManager = adoptProblem->refinementManager;
        coarseningManager = adoptProblem->coarseningManager;
      }
    }

212
    if (refinementManager == nullptr || coarseningManager == nullptr) 
213
214
      WARNING("no refinement-/coarseningmanager created\n");

215
    // === create fespace ===
216
    if (feSpaces.size() != 0) {
217
218
219
      WARNING("feSpaces already created\n");
    } else {
      if (initFlag.isSet(INIT_FE_SPACE) || 
220
	  (initFlag.isSet(INIT_SYSTEM) && !adoptFlag.isSet(INIT_FE_SPACE)))
221
	createFeSpace(nullptr);
222

223
224
      if (adoptProblem &&
	  (adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) {
225
	feSpaces = adoptProblem->getFeSpaces();
226
	componentSpaces = adoptProblem->componentSpaces;
Thomas Witkowski's avatar
Thomas Witkowski committed
227
	traverseInfo = adoptProblem->traverseInfo;
228
229
230

	// If the adopt problem has fewer components than this problem, but only one
	// fe space for all component, than scal up the componentSpaces array.
Thomas Witkowski's avatar
Thomas Witkowski committed
231
232
	if (adoptProblem->getNumComponents() < nComponents) {
	  TEST_EXIT(feSpaces.size() == 1)("Daran muss ich noch arbeiten!\n");
233
	  
Thomas Witkowski's avatar
Thomas Witkowski committed
234
	  componentSpaces.resize(nComponents);
Thomas Witkowski's avatar
Thomas Witkowski committed
235
	  for (int i = adoptProblem->getNumComponents(); i < nComponents; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
236
	    componentSpaces[i] = componentSpaces[0];
237
238
	}

239
240
241
      }
    }

242
    if (feSpaces.size() == 0) 
243
244
245
      WARNING("no feSpace created\n");

    // === create system ===
246
    if (initFlag.isSet(INIT_SYSTEM)) 
247
      createMatricesAndVectors();
248
    
249
    if (adoptProblem && adoptFlag.isSet(INIT_SYSTEM)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
250
      solution = adoptProblem->getSolution();
251
      rhs = adoptProblem->getRhs();
Thomas Witkowski's avatar
Thomas Witkowski committed
252
      systemMatrix = adoptProblem->getSystemMatrix();
253
254
255
    }

    // === create solver ===
Thomas Witkowski's avatar
Thomas Witkowski committed
256
    if (solver) {
257
258
      WARNING("solver already created\n");
    } else {
259
      if (initFlag.isSet(INIT_SOLVER))
260
	createSolver();
261
      
262
      if (adoptProblem && adoptFlag.isSet(INIT_SOLVER)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
263
264
	TEST_EXIT(!solver)("solver already created\n");
	solver = adoptProblem->getSolver();
265
266
267
      }
    }

268
// #ifndef HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
269
    if (!solver) 
270
      WARNING("no solver created\n");
271
// #endif
272
273

    // === create estimator ===
Thomas Witkowski's avatar
Thomas Witkowski committed
274
    if (initFlag.isSet(INIT_ESTIMATOR))
275
      createEstimator();
Thomas Witkowski's avatar
Thomas Witkowski committed
276
277

    if (adoptProblem && adoptFlag.isSet(INIT_ESTIMATOR))
278
      estimator = adoptProblem->getEstimators();
279
280

    // === create marker ===
Thomas Witkowski's avatar
Thomas Witkowski committed
281
    if (initFlag.isSet(INIT_MARKER))
282
      createMarker();
Thomas Witkowski's avatar
Thomas Witkowski committed
283

Thomas Witkowski's avatar
Thomas Witkowski committed
284
    if (adoptProblem && adoptFlag.isSet(INIT_MARKER))
285
      marker = adoptProblem->getMarkers();
286
287

    // === create file writer ===
Thomas Witkowski's avatar
Thomas Witkowski committed
288
    if (initFlag.isSet(INIT_FILEWRITER))
289
290
291
292
293
      createFileWriter();
    
    // === read serialization and init mesh ===
    
    // There are two possiblities where the user can define a serialization
294
295
296
297
    // to be read from disk. Either by providing the parameter -rs when 
    // executing the program or in the init file. The -rs parameter is always
    // checked first, because it can be added automatically when rescheduling
    // the program before timeout of the runqueue.
298
299

    int readSerialization = 0;
300
301
    string serializationFilename = "";
    Parameters::get("argv->rs", serializationFilename);
302
303
304
305
306
307

    // If the parameter -rs is set, we do nothing here, because the problem will be
    // deserialized in the constructor of a following AdaptInstationary initialization.
    if (!serializationFilename.compare("")) {
      int readSerializationWithAdaptInfo = 0;

308
309
310
      Parameters::get(name + "->input->read serialization", readSerialization);
      Parameters::get(name + "->input->serialization with adaptinfo",
		      readSerializationWithAdaptInfo);
311
312
313
314
315
316

      // The serialization file is only read, if the adaptInfo part should not be used.
      // If the adaptInfo part should be also read, the serialization file will be read
      // in the constructor of the AdaptInstationary problem, because we do not have here
      // the adaptInfo object.
      if (readSerialization && !readSerializationWithAdaptInfo) {
317
318
	Parameters::get(name + "->input->serialization filename", 
			serializationFilename);
319
320
	TEST_EXIT(serializationFilename != "")("no serialization file\n");

321
322
	// If AMDiS is compiled for parallel computations, the deserialization is
	// controlled by the parallel problem object.
323
#ifndef HAVE_PARALLEL_DOMAIN_AMDIS
324
	MSG("Deserialization from file: %s\n", serializationFilename.c_str());
325
	ifstream in(serializationFilename.c_str());
326
327
328
329
330
331

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

332
333
	deserialize(in);
	in.close();
334
#endif
335
336

	deserialized = true;
337
      } else {
338
	int globalRefinements = 0;
339

340
	// If AMDiS is compiled for parallel computations, the global refinements are
341
342
343
	// ignored here. Later, each rank will add the global refinements to its 
	// private mesh.
#ifndef HAVE_PARALLEL_DOMAIN_AMDIS
344
345
	Parameters::get(meshes[0]->getName() + "->global refinements", 
			globalRefinements);
346
#endif
347

348
	bool initMesh = initFlag.isSet(INIT_MESH);
349

Thomas Witkowski's avatar
Thomas Witkowski committed
350
351
	// Initialize the meshes if there is no serialization file.
	for (int i = 0; i < static_cast<int>(meshes.size()); i++)
352
	  if (initMesh && meshes[i] && !(meshes[i]->isInitialized()))
Thomas Witkowski's avatar
Thomas Witkowski committed
353
	    meshes[i]->initialize();	    	
354
355

	// === read value file and use it for the mesh values ===
356
357
	string valueFilename("");
	Parameters::get(meshes[0]->getName() + "->value file name", valueFilename); 
358
	if (valueFilename.length())
359
	  io::ValueReader::readValue(valueFilename,
360
361
362
363
364
				 meshes[0],
				 solution->getDOFVector(0),
				 meshes[0]->getMacroFileInfo());

	// === do global refinements ===
365

Thomas Witkowski's avatar
Thomas Witkowski committed
366
	for (unsigned int i = 0; i < meshes.size(); i++)
367
	  if (initMesh && meshes[i])
Thomas Witkowski's avatar
Thomas Witkowski committed
368
	    refinementManager->globalRefine(meshes[i], globalRefinements);	
369
370
371
372
373
374
      }
    }

    doOtherStuff();
  }

375

376
  void ProblemStatSeq::createMesh() 
377
  {
378
    FUNCNAME("ProblemStat::createMesh()");
379

Thomas Witkowski's avatar
Thomas Witkowski committed
380
    componentMeshes.resize(nComponents);
381
382
383
384
385
386
    map<int, Mesh*> meshForRefinementSet;

    string meshName("");
    Parameters::get(name + "->mesh", meshName);
    TEST_EXIT(meshName != "")("No mesh name specified for \"%s->mesh\"!\n", 
			      name.c_str());
387
388

    int dim = 0;
389
390
391
    Parameters::get(name + "->dim", dim);
    TEST_EXIT(dim)("No problem dimension specified for \"%s->dim\"!\n",
		   name.c_str());
392

393
    for (int i = 0; i < nComponents; i++) {
394
      int refSet = -1;
395
396
      Parameters::get(name + "->refinement set[" + lexical_cast<string>(i) + "]", 
		      refSet);
397
      if (refSet < 0)
398
	refSet = 0;
399
      
400
      if (meshForRefinementSet[refSet] == nullptr) {
Thomas Witkowski's avatar
Thomas Witkowski committed
401
	Mesh *newMesh = new Mesh(meshName, dim);
402
	meshForRefinementSet[refSet] = newMesh;
Thomas Witkowski's avatar
Thomas Witkowski committed
403
404
	meshes.push_back(newMesh);
	nMeshes++;
405
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
406
      componentMeshes[i] = meshForRefinementSet[refSet];
407
    }
408
409
410
411
412
413
414
415
416
417
418
  }


  void ProblemStatSeq::createRefCoarseManager() 
  {
    FUNCNAME("ProblemStat::createRefCoarseManager()");

    int dim = 0;
    Parameters::get(name + "->dim", dim);
    TEST_EXIT(dim)("No problem dimension specified for \"%s->dim\"!\n",
                   name.c_str());
419
420

    switch (dim) {
421
    case 1:
Thomas Witkowski's avatar
Thomas Witkowski committed
422
      refinementManager = new RefinementManager1d();
423
      coarseningManager = new CoarseningManager1d();
424
425
      break;
    case 2:
Thomas Witkowski's avatar
Thomas Witkowski committed
426
      refinementManager = new RefinementManager2d();
427
      coarseningManager = new CoarseningManager2d();
428
429
      break;
    case 3:
Thomas Witkowski's avatar
Thomas Witkowski committed
430
      refinementManager = new RefinementManager3d();
431
      coarseningManager = new CoarseningManager3d();
432
433
434
435
436
437
      break;
    default:
      ERROR_EXIT("invalid dim!\n");
    }
  }

438

439
  void ProblemStatSeq::createFeSpace(DOFAdmin *admin)
440
  {
441
    FUNCNAME("ProblemStat::createFeSpace()");
442

443
    map<pair<Mesh*, string>, FiniteElemSpace*> feSpaceMap;
444
    int dim = -1;
445
    Parameters::get(name + "->dim", dim);
446
    TEST_EXIT(dim != -1)("no problem dimension specified!\n");
447

448
    componentSpaces.resize(nComponents, nullptr);
Thomas Witkowski's avatar
Thomas Witkowski committed
449
    traverseInfo.resize(nComponents);
450

451
    for (int i = 0; i < nComponents; i++) {
452
      TEST_EXIT(componentSpaces[i] == nullptr)("feSpace already created\n");
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
      string componentString = "[" + boost::lexical_cast<string>(i) + "]";
      
      string feSpaceName = "";
      string initFileStr = name + "->feSpace" + componentString;
      Parameters::get(initFileStr, feSpaceName);
      
      // synonym for "feSpace"
      if (feSpaceName.size() == 0) {
	initFileStr = name + "->finite element space" + componentString;
	Parameters::get(initFileStr, feSpaceName);
      }
      
      // for backward compatibility also accept the old syntax
      if (feSpaceName.size() == 0) {
	int degree = 1;
	initFileStr = name + "->polynomial degree" + componentString;
	Parameters::get(initFileStr, degree);
	TEST_EXIT(degree > 0)
	  ("Poynomial degree in component %d must be larger than zero!\n", i);
	  
	feSpaceName = "Lagrange" + boost::lexical_cast<string>(degree);
      }
475

476
477
478
479
480
481
482
      if (feSpaceMap[pair<Mesh*, string>(componentMeshes[i], feSpaceName)] == NULL) {      
	BasisFunctionCreator *basisFctCreator = 
	dynamic_cast<BasisFunctionCreator*>(CreatorMap<BasisFunction>::getCreator(feSpaceName, initFileStr));
	TEST_EXIT(basisFctCreator)
	  ("No valid basisfunction type found in parameter \"%s\"\n", initFileStr.c_str());
	basisFctCreator->setDim(dim);
	
483
	FiniteElemSpace *newFeSpace = 
484
485
486
487
	  FiniteElemSpace::provideFeSpace(admin, basisFctCreator->create(),
					  componentMeshes[i], "FeSpace" + componentString + " (" + feSpaceName + ")");

	feSpaceMap[pair<Mesh*, string>(componentMeshes[i], feSpaceName)] = newFeSpace;
488
	feSpaces.push_back(newFeSpace);
489
      }
490

491
      componentSpaces[i] = feSpaceMap[pair<Mesh*, string>(componentMeshes[i], feSpaceName)];
Thomas Witkowski's avatar
Thomas Witkowski committed
492
493
494
    }

    for (int i = 0; i < nComponents; i++) {
495
      for (int j = 0; j < nComponents; j++)
496
	traverseInfo.getMatrix(i, j).setFeSpace(componentSpaces[i], componentSpaces[j]);
497
      
498
      traverseInfo.getVector(i).setFeSpace(componentSpaces[i]);
499
500
501
    }

    // create dof admin for vertex dofs if neccessary
Thomas Witkowski's avatar
Thomas Witkowski committed
502
    for (int i = 0; i < static_cast<int>(meshes.size()); i++) {
503
      if (meshes[i]->getNumberOfDofs(VERTEX) == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
504
	DimVec<int> ln_dof(meshes[i]->getDim(), DEFAULT_VALUE, 0);
505
	ln_dof[VERTEX] = 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
506
	meshes[i]->createDOFAdmin("vertex dofs", ln_dof);      
507
508
509
510
      }
    }
  }

511

512
  void ProblemStatSeq::createMatricesAndVectors()
513
514
515
  {
    // === create vectors and system matrix ===

Thomas Witkowski's avatar
Thomas Witkowski committed
516
    systemMatrix = new Matrix<DOFMatrix*>(nComponents, nComponents);
517
    systemMatrix->set(nullptr);
Thomas Witkowski's avatar
Thomas Witkowski committed
518
519
    rhs = new SystemVector("rhs", componentSpaces, nComponents);
    solution = new SystemVector("solution", componentSpaces, nComponents);
520

Thomas Witkowski's avatar
Thomas Witkowski committed
521
    for (int i = 0; i < nComponents; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
522
      (*systemMatrix)[i][i] = new DOFMatrix(componentSpaces[i], 
Thomas Witkowski's avatar
Thomas Witkowski committed
523
					    componentSpaces[i], "A_ii");
Thomas Witkowski's avatar
Thomas Witkowski committed
524
      (*systemMatrix)[i][i]->setCoupleMatrix(false);
Thomas Witkowski's avatar
Thomas Witkowski committed
525

526
      //set parallel synchronization later in ParalleProblemStat
527
      rhs->setDOFVector(i, new DOFVector<double>(componentSpaces[i], 
528
						 "rhs[" + lexical_cast<string>(i) + "]", false));
529

530
      //set parallel synchronization later in ParalleProblemStat
Thomas Witkowski's avatar
Thomas Witkowski committed
531
      solution->setDOFVector(i, new DOFVector<double>(componentSpaces[i],
532
      						      componentNames[i], false));
Thomas Witkowski's avatar
Thomas Witkowski committed
533
534
      solution->getDOFVector(i)->setCoarsenOperation(COARSE_INTERPOL);
      solution->getDOFVector(i)->set(0.0);
535
536
537
    }
  }

538

539
  void ProblemStatSeq::createSolver()
540
  {
541
    FUNCNAME("ProblemStat::createSolver()");
542
    
Praetorius, Simon's avatar
Praetorius, Simon committed
543
    // definition of standard-backends
544
545
546
547
#if defined HAVE_PARALLEL_PETSC
    string backend("p_petsc");
#elif defined HAVE_PARALLEL_MTL
    string backend("p_mtl");
548
#elif defined HAVE_SEQ_PETSC
549
550
551
552
553
554
555
556
    string backend("petsc");
#else
    string backend("mtl");
#endif
    
    // === read backend-name ===
    string initFileStr = name + "->solver";    
    Parameters::get(initFileStr + "->backend", backend);
557

558
    // === read solver-name ===
559
    string solverType("0");
560
    Parameters::get(initFileStr, solverType);
561
562
563
564
565
566
567
    
    if (backend != "0" && backend != "no" && backend != "")
      solverType = backend + "_" + solverType;
    
    // === create solver ===
    LinearSolverCreator *solverCreator = 
      dynamic_cast<LinearSolverCreator*>(CreatorMap<LinearSolver>::getCreator(solverType, initFileStr));
568
569
    TEST_EXIT(solverCreator)
      ("No valid solver type found in parameter \"%s\"\n", initFileStr.c_str());
570
    solverCreator->setName(initFileStr);
Thomas Witkowski's avatar
Thomas Witkowski committed
571
    solver = solverCreator->create();
572
573
  }

574

575
  void ProblemStatSeq::createEstimator()
576
  {
577
    FUNCNAME("ProblemStat::createEstimator()");
578
579

    // create and set leaf data prototype
Thomas Witkowski's avatar
Thomas Witkowski committed
580
    for (unsigned int i = 0; i < meshes.size(); i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
581
      meshes[i]->setElementDataPrototype
Thomas Witkowski's avatar
Thomas Witkowski committed
582
	(new LeafDataEstimatableVec(new LeafDataCoarsenableVec));
583

Thomas Witkowski's avatar
Thomas Witkowski committed
584
    for (int i = 0; i < nComponents; i++) {
585
      TEST_EXIT(estimator[i] == nullptr)("estimator already created\n");
586
587
      string estName = 
	name + "->estimator[" + boost::lexical_cast<string>(i) + "]";
588
589

      // === create estimator ===
590
591
      string estimatorType("0");
      Parameters::get(estName, estimatorType);
592
      EstimatorCreator *estimatorCreator = 
593
	dynamic_cast<EstimatorCreator*>(CreatorMap<Estimator>::getCreator(estimatorType, estName));
Thomas Witkowski's avatar
Thomas Witkowski committed
594
      if (estimatorCreator) {
595
596
	estimatorCreator->setName(estName);
	estimatorCreator->setRow(i);
597
	estimatorCreator->setSolution(solution->getDOFVector(i));
Thomas Witkowski's avatar
Thomas Witkowski committed
598
	estimator[i] = estimatorCreator->create();
599
600
601
      }


Thomas Witkowski's avatar
Thomas Witkowski committed
602
603
604
605
      if (estimator[i])
	for (int j = 0; j < nComponents; j++)
	  estimator[i]->addSystem((*systemMatrix)[i][j], 
				  solution->getDOFVector(j), 
Praetorius, Simon's avatar
Praetorius, Simon committed
606
				  rhs->getDOFVector(j));      // TODO: hier eventuell (i) statt (j) ???
607
608
609
    }
  }

610

611
  void ProblemStatSeq::createMarker()
612
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
613
    int nMarkersCreated = 0;
614

615
    for (int i = 0; i < nComponents; i++) {
616
      marker[i] = Marker::createMarker
617
	(name + "->marker[" + boost::lexical_cast<string>(i) + "]", i);
618

619
      if (marker[i]) {
Thomas Witkowski's avatar
Thomas Witkowski committed
620
621
622
623
	nMarkersCreated++;

	// If there is more than one marker, and all components are defined
	// on the same mesh, the maximum marking has to be enabled.
624
 	if (nMarkersCreated > 1 && nMeshes == 1)
Thomas Witkowski's avatar
Thomas Witkowski committed
625
 	  marker[i]->setMaximumMarking(true);
626
627
628
629
      }
    }
  }

630

631
  void ProblemStatSeq::createFileWriter()
632
  {
633
    FUNCNAME("ProblemStat::createFileWriter()");
634
635
  
    // Create one filewriter for all components of the problem
636
637
638
    string numberedName  = name + "->output";
    string filename = "";
    Parameters::get(numberedName + "->filename", filename);
639
640

    if (filename != "") {
641
      vector< DOFVector<double>* > solutionList(nComponents);
642

643
      for (int i = 0; i < nComponents; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
644
	TEST_EXIT(componentMeshes[0] == componentMeshes[i])
645
646
	  ("All Meshes have to be equal to write a vector file.\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
647
	solutionList[i] = solution->getDOFVector(i);
648
649
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
650
      fileWriters.push_back(new FileWriter(numberedName,
651
652
					   componentMeshes[0],
					   solutionList));
653
654
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
655
    // Create own filewriters for each component of the problem
656
    for (int i = 0; i < nComponents; i++) {
657
      numberedName = name + "->output[" + boost::lexical_cast<string>(i) + "]";
658
      filename = "";
659
      Parameters::get(numberedName + "->filename", filename);
660

661
      if (filename != "")
Thomas Witkowski's avatar
Thomas Witkowski committed
662
	fileWriters.push_back(new FileWriter(numberedName, 
663
					     componentMeshes[i], 
664
					     solution->getDOFVector(i)));      
665
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
666
    
667
    // create a filewrite for groups of components to write vector-valued output
Praetorius, Simon's avatar
Praetorius, Simon committed
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
    int nVectors = 0;
    Parameters::get(name + "->output->num vectors", nVectors);
    if (nVectors > 0) {
      for (int j = 0; j < nVectors; j++) {
	numberedName = name + "->output->vector[" + boost::lexical_cast<string>(j) + "]";
	
	filename = "";
	Parameters::get(numberedName + "->filename", filename);
	std::string componentName = "";
	Parameters::get(numberedName + "->name", componentName);
	std::vector<std::string> names; 
	if (name != "")
	  names.push_back(componentName);
	std::vector<int> comp;
	Parameters::get(numberedName + "->components", comp);
	
	if (filename != "" && comp.size() > 0) {
	  // Create own filewriters for each component of the problem
	  std::vector<DOFVector<double>*> vectors;
	  for (size_t i = 0; i < comp.size(); i++)
	    vectors.push_back(solution->getDOFVector(comp[i]));
	  
	  fileWriters.push_back(new FileWriter(numberedName, 
					      componentMeshes[comp[0]], 
					      vectors,
					      names
					      ));      
	  
	}
      }
    }
699
700

    int writeSerialization = 0;
701
    Parameters::get(name + "->output->write serialization", writeSerialization);
702
    if (writeSerialization)
703
      fileWriters.push_back(new Serializer<ProblemStatSeq>(this));
704
705
  }

706

707
  void ProblemStatSeq::doOtherStuff()
708
  {}
709

710

711
712
713
  void ProblemStatSeq::solve(AdaptInfo *adaptInfo, 	       
			     bool createMatrixData,
			     bool storeMatrixData)
714
715
716
  {
    FUNCNAME("Problem::solve()");

Thomas Witkowski's avatar
Thomas Witkowski committed
717
    if (!solver) {
718
719
720
721
      WARNING("no solver\n");
      return;
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
722
    Timer t;
723
724
    solver->solveSystem(solverMatrix, *solution, *rhs, 
			createMatrixData, storeMatrixData);
725

Thomas Witkowski's avatar
Thomas Witkowski committed
726
    INFO(info, 8)("solution of discrete system needed %.5f seconds\n", 
Praetorius, Simon's avatar
Praetorius, Simon committed
727
		  t.elapsed());
728
    solutionTime = t.elapsed();
729

730
    adaptInfo->setSolverIterations(solver->getIterations());
Thomas Witkowski's avatar
Thomas Witkowski committed
731
732
733
    adaptInfo->setMaxSolverIterations(solver->getMaxIterations());
    adaptInfo->setSolverTolerance(solver->getTolerance());
    adaptInfo->setSolverResidual(solver->getResidual());
734
735
  }

736

737
  void ProblemStatSeq::estimate(AdaptInfo *adaptInfo) 
738
  {
739
    FUNCNAME("ProblemStat::estimate()");
740

Praetorius, Simon's avatar
Praetorius, Simon committed
741
    Timer t;
742

743
    if (computeExactError) {
Thomas Witkowski's avatar
Thomas Witkowski committed
744
      computeError(adaptInfo);
745
746
    } else {
      for (int i = 0; i < nComponents; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
747
	Estimator *scalEstimator = estimator[i];
748
749
	
	if (scalEstimator) {
750
751
	  traverseInfo.updateStatus();
	  scalEstimator->setTraverseInfo(traverseInfo);
752
	  scalEstimator->estimate(adaptInfo->getTimestep());
753

754
755
	  adaptInfo->setEstSum(scalEstimator->getErrorSum(), i);
	  adaptInfo->setEstMax(scalEstimator->getErrorMax(), i);
756
757
758
759
760

	  if (adaptInfo->getRosenbrockMode() == false) {
	    adaptInfo->setTimeEstSum(scalEstimator->getTimeEst(), i);
	    adaptInfo->setTimeEstMax(scalEstimator->getTimeEstMax(), i);
	  }
761
	}
762
763
764
      }
    }

765
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
766
    MPI::COMM_WORLD.Barrier();
767
#endif
Praetorius, Simon's avatar
Praetorius, Simon committed
768
769
    INFO(info, 8)("estimation of the error needed %.5f seconds\n", 
		  t.elapsed());
770
771
  }

772

773
  Flag ProblemStatSeq::markElements(AdaptInfo *adaptInfo) 
774
  {
775
    FUNCNAME_DBG("ProblemStat::markElements()");
776

777
    TEST_EXIT_DBG(static_cast<unsigned int>(nComponents) == marker.size())
Thomas Witkowski's avatar
Thomas Witkowski committed
778
779
      ("Wrong number of markers!\n");

780
    Flag markFlag = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
781
    for (int i = 0; i < nComponents; i++)
782
      if (marker[i])
Thomas Witkowski's avatar
Thomas Witkowski committed
783
	markFlag |= marker[i]->markMesh(adaptInfo, componentMeshes[i]);
784
    
785
786
787
    return markFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
788

789
  Flag ProblemStatSeq::refineMesh(AdaptInfo *adaptInfo) 
790
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
791
    int nMeshes = static_cast<int>(meshes.size());
792
    Flag refineFlag = 0;
793
794
795
796
    for (int i = 0; i < nMeshes; i++)
      if (adaptInfo->isRefinementAllowed(i))
	refineFlag |= refinementManager->refineMesh(meshes[i]);

797
798
799
    return refineFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
800

801
  Flag ProblemStatSeq::coarsenMesh(AdaptInfo *adaptInfo) 
802
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
803
    int nMeshes = static_cast<int>(meshes.size());
Peter Gottschling's avatar