ParallelCoarseSpaceSolver.cc 11.3 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


13
#include "parallel/ParallelCoarseSpaceSolver.h"
14
#include "parallel/ParallelDofMapping.h"
15
#include "parallel/MatrixNnzStructure.h"
16
17
18
19
20

namespace AMDiS {

  using namespace std;

21
22
23
24
  ParallelCoarseSpaceSolver::ParallelCoarseSpaceSolver(string name)
    : OEMSolver(name), 
      initFileStr(name),
      interiorMap(NULL),
25
      lastMeshNnz(-1),
26
27
28
29
30
      alwaysCreateNnzStructure(false),
      meshDistributor(NULL),
      subdomainLevel(0),
      rStartInterior(0),
      nGlobalOverallInterior(0)
31
  {
32
33
34
35
36
    Parameters::get("parallel->always create nnz structure", 
		    alwaysCreateNnzStructure);
  }


37
  void ParallelCoarseSpaceSolver::setCoarseSpaceDofMapping(ParallelDofMapping *coarseDofs, 
38
							   int component)
39
  {
40
    FUNCNAME("ParallelCoarseSpaceSolver::setCoarseSpaceDofMapping()");
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56

    TEST_EXIT_DBG(coarseDofs)("Should not happen!\n");

    if (component == -1) {
      // === Set coarse space for all components. ===

      coarseSpaceMap.clear();

      int nComponents = coarseDofs->getNumberOfComponents();
      for (int i = 0; i < nComponents; i++)
	coarseSpaceMap[i] = coarseDofs;
    } else {
      // === Set coarse space for just one component. ===

      coarseSpaceMap[component] = coarseDofs;
    }
Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
57
58
59
60

    if (find(uniqueCoarseMap.begin(), uniqueCoarseMap.end(), coarseDofs) == 
	uniqueCoarseMap.end())
      uniqueCoarseMap.push_back(coarseDofs);
61
62
63
  }


64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
  void ParallelCoarseSpaceSolver::setMeshDistributor(MeshDistributor *md)
  {
    FUNCNAME("ParallelCoarseSpaceSolver::setMeshDistributor()");

    meshDistributor = md;
    mpiCommGlobal = meshDistributor->getMpiComm();
    mpiCommLocal = meshDistributor->getLocalMpiComm();
  }


