ProblemStat.cc 50.6 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"
Thomas Witkowski's avatar
Thomas Witkowski committed
44
#include "Flag.h"
45
#include "est/Estimator.h"
46
47
#include "io/VtkWriter.h"
#include "io/ValueReader.h"
48
#include "ProblemStatDbg.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
49
#include "Debug.h"
50
51
52

namespace AMDiS {

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
114
115
116

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

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

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

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

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


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

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

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

	// 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
180
	if (adoptProblem->getNumComponents() < nComponents) {
Thomas Witkowski's avatar
Thomas Witkowski committed
181
	  TEST_EXIT(meshes.size() == 1)("Daran muss ich noch arbeiten!\n");
182
	  
Thomas Witkowski's avatar
Thomas Witkowski committed
183
	  componentMeshes.resize(nComponents);
184
	  for (int i = adoptProblem->getNumComponents(); i < nComponents; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
185
	    componentMeshes[i] = componentMeshes[0];
186
	}
187
188
189
      }
    }

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

193
    // === create refinement/corasening-manager ===
194
    if (refinementManager != nullptr && coarseningManager != nullptr) { 
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
      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;
      }
    }

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

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

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

	// 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
230
231
	if (adoptProblem->getNumComponents() < nComponents) {
	  TEST_EXIT(feSpaces.size() == 1)("Daran muss ich noch arbeiten!\n");
232
	  
Thomas Witkowski's avatar
Thomas Witkowski committed
233
	  componentSpaces.resize(nComponents);
Thomas Witkowski's avatar
Thomas Witkowski committed
234
	  for (int i = adoptProblem->getNumComponents(); i < nComponents; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
235
	    componentSpaces[i] = componentSpaces[0];
236
237
	}

238
239
240
      }
    }

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

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

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

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

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

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

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

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

    // === create file writer ===
Thomas Witkowski's avatar
Thomas Witkowski committed
287
    if (initFlag.isSet(INIT_FILEWRITER))
288
289
290
291
292
      createFileWriter();
    
    // === read serialization and init mesh ===
    
    // There are two possiblities where the user can define a serialization
293
294
295
296
    // 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.
297
298

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

    // 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;

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

      // 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) {
316
317
	Parameters::get(name + "->input->serialization filename", 
			serializationFilename);
318
319
	TEST_EXIT(serializationFilename != "")("no serialization file\n");

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

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

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

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

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

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

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

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

	// === do global refinements ===
364

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

    doOtherStuff();
  }

374

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

Thomas Witkowski's avatar
Thomas Witkowski committed
379
    componentMeshes.resize(nComponents);
380
381
382
383
384
385
    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());
386
387

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

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


  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());
418
419

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

437

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

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

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

450
    for (int i = 0; i < nComponents; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
451
      int degree = 1;
452
453
      Parameters::get(name + "->polynomial degree[" + 
		      boost::lexical_cast<string>(i) + "]", degree);
Thomas Witkowski's avatar
Thomas Witkowski committed
454
455
      TEST_EXIT(degree > 0)
	("Poynomial degree in component %d must be larger than zero!\n", i);
456
      TEST_EXIT(componentSpaces[i] == nullptr)("feSpace already created\n");
457

458
      if (feSpaceMap[pair<Mesh*, int>(componentMeshes[i], degree)] == nullptr) {
Thomas Witkowski's avatar
Thomas Witkowski committed
459
460
461
	stringstream s;
	s << name << "->feSpace[" << i << "]";

462
463
	FiniteElemSpace *newFeSpace = 
	  FiniteElemSpace::provideFeSpace(admin, Lagrange::getLagrange(dim, degree),
Thomas Witkowski's avatar
Thomas Witkowski committed
464
					  componentMeshes[i], s.str());
465
	feSpaceMap[pair<Mesh*, int>(componentMeshes[i], degree)] = newFeSpace;
466
	feSpaces.push_back(newFeSpace);
467
      }
468

469
      componentSpaces[i] = feSpaceMap[pair<Mesh*, int>(componentMeshes[i], degree)];
Thomas Witkowski's avatar
Thomas Witkowski committed
470
471
472
    }

    for (int i = 0; i < nComponents; i++) {
473
      for (int j = 0; j < nComponents; j++)
474
	traverseInfo.getMatrix(i, j).setFeSpace(componentSpaces[i], componentSpaces[j]);
475
      
476
      traverseInfo.getVector(i).setFeSpace(componentSpaces[i]);
477
478
479
    }

    // create dof admin for vertex dofs if neccessary
