ProblemStat.cc 51.5 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
527
      rhs->setDOFVector(i, new DOFVector<double>(componentSpaces[i], 
						 "rhs[" + lexical_cast<string>(i) + "]"));
528

Thomas Witkowski's avatar
Thomas Witkowski committed
529
      solution->setDOFVector(i, new DOFVector<double>(componentSpaces[i],
530
      						      componentNames[i]));
Thomas Witkowski's avatar
Thomas Witkowski committed
531
532
      solution->getDOFVector(i)->setCoarsenOperation(COARSE_INTERPOL);
      solution->getDOFVector(i)->set(0.0);
533
534
535
    }
  }

536

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

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

572

573
  void ProblemStatSeq::createEstimator()
574
  {
575
    FUNCNAME("ProblemStat::createEstimator()");
576
577

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

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

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


Thomas Witkowski's avatar
Thomas Witkowski committed
600
601
602
603
      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
604
				  rhs->getDOFVector(j));      // TODO: hier eventuell (i) statt (j) ???
605
606
607
    }
  }

608

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

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

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

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

628

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

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

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

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

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

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

659
      if (filename != "")
Thomas Witkowski's avatar
Thomas Witkowski committed
660
	fileWriters.push_back(new FileWriter(numberedName, 
661
					     componentMeshes[i], 
662
					     solution->getDOFVector(i)));      
663
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
664
    
665
    // create a filewrite for groups of components to write vector-valued output
Praetorius, Simon's avatar
Praetorius, Simon committed
666
667
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
    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
					      ));      
	  
	}
      }
    }
697
698

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

704

705
  void ProblemStatSeq::doOtherStuff()
706
  {}
707

708

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

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

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

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

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

734

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

Praetorius, Simon's avatar
Praetorius, Simon committed
739
    Timer t;
740

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

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

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

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

770

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
786

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

795
796
797
    return refineFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
798

799
  Flag ProblemStatSeq::coarsenMesh(AdaptInfo *adaptInfo) 
800
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
801
    int nMeshes = static_cast<int>(meshes.size());
802
    Flag coarsenFlag = 0;
803
804
    for (int i = 0; i < nMeshes; i++)
      if (adaptInfo->isCoarseningAllowed(i))
Thomas Witkowski's avatar
Thomas Witkowski committed