  void ParallelCoarseSpaceSolver::setMeshDistributor(MeshDistributor *md,
						     MPI::Intracomm mpiComm0,
						     MPI::Intracomm mpiComm1)
  {
    FUNCNAME("ParallelCoarseSpaceSolver::setMeshDistributor()");

    meshDistributor = md;
    mpiCommGlobal = mpiComm0;
    mpiCommLocal = mpiComm1; 
  }


86
  void ParallelCoarseSpaceSolver::prepare()
87
  {
88
    FUNCNAME("ParallelCoarseSpaceSolver:prepare()");
89

90
91
    TEST_EXIT(uniqueCoarseMap.size() <= 2)
      ("Not yet implemented for more than two coarse spaces!\n");
92
93
94

    // === Create pointers to PETSc matrix and vector objects. ===

95
96
97
98
99
    int nCoarseMap = uniqueCoarseMap.size();
    mat.resize(nCoarseMap + 1);
    for (int i = 0; i < nCoarseMap + 1; i++)
      mat[i].resize(nCoarseMap + 1);

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
100
101
    vecSol.resize(nCoarseMap + 1);
    vecRhs.resize(nCoarseMap + 1);
102
   
103
104
    // === Create map from component number to its coarse space map. ===

105
    componentIthCoarseMap.resize(coarseSpaceMap.size());
Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
106
    for (unsigned int i = 0; i < coarseSpaceMap.size(); i++) {
107
108
109
110
111
112
113
114
115
116
117
      bool found = false;
      for (int j = 0; j < nCoarseMap; j++) {
	if (coarseSpaceMap[i] == uniqueCoarseMap[j]) {
	  componentIthCoarseMap[i] = j;
	  found = true;
	  break;
	}
      }
      
      TEST_EXIT_DBG(found)("Should not happen!\n");
    }
118
119
120
  }


121
  void ParallelCoarseSpaceSolver::createMatVec(Matrix<DOFMatrix*>& seqMat)
122
  {
123
    FUNCNAME("ParallelCoarseSpaceSolver::createMatVec()");
124

125
126
127
128
129
130
131
132
    // === Prepare coarse space information and generate the correct number ===
    // === of empty PETSc matrix and vector objects.                        ===

    prepare();

    
    // === Update subdomain data (mostly required for multi-level methods) ===

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
133
    updateSubdomainData();
134

135

136
137
138
    // === If required, recompute non zero structure of the matrix. ===

    bool localMatrix = (coarseSpaceMap.size() && subdomainLevel == 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
139

140
    if (checkMeshChange()) {
141
      int nMat = uniqueCoarseMap.size() + 1;
142
143
      nnz.resize(nMat);
      for (int i = 0; i < nMat; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
144
	nnz[i].resize(nMat);
145
146
147
148
149
150
151
152
153
154
155
156
157
158
	for (int j = 0; j < nMat; j++)
	  nnz[i][j].clear();
      }

      nnz[0][0].create(seqMat, mpiCommGlobal, *interiorMap,
		       (coarseSpaceMap.size() == 0 ? &(meshDistributor->getPeriodicMap()) : NULL),
		       meshDistributor->getElementObjectDb(),
		       localMatrix);

      for (int i = 0; i < nMat; i++) {
	for (int j = 0; j < nMat; j++) {
	  if (i == 0 && j == 0)
	    continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
159
160
161
162
163
164
165
	  ParallelDofMapping &rowMap = 
	    (i == 0 ? *interiorMap : *(uniqueCoarseMap[i - 1]));
	  ParallelDofMapping &colMap =
	    (j == 0 ? *interiorMap : *(uniqueCoarseMap[j - 1]));

	  nnz[i][j].create(seqMat, mpiCommGlobal, rowMap, colMap, NULL,
			   meshDistributor->getElementObjectDb());
166
167
168
169
170
	}
      }
    }


Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
171
172
    // === Create PETSc matrices and PETSc vectors with the computed  ===
    // === nnz data structure.                                        ===
173
    
Thomas Witkowski's avatar
Thomas Witkowski committed
174
175
    int nRankInteriorRows = interiorMap->getRankDofs();
    int nOverallInteriorRows = interiorMap->getOverallDofs();
176
177

    MPI::Intracomm mpiGlobal = meshDistributor->getMeshLevelData().getMpiComm(subdomainLevel);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
178
    
179
    if (localMatrix) {
Thomas Witkowski's avatar
Thomas Witkowski committed
180
      MatCreateSeqAIJ(mpiCommLocal, nRankInteriorRows, nRankInteriorRows,
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
181
		      0, nnz[0][0].dnnz,
182
183
		      &mat[0][0]);
    } else {
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
      if (subdomainLevel == 0) {
	MatCreateAIJ(mpiGlobal, nRankInteriorRows, nRankInteriorRows, 
		     nOverallInteriorRows, nOverallInteriorRows,
		     0, nnz[0][0].dnnz, 0, nnz[0][0].onnz,		   
		     &mat[0][0]);    
      } else {
	MSG("WARNING: USE MULTILEVEL, MAKE THIS GENERAL!\n");
	// TODO: MAKE NNZ CALCULATION GENERAL (see also creation of coarse space
	// matrices some lines below)

	MatCreateAIJ(mpiGlobal, nRankInteriorRows, nRankInteriorRows, 
		     nOverallInteriorRows, nOverallInteriorRows,
		     300, PETSC_NULL, 300, PETSC_NULL, &mat[0][0]);    

      }
199
    }
200
    MatSetOption(mat[0][0], MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
201

202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221

    if (subdomainLevel == 0) {
      VecCreateMPI(mpiGlobal, nRankInteriorRows, nOverallInteriorRows, 
		   &vecSol[0]);
      VecCreateMPI(mpiGlobal, nRankInteriorRows, nOverallInteriorRows, 
		   &vecRhs[0]);
    } else {
      MSG("WARNING: USE MULTILEVEL AND THIS IS A VERY BAD HACK, MAKE GENERAL!!\n");

      MPI::Intracomm mpiComm = 
	meshDistributor->getMeshLevelData().getMpiComm(subdomainLevel - 1);
      
      int nAllRows = nRankInteriorRows;
      mpi::globalAdd(nAllRows);
      VecCreateMPI(mpiComm, nRankInteriorRows, nAllRows, 
		   &vecSol[0]);
      VecCreateMPI(mpiComm, nRankInteriorRows, nAllRows, 
		   &vecRhs[0]);
    }

222

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
223
224
225
226
    int nCoarseMap = uniqueCoarseMap.size();
    for (int i = 0; i < nCoarseMap; i++) {
      ParallelDofMapping* cMap = uniqueCoarseMap[i];
      
Thomas Witkowski's avatar
Thomas Witkowski committed
227
228
      int nRankCoarseRows = cMap->getRankDofs();
      int nOverallCoarseRows = cMap->getOverallDofs();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
229
      
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
      if (subdomainLevel == 0) {
	MatCreateAIJ(mpiCommGlobal,
		     nRankCoarseRows, nRankCoarseRows,
		     nOverallCoarseRows, nOverallCoarseRows,
		     0, nnz[i + 1][i + 1].dnnz, 0, nnz[i + 1][i + 1].onnz,
		     &mat[i + 1][i + 1]);
      } else {
	MSG("WARNING: USE MULTILEVEL AND THIS IS A VERY BAD HACK, MAKE GENERAL!!\n");
	
	MPI::Intracomm mpiComm = 
	  meshDistributor->getMeshLevelData().getMpiComm(subdomainLevel - 1);

	MatCreateAIJ(mpiComm,
		     nRankCoarseRows, nRankCoarseRows,
		     nOverallCoarseRows, nOverallCoarseRows,
		     300, PETSC_NULL, 300, PETSC_NULL, 
		     &mat[i + 1][i + 1]);
      }
248
      MatSetOption(mat[i + 1][i + 1], MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
249
250
      cMap->createVec(vecSol[i + 1]);
      cMap->createVec(vecRhs[i + 1]);
Thomas Witkowski's avatar
Thomas Witkowski committed
251
252
    }

253

Thomas Witkowski's avatar
Thomas Witkowski committed
254
    for (int i = 0; i < nCoarseMap + 1; i++) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
255
      for (int j = 0; j < nCoarseMap + 1; j++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
256
257
258
259
260
261
262
263
	if (i == j)
	  continue;

	int nRowsRankMat = (i == 0 ? nRankInteriorRows : uniqueCoarseMap[i - 1]->getRankDofs());
	int nRowsOverallMat = (i == 0 ? nOverallInteriorRows : uniqueCoarseMap[i - 1]->getOverallDofs());

	int nColsRankMat = (j == 0 ? nRankInteriorRows : uniqueCoarseMap[j - 1]->getRankDofs());
	int nColsOverallMat = (j == 0 ? nOverallInteriorRows : uniqueCoarseMap[j - 1]->getOverallDofs());
264

265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
	if (subdomainLevel == 0) {  
	  MatCreateAIJ(mpiCommGlobal,
		       nRowsRankMat, nColsRankMat,
		       nRowsOverallMat, nColsOverallMat,
		       0, nnz[i][j].dnnz, 0, nnz[i][j].onnz,
		       &mat[i][j]);	  
	  MatCreateAIJ(mpiCommGlobal,
		       nColsRankMat, nRowsRankMat,
		       nColsOverallMat, nRowsOverallMat,
		       0, nnz[j][i].dnnz, 0, nnz[j][i].onnz,
		       &mat[j][i]);
	} else {
	  MSG("WARNING: USE MULTILEVEL AND THIS IS A VERY BAD HACK, MAKE GENERAL!!\n");	  

	  if (i == 0) {
	    nRowsOverallMat = nRowsRankMat;
	    mpi::globalAdd(nRowsOverallMat);
	  }

	  if (j == 0) {
	    nColsOverallMat = nColsRankMat;
	    mpi::globalAdd(nColsOverallMat);
	  }

	  MatCreateAIJ(mpiCommGlobal,
		       nRowsRankMat, nColsRankMat,
		       nRowsOverallMat, nColsOverallMat,
		       300, PETSC_NULL, 300, PETSC_NULL,
		       &mat[i][j]);	  
	  MatCreateAIJ(mpiCommGlobal,
		       nColsRankMat, nRowsRankMat,
		       nColsOverallMat, nRowsOverallMat,
		       300, PETSC_NULL, 300, PETSC_NULL,
		       &mat[j][i]);
	}
300
      }
301
    }
302
303
304
  }


305
  void ParallelCoarseSpaceSolver::matDestroy()
306
  {
307
    FUNCNAME("ParallelCoarseSpaceSolver::matDestroy()");
308
309
310
311
312
313
314
315

    int nMatrix = mat.size();
    for (int i = 0; i < nMatrix; i++)
      for (int j = 0; j < nMatrix; j++)
	MatDestroy(&mat[i][j]);
  }


316
  void ParallelCoarseSpaceSolver::vecDestroy()
317
  {
318
    FUNCNAME("ParallelCoarseSpaceSolver::vecDestroy()");
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
319
320
321
322
323
324
325
326
327

    int nVec = vecSol.size();
    for (int i = 0; i < nVec; i++) {
      VecDestroy(&vecSol[i]);
      VecDestroy(&vecRhs[i]);
    }
  }


328
  void ParallelCoarseSpaceSolver::matAssembly()
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
329
  {
330
    FUNCNAME("ParallelCoarseSpaceSolver::matAssembly()");
331
332
333
334
335
336
337
338
339
340

    int nMatrix = mat.size();
    for (int i = 0; i < nMatrix; i++) {
      for (int j = 0; j < nMatrix; j++) {
	MatAssemblyBegin(mat[i][j], MAT_FINAL_ASSEMBLY);
	MatAssemblyEnd(mat[i][j], MAT_FINAL_ASSEMBLY);  
      }
    }
  }

341

342
  void ParallelCoarseSpaceSolver::vecRhsAssembly()
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
343
  {
344
    FUNCNAME("ParallelCoarseSpaceSolver::vecRhsAssembly()");
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
345
346
347
348
349
350
351
352
353

    int nVec = vecRhs.size();
    for (int i = 0; i < nVec; i++) {
      VecAssemblyBegin(vecRhs[i]);
      VecAssemblyEnd(vecRhs[i]);
    }
  }


354
  void ParallelCoarseSpaceSolver::vecSolAssembly()
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
355
  {
356
    FUNCNAME("ParallelCoarseSpaceSolver::vecSolAssembly()");
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
357
358
359
360
361
362
363
364
365

    int nVec = vecRhs.size();
    for (int i = 0; i < nVec; i++) {
      VecAssemblyBegin(vecSol[i]);
      VecAssemblyEnd(vecSol[i]);
    }
  }


366
  bool ParallelCoarseSpaceSolver::checkMeshChange()
367
  {
368
    FUNCNAME("ParallelCoarseSpaceSolver::checkMeshChange()");
369
370
371
372
373
374
375
376
377
378
379
380
381
382

    int recvAllValues = 0;
    int sendValue = 
      static_cast<int>(meshDistributor->getLastMeshChangeIndex() != lastMeshNnz);
    mpiCommGlobal.Allreduce(&sendValue, &recvAllValues, 1, MPI_INT, MPI_SUM);

    if (recvAllValues != 0 || alwaysCreateNnzStructure) {
      lastMeshNnz = meshDistributor->getLastMeshChangeIndex();
      return true;
    }

    return false;
  }

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
383

384
  void ParallelCoarseSpaceSolver::updateSubdomainData()
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
385
  {
386
    FUNCNAME("ParallelCoarseSpaceSolver::updateSubdomainData()");
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406

    if (mpiCommLocal.Get_size() == 1) {
      rStartInterior = 0;
      nGlobalOverallInterior = interiorMap->getOverallDofs();
    } else {
      int groupRowsInterior = 0;
      if (mpiCommLocal.Get_rank() == 0)
	groupRowsInterior = interiorMap->getOverallDofs();

      mpi::getDofNumbering(mpiCommGlobal, groupRowsInterior,
			   rStartInterior, nGlobalOverallInterior);

      int tmp = 0;
      if (mpiCommLocal.Get_rank() == 0)
	tmp = rStartInterior;

      mpiCommLocal.Allreduce(&tmp, &rStartInterior, 1, MPI_INT, MPI_SUM);
    }
  }

407
}