Thomas Witkowski's avatar
Thomas Witkowski committed
480
    for (int i = 0; i < static_cast<int>(meshes.size()); i++) {
481
      if (meshes[i]->getNumberOfDofs(VERTEX) == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
482
	DimVec<int> ln_dof(meshes[i]->getDim(), DEFAULT_VALUE, 0);
483
	ln_dof[VERTEX] = 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
484
	meshes[i]->createDOFAdmin("vertex dofs", ln_dof);      
485
486
487
488
      }
    }
  }

489

490
  void ProblemStatSeq::createMatricesAndVectors()
491
492
493
  {
    // === create vectors and system matrix ===

Thomas Witkowski's avatar
Thomas Witkowski committed
494
    systemMatrix = new Matrix<DOFMatrix*>(nComponents, nComponents);
495
    systemMatrix->set(nullptr);
Thomas Witkowski's avatar
Thomas Witkowski committed
496
497
    rhs = new SystemVector("rhs", componentSpaces, nComponents);
    solution = new SystemVector("solution", componentSpaces, nComponents);
498

Thomas Witkowski's avatar
Thomas Witkowski committed
499
    for (int i = 0; i < nComponents; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
500
      (*systemMatrix)[i][i] = new DOFMatrix(componentSpaces[i], 
Thomas Witkowski's avatar
Thomas Witkowski committed
501
					    componentSpaces[i], "A_ii");
Thomas Witkowski's avatar
Thomas Witkowski committed
502
      (*systemMatrix)[i][i]->setCoupleMatrix(false);
Thomas Witkowski's avatar
Thomas Witkowski committed
503

504
505
      rhs->setDOFVector(i, new DOFVector<double>(componentSpaces[i], 
						 "rhs[" + lexical_cast<string>(i) + "]"));
506

Thomas Witkowski's avatar
Thomas Witkowski committed
507
      solution->setDOFVector(i, new DOFVector<double>(componentSpaces[i],
508
      						      componentNames[i]));
Thomas Witkowski's avatar
Thomas Witkowski committed
509
510
      solution->getDOFVector(i)->setCoarsenOperation(COARSE_INTERPOL);
      solution->getDOFVector(i)->set(0.0);
511
512
513
    }
  }

514

515
  void ProblemStatSeq::createSolver()
516
  {
517
    FUNCNAME("ProblemStat::createSolver()");
518
    
Praetorius, Simon's avatar
Praetorius, Simon committed
519
    // definition of standard-backends
520
521
522
523
#if defined HAVE_PARALLEL_PETSC
    string backend("p_petsc");
#elif defined HAVE_PARALLEL_MTL
    string backend("p_mtl");
524
#elif defined HAVE_SEQ_PETSC
525
526
527
528
529
530
531
532
    string backend("petsc");
#else
    string backend("mtl");
#endif
    
    // === read backend-name ===
    string initFileStr = name + "->solver";    
    Parameters::get(initFileStr + "->backend", backend);
533

534
    // === read solver-name ===
535
    string solverType("0");
536
    Parameters::get(initFileStr, solverType);
537
538
539
540
541
542
543
    
    if (backend != "0" && backend != "no" && backend != "")
      solverType = backend + "_" + solverType;
    
    // === create solver ===
    LinearSolverCreator *solverCreator = 
      dynamic_cast<LinearSolverCreator*>(CreatorMap<LinearSolver>::getCreator(solverType, initFileStr));
544
545
    TEST_EXIT(solverCreator)
      ("No valid solver type found in parameter \"%s\"\n", initFileStr.c_str());
546
    solverCreator->setName(initFileStr);
Thomas Witkowski's avatar
Thomas Witkowski committed
547
    solver = solverCreator->create();
548
549
  }

