PetscSolverGlobalBlockMatrix.cc 9.11 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
//
// 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.


#include "parallel/PetscSolverGlobalBlockMatrix.h"
#include "parallel/StdMpi.h"
#include "parallel/MpiHelper.h"

namespace AMDiS {

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
19
  void PetscSolverGlobalBlockMatrix::fillPetscMatrix(Matrix<DOFMatrix*> *seqMat)
20
21
  {
    FUNCNAME("PetscSolverGlobalBlockMatrix::fillPetscMatrix()");
22

23
    TEST_EXIT_DBG(meshDistributor)("No mesh distributor object defined!\n");
24
    TEST_EXIT_DBG(interiorMap)("No parallel mapping object defined!\n");
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
25
26
    TEST_EXIT_DBG(seqMat)("No DOF matrix defined!\n");

27
    double wtime = MPI::Wtime();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
28

29
    prepare();
30

31
    const FiniteElemSpace *feSpace = componentSpaces[0];
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
32
    nComponents = seqMat->getNumRows();
33
34
    int nRankRows = (*interiorMap)[feSpace].nRankDofs;
    int nOverallRows = (*interiorMap)[feSpace].nOverallDofs;
35
36
37
38
39

#if (DEBUG != 0)
    MSG("Fill petsc matrix 1 needed %.5f seconds\n", MPI::Wtime() - wtime);
#endif

40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
    if (nBlocks == -1) {
      nBlocks = nComponents;
      for (int i = 0; i < nBlocks; i++)
	componentInBlock[i] = i;
    }

    vector<int> compNthInBlock(nComponents, 0);
    vector<int> blockSize(nBlocks, 0);

    for (int i = 0; i < nComponents; i++) {
      compNthInBlock[i] = blockSize[componentInBlock[i]];
      blockSize[componentInBlock[i]]++;
    }

    nestMat.resize(nBlocks * nBlocks);
55
56
57

    // === Transfer values from DOF matrices to the PETSc matrix. === 

58
59
    for (int i = 0; i < nBlocks; i++)
      for (int j = 0; j < nBlocks; j++)
60
	MatCreateAIJ(domainComm,
61
62
		     nRankRows * blockSize[i], nRankRows * blockSize[j],
		     nOverallRows * blockSize[i], nOverallRows * blockSize[j],
63
64
		     300 * blockSize[i], PETSC_NULL, 
		     300 * blockSize[j], PETSC_NULL,
65
		     &(nestMat[i * nBlocks + j]));
66
			
67
68
    for (int i = 0; i < nComponents; i++)
      for (int j = 0; j < nComponents; j++)
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
69
	if ((*seqMat)[i][j]) {
70
	  int idx = componentInBlock[i] * nBlocks + componentInBlock[j];
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
71
	  setDofMatrix(nestMat[idx], (*seqMat)[i][j], 
72
		       compNthInBlock[i], compNthInBlock[j]);
73
74
	}

75
76
77
78
79
80
81
82
83
84
85
    for (int i = 0; i < nBlocks; i++) {
      for (int j = 0; j < nBlocks; j++) {
	int idx = i * nBlocks + j;
	if (nestMat[idx]) {
	  MatAssemblyBegin(nestMat[idx], MAT_FINAL_ASSEMBLY);
	  MatAssemblyEnd(nestMat[idx], MAT_FINAL_ASSEMBLY);
	}
      }
    }	  
	

86
    MatCreateNest(domainComm, nBlocks, PETSC_NULL, nBlocks, PETSC_NULL,
87
		  &(nestMat[0]), &getMatInterior());
88
89
90
91
92

#if (DEBUG != 0)
    MSG("Fill petsc matrix 2 needed %.5f seconds\n", MPI::Wtime() - wtime);
#endif

93
    matAssembly();
94

95
96
    
    /// initPreconditioner(...)
97
98
99
100
101
102

    // === Init PETSc solver and preconditioner objects. ===

    initSolver(kspInterior);
    KSPGetPC(kspInterior, &pcInterior);
    initPreconditioner(pcInterior);
103

104
105
106
107
108
109
110
111
112
113
114
    MSG("Fill petsc matrix needed %.5f seconds\n", MPI::Wtime() - wtime);
  }


