ParallelDomainVec.cc 4.86 KB
Newer Older
1
#include "boost/lexical_cast.hpp"
2
3
4
#include "ParallelDomainVec.h"
#include "ProblemVec.h"
#include "ProblemInstat.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
5
#include "SystemVector.h"
6
#include "Serializer.h"
7
#include "BoundaryManager.h"
8
9
10

namespace AMDiS {

11
12
  using boost::lexical_cast;

13
  ParallelDomainVec::ParallelDomainVec(ProblemVec *problem,
14
				       ProblemInstatVec *problemInstat)
15
    : ParallelDomainBase(problem, 
16
17
18
19
20
21
22
23
24
25
			 problemInstat, 
			 problem->getFESpace(0),
			 problem->getRefinementManager(0)),
      probVec(problem)
  {
    FUNCNAME("ParallelDomainVec::ParallelDomainVec()");

    info = probVec->getInfo();
    nComponents = probVec->getNumComponents();

26
27
    // Test if all fe spaces are equal. Yet, different fe spaces are not supported for
    // domain parallelization.
28
29
30
31
    const FiniteElemSpace *fe = probVec->getFESpace(0);
    for (int i = 0; i < nComponents; i++)
      TEST_EXIT(fe == probVec->getFESpace(i))
	("Parallelization does not supported different FE spaces!\n");
32
33
34
35
36
37

    // Create parallel serialization file writer, if needed.
    int writeSerialization = 0;
    GET_PARAMETER(0, name + "->output->write serialization", "%d", &writeSerialization);
    if (writeSerialization)
      problem->getFileWriterList().push_back(new Serializer<ParallelDomainVec>(this));
38
39
40
41
42
43
44
45
46
47
48
49

    int readSerialization = 0;
    GET_PARAMETER(0, name + "->input->read serialization", "%d", &readSerialization);
    if (readSerialization) {
      std::string filename = "";
      GET_PARAMETER(0, name + "->input->serialization filename", &filename);
      filename += ".p" + lexical_cast<std::string>(mpiRank);
      MSG("Start serialization with %s\n", filename.c_str());
      std::ifstream in(filename.c_str());
      deserialize(in);
      in.close();
    }
50
51
52
53
54
55
56
57
58
  }


  void ParallelDomainVec::initParallelization(AdaptInfo *adaptInfo)
  {
    FUNCNAME("ParallelDomainVec::initParallelization()");

    ParallelDomainBase::initParallelization(adaptInfo);

Thomas Witkowski's avatar
Thomas Witkowski committed
59
    for (int i = 0; i < nComponents; i++) {
60
      for (int j = 0; j < nComponents; j++)
61
 	if (probVec->getSystemMatrix(i, j))
62
 	  probVec->getSystemMatrix(i, j)->setRankDofs(isRankDof);      
63
64
65

      TEST_EXIT_DBG(probVec->getRHS()->getDOFVector(i))("No rhs vector!\n");
      TEST_EXIT_DBG(probVec->getSolution()->getDOFVector(i))("No solution vector!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
66
67
68
69

      probVec->getRHS()->getDOFVector(i)->setRankDofs(isRankDof);
      probVec->getSolution()->getDOFVector(i)->setRankDofs(isRankDof);
    }
70

71
72
73
    // === Remove periodic boundary conditions in sequential problem definition. ===

    // Remove periodic boundaries in boundary manager on matrices and vectors.
74
    for (int i = 0; i < nComponents; i++) {
75
76
77
78
79
      for (int j = 0; j < nComponents; j++) {
	DOFMatrix* mat = probVec->getSystemMatrix(i, j);
 	if (mat && mat->getBoundaryManager())
	  removeBoundaryCondition(const_cast<BoundaryIndexMap&>(mat->getBoundaryManager()->getBoundaryConditionMap()));
      }
80
81
82

      if (probVec->getSolution()->getDOFVector(i)->getBoundaryManager())
	removeBoundaryCondition(const_cast<BoundaryIndexMap&>(probVec->getSolution()->getDOFVector(i)->getBoundaryManager()->getBoundaryConditionMap()));
83
84
85

      if (probVec->getRHS()->getDOFVector(i)->getBoundaryManager())
	removeBoundaryCondition(const_cast<BoundaryIndexMap&>(probVec->getRHS()->getDOFVector(i)->getBoundaryManager()->getBoundaryConditionMap()));
86
87
88
89

      std::cout << "TEST SOL " << probVec->getSolution()->getDOFVector(i)  << ": " << probVec->getSolution()->getDOFVector(i)->getCoarsenOperation() << std::endl;

      std::cout << "TEST RHS " << probVec->getRHS()->getDOFVector(i) << ": " << probVec->getRHS()->getDOFVector(i)->getCoarsenOperation() << std::endl;
90
91
    }

92
    // Remove periodic boundaries on elements in mesh.
93
94
95
96
97
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh,  -1, Mesh::CALL_EVERY_EL_PREORDER);
    while (elInfo) {
      elInfo->getElement()->deleteElementData(PERIODIC);
      elInfo = stack.traverseNext(elInfo);
98
    }    
99
100
  }

101
102
103
104
105
  void ParallelDomainVec::exitParallelization(AdaptInfo *adaptInfo)
  {
    ParallelDomainBase::exitParallelization(adaptInfo);
  }

106

107
108
109
110
111
112
113
114
115
116
  void ParallelDomainVec::solve()
  {
    FUNCNAME("ParallelDomainVec::solve()");

#ifdef _OPENMP
    double wtime = omp_get_wtime();
#endif
    clock_t first = clock();

    fillPetscMatrix(probVec->getSystemMatrix(), probVec->getRHS());      
117
    solvePetscMatrix(*(probVec->getSolution()));
118
119
120
121
122
123
124
125

#ifdef _OPENMP
    INFO(info, 8)("solution of discrete system needed %.5f seconds system time / %.5f seconds wallclock time\n",
		   TIME_USED(first, clock()),
		   omp_get_wtime() - wtime);
#else
    INFO(info, 8)("solution of discrete system needed %.5f seconds\n",
		   TIME_USED(first, clock()));
126
#endif    
127
128
  }

129
130
131
132
133
134
135
136
137
138
139
140

  void ParallelDomainVec::removeBoundaryCondition(BoundaryIndexMap& boundaryMap)
  {
    BoundaryIndexMap::iterator it = boundaryMap.begin();
    while (it != boundaryMap.end()) {
      if (it->second->isPeriodic())
	boundaryMap.erase(it++);
      else
	++it;      
    }    
  }

141
}