550

551
  void ProblemStatSeq::createEstimator()
552
  {
553
    FUNCNAME("ProblemStat::createEstimator()");
554
555

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

Thomas Witkowski's avatar
Thomas Witkowski committed
560
    for (int i = 0; i < nComponents; i++) {
561
      TEST_EXIT(estimator[i] == nullptr)("estimator already created\n");
562
563
      string estName = 
	name + "->estimator[" + boost::lexical_cast<string>(i) + "]";
564
565

      // === create estimator ===
566
567
      string estimatorType("0");
      Parameters::get(estName, estimatorType);
568
      EstimatorCreator *estimatorCreator = 
569
	dynamic_cast<EstimatorCreator*>(CreatorMap<Estimator>::getCreator(estimatorType, estName));
Thomas Witkowski's avatar
Thomas Witkowski committed
570
      if (estimatorCreator) {
571
572
	estimatorCreator->setName(estName);
	estimatorCreator->setRow(i);
573
	estimatorCreator->setSolution(solution->getDOFVector(i));
Thomas Witkowski's avatar
Thomas Witkowski committed
574
	estimator[i] = estimatorCreator->create();
575
576
577
      }


Thomas Witkowski's avatar
Thomas Witkowski committed
578
579
580
581
      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
582
				  rhs->getDOFVector(j));      // TODO: hier eventuell (i) statt (j) ???
583
584
585
    }
  }

586

587
  void ProblemStatSeq::createMarker()
588
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
589
    int nMarkersCreated = 0;
590

591
    for (int i = 0; i < nComponents; i++) {
592
      marker[i] = Marker::createMarker
593
	(name + "->marker[" + boost::lexical_cast<string>(i) + "]", i);
594

595
      if (marker[i]) {
Thomas Witkowski's avatar
Thomas Witkowski committed
596
597
598
599
	nMarkersCreated++;

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

606

607
  void ProblemStatSeq::createFileWriter()
608
  {
609
    FUNCNAME("ProblemStat::createFileWriter()");
610
611
  
    // Create one filewriter for all components of the problem
612
613
614
    string numberedName  = name + "->output";
    string filename = "";
    Parameters::get(numberedName + "->filename", filename);
615
616

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
623
	solutionList[i] = solution->getDOFVector(i);
624
625
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
626
      fileWriters.push_back(new FileWriter(numberedName,
627
628
					   componentMeshes[0],
					   solutionList));
629
630
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
631
    // Create own filewriters for each component of the problem
632
    for (int i = 0; i < nComponents; i++) {
633
      numberedName = name + "->output[" + boost::lexical_cast<string>(i) + "]";
634
      filename = "";
635
      Parameters::get(numberedName + "->filename", filename);
636

637
      if (filename != "")
Thomas Witkowski's avatar
Thomas Witkowski committed
638
	fileWriters.push_back(new FileWriter(numberedName, 
639
					     componentMeshes[i], 
640
					     solution->getDOFVector(i)));      
641
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
    
    // create a filewrite for groups of components to write vector-values output
    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
					      ));      
	  
	}
      }
    }
675
676

    int writeSerialization = 0;
677
    Parameters::get(name + "->output->write serialization", writeSerialization);
678
    if (writeSerialization)
679
      fileWriters.push_back(new Serializer<ProblemStatSeq>(this));
680
681
  }

682

683
  void ProblemStatSeq::doOtherStuff()
684
  {}
685

686