  void PetscSolverGlobalBlockMatrix::fillPetscRhs(SystemVector *vec)
  {
    FUNCNAME("PetscSolverGlobalBlockMatrix::fillPetscRhs()");

    TEST_EXIT_DBG(vec)("NO DOF vector defined!\n");

    nComponents = vec->getSize();
115
    const FiniteElemSpace *feSpace = componentSpaces[0];
116
117
    int nRankRows = (*interiorMap)[feSpace].nRankDofs;
    int nOverallRows = (*interiorMap)[feSpace].nOverallDofs;
118

119
    nestVec.resize(nComponents);
120

121
    for (int i = 0; i < nComponents; i++) {
122
      VecCreateMPI(domainComm, nRankRows, nOverallRows, &(nestVec[i]));
123

124
125
126
127
128
129
      setDofVector(nestVec[i], vec->getDOFVector(i));
      
      VecAssemblyBegin(nestVec[i]);
      VecAssemblyEnd(nestVec[i]);
    }

130
    VecCreateNest(domainComm, nComponents, PETSC_NULL, 
131
		  &(nestVec[0]), &(getVecRhsInterior()));
132

133
    vecRhsAssembly();
134
135
136
  }


137
138
139
140
  void PetscSolverGlobalBlockMatrix::initSolver(KSP &ksp)
  {
    FUNCNAME("PetscSolverGlobalBlockMatrix::initSolver()");

141
    KSPCreate(domainComm, &ksp);
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
    KSPSetOperators(ksp, getMatInterior(), getMatInterior(), 
		    SAME_NONZERO_PATTERN); 
    KSPSetOptionsPrefix(ksp, kspPrefix.c_str());
    KSPSetFromOptions(ksp);
  }


  void PetscSolverGlobalBlockMatrix::exitSolver(KSP ksp)
  {
    FUNCNAME("PetscSolverGlobalBlockMatrix::exitSolver()");

    KSPDestroy(&ksp);
  }


  void PetscSolverGlobalBlockMatrix::initPreconditioner(PC pc)
  {
    FUNCNAME("PetscSolverGlobalBlockMatrix::initPreconditioner()");

    PCSetFromOptions(pc);
  }


  void PetscSolverGlobalBlockMatrix::exitPreconditioner(PC pc)
  {
    FUNCNAME("PetscSolverGlobalBlockMatrix::exitPreconditioner()");
  }


171
172
173
174
175
  void PetscSolverGlobalBlockMatrix::solvePetscMatrix(SystemVector &vec, 
						      AdaptInfo *adaptInfo)
  {
    FUNCNAME("PetscSolverGlobalBlockMatrix::solvePetscMatrix()");

176
    const FiniteElemSpace *feSpace = componentSpaces[0];
177
    VecDuplicate(getVecRhsInterior(), &petscSolVec);
178
179
180
181
182
183
184
185
186
    
    for (int i = 0; i < vec.getSize(); i++)
    {
      Vec tmp;
      VecNestGetSubVec(petscSolVec, i, &tmp);
      setDofVector(tmp, vec.getDOFVector(i));
      VecAssemblyBegin(tmp); 
      VecAssemblyEnd(tmp);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
187

188
    // PETSc.
189
    solve(getVecRhsInterior(), petscSolVec);
190
191
192
193

    // === Transfere values from PETSc's solution vectors to the DOF vectors. ===
    for (int i = 0; i < nComponents; i++) {
      DOFVector<double> &dofvec = *(vec.getDOFVector(i));
194
195

      Vec tmp;
Thomas Witkowski's avatar
Thomas Witkowski committed
196
      VecNestGetSubVec(petscSolVec, i, &tmp);
197

198
      int nRankDofs = (*interiorMap)[feSpace].nRankDofs;
199
200
201
      PetscScalar *vecPointer;
      VecGetArray(tmp, &vecPointer);

202
      DofMap& d = (*interiorMap)[feSpace].getMap();
203
204
205
      for (DofMap::iterator it = d.begin(); it != d.end(); ++it)
	if (it->second.local != -1)
	  dofvec[it->first] = vecPointer[it->second.local];
206

207
208
      VecRestoreArray(tmp, &vecPointer);
    }
209

210
    VecDestroy(&petscSolVec);
211
212
213
214
215
216
217
218
219
220
221

    // === Synchronize DOFs at common DOFs, i.e., DOFs that correspond to ===
    // === more than one partition.                                       ===
    meshDistributor->synchVector(vec);
  }


  void PetscSolverGlobalBlockMatrix::destroyMatrixData()
  {
    FUNCNAME("PetscSolverGlobalBlockMatrix::destroyMatrixData()");

222
223
224
225
    for (unsigned int i = 0; i < nestMat.size(); i++)
      if (nestMat[i] != PETSC_NULL)
	MatDestroy(&(nestMat[i]));

226
    matDestroy();
227

228
229
230
    exitPreconditioner(pcInterior);

    exitSolver(kspInterior);
231
232
233
  }


234
235
236
237
  void PetscSolverGlobalBlockMatrix::destroyVectorData()
  {
    FUNCNAME("PetscSolverGlobalBlockMatrix::destroyVectorData()");

238
239
240
    for (unsigned int i = 0; i < nestVec.size(); i++)
      VecDestroy(&(nestVec[i]));
      
241
    vecDestroy();
242
243
244
  }


245
  void PetscSolverGlobalBlockMatrix::setDofMatrix(Mat& petscMat, 
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
246
						  DOFMatrix* seqMat,
247
248
						  int dispRowBlock, 
						  int dispColBlock)
249
250
251
252
  {
    FUNCNAME("PetscSolverGlobalBlockMatrix::setDofMatrix()");

    TEST_EXIT(petscMat)("No PETSc matrix!\n");
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
253
    TEST_EXIT(seqMat)("No DOFMatrix!\n");
254

255
    const FiniteElemSpace *feSpace = componentSpaces[0];
256

257
258
259
260
    using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
    namespace traits = mtl::traits;
    typedef DOFMatrix::base_matrix_type Matrix;

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
261
262
    traits::col<Matrix>::type col(seqMat->getBaseMatrix());
    traits::const_value<Matrix>::type value(seqMat->getBaseMatrix());
263
264
265
266

    typedef traits::range_generator<row, Matrix>::type cursor_type;
    typedef traits::range_generator<nz, cursor_type>::type icursor_type;

267
268
    int dispRowIndex = (*interiorMap)[feSpace].nRankDofs * dispRowBlock;
    int dispColIndex = (*interiorMap)[feSpace].nRankDofs * dispColBlock;
269

270
271
    vector<int> cols;
    vector<double> values;
272
273
    cols.reserve(3000);
    values.reserve(3000);
274
275
276
277
    
    // === Traverse all rows of the dof matrix and insert row wise the values ===
    // === to the PETSc matrix.                                               ===

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
278
279
    for (cursor_type cursor = begin<row>(seqMat->getBaseMatrix()), 
	   cend = end<row>(seqMat->getBaseMatrix()); cursor != cend; ++cursor) {
280
281

      // Global index of the current row DOF.
282
      int rowIndex = (*interiorMap)[feSpace][*cursor].global + dispRowIndex;
283
284
285
286
287
288
289

      cols.clear();
      values.clear();
      
      for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); 
	   icursor != icend; ++icursor) {	
	// Global index of the current column index.
290
	int colIndex = (*interiorMap)[feSpace][col(*icursor)].global + dispColIndex;
291
292
293
294
295
296
297
298
299
	
	// Ignore all zero entries, expect it is a diagonal entry.
	if (value(*icursor) == 0.0 && rowIndex != colIndex)
	  continue;
	
	// Calculate the exact position of the column index in the PETSc matrix.
	cols.push_back(colIndex);
	values.push_back(value(*icursor));
      }
300
301
302

      MatSetValues(petscMat, 1, &rowIndex, cols.size(), 
 		   &(cols[0]), &(values[0]), ADD_VALUES);	
303
304
305
306
307
    }
  }
  
  
  void PetscSolverGlobalBlockMatrix::setDofVector(Vec& petscVec, 
308
						  DOFVector<double>* vec)
309
310
311
  {
    FUNCNAME("PetscSolverGlobalBlockMatrix::setDofVector()");

312
    const FiniteElemSpace *feSpace = componentSpaces[0];
313

314
315
316
    // Traverse all used DOFs in the dof vector.
    DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
    for (dofIt.reset(); !dofIt.end(); ++dofIt) {
317
      int index = (*interiorMap)[feSpace][dofIt.getDOFIndex()].global;
318
319
320
321
322
323
324
      double value = *dofIt;

      VecSetValues(petscVec, 1, &index, &value, ADD_VALUES);
    }
  }

}