687
688
689
  void ProblemStatSeq::solve(AdaptInfo *adaptInfo, 	       
			     bool createMatrixData,
			     bool storeMatrixData)
690
691
692
  {
    FUNCNAME("Problem::solve()");

Thomas Witkowski's avatar
Thomas Witkowski committed
693
    if (!solver) {
694
695
696
697
      WARNING("no solver\n");
      return;
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
698
    Timer t;
699
700
    solver->solveSystem(solverMatrix, *solution, *rhs, 
			createMatrixData, storeMatrixData);
701

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

706
    adaptInfo->setSolverIterations(solver->getIterations());
Thomas Witkowski's avatar
Thomas Witkowski committed
707
708
709
    adaptInfo->setMaxSolverIterations(solver->getMaxIterations());
    adaptInfo->setSolverTolerance(solver->getTolerance());
    adaptInfo->setSolverResidual(solver->getResidual());
710
711
  }

712

713
  void ProblemStatSeq::estimate(AdaptInfo *adaptInfo) 
714
  {
715
    FUNCNAME("ProblemStat::estimate()");
716

Praetorius, Simon's avatar
Praetorius, Simon committed
717
    Timer t;
718

719
    if (computeExactError) {
Thomas Witkowski's avatar
Thomas Witkowski committed
720
      computeError(adaptInfo);
721
722
    } else {
      for (int i = 0; i < nComponents; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
723
	Estimator *scalEstimator = estimator[i];
724
725
	
	if (scalEstimator) {
726
727
	  traverseInfo.updateStatus();
	  scalEstimator->setTraverseInfo(traverseInfo);
728
	  scalEstimator->estimate(adaptInfo->getTimestep());
729

730
731
	  adaptInfo->setEstSum(scalEstimator->getErrorSum(), i);
	  adaptInfo->setEstMax(scalEstimator->getErrorMax(), i);
732
733
734
735
736

	  if (adaptInfo->getRosenbrockMode() == false) {
	    adaptInfo->setTimeEstSum(scalEstimator->getTimeEst(), i);
	    adaptInfo->setTimeEstMax(scalEstimator->getTimeEstMax(), i);
	  }
737
	}
738
739
740
      }
    }

741
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
742
    MPI::COMM_WORLD.Barrier();
743
#endif
Praetorius, Simon's avatar
Praetorius, Simon committed
744
745
    INFO(info, 8)("estimation of the error needed %.5f seconds\n", 
		  t.elapsed());
746
747
  }

748

749
  Flag ProblemStatSeq::markElements(AdaptInfo *adaptInfo) 
750
  {
751
    FUNCNAME_DBG("ProblemStat::markElements()");
752

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

756
    Flag markFlag = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
757
    for (int i = 0; i < nComponents; i++)
758
      if (marker[i])
Thomas Witkowski's avatar
Thomas Witkowski committed
759
	markFlag |= marker[i]->markMesh(adaptInfo, componentMeshes[i]);
760
    
761
762
763
    return markFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
764

765
  Flag ProblemStatSeq::refineMesh(AdaptInfo *adaptInfo) 
766
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
767
    int nMeshes = static_cast<int>(meshes.size());
768
    Flag refineFlag = 0;
769
770
771
772
    for (int i = 0; i < nMeshes; i++)
      if (adaptInfo->isRefinementAllowed(i))
	refineFlag |= refinementManager->refineMesh(meshes[i]);

773
774
775
    return refineFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
776

777
  Flag ProblemStatSeq::coarsenMesh(AdaptInfo *adaptInfo) 
778
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
779
    int nMeshes = static_cast<int>(meshes.size());
780
    Flag coarsenFlag = 0;
781
782
    for (int i = 0; i < nMeshes; i++)
      if (adaptInfo->isCoarseningAllowed(i))
Thomas Witkowski's avatar
Thomas Witkowski committed
783
	coarsenFlag |= coarseningManager->coarsenMesh(meshes[i]);
784
785
786
787

    return coarsenFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
788