PetscSolverGlobalBlockMatrix.cc 7.19 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
//
// 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 {

  void PetscSolverGlobalBlockMatrix::fillPetscMatrix(Matrix<DOFMatrix*> *mat)
  {
    FUNCNAME("PetscSolverGlobalBlockMatrix::fillPetscMatrix()");

    TEST_EXIT_DBG(meshDistributor)("No mesh distributor object defined!\n");
    TEST_EXIT_DBG(mat)("No DOF matrix defined!\n");

    double wtime = MPI::Wtime();
27
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
28
    nComponents = mat->getNumRows();
29
30
    int nRankRows = meshDistributor->getNumberRankDofs(feSpace);
    int nOverallRows = meshDistributor->getNumberOverallDofs(feSpace);
31
32
33
34
35

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

36
    nestMat.resize(nComponents * nComponents);
37
38
39
40
41
42
43
44
45
46

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

    for (int i = 0; i < nComponents; i++)
      for (int j = 0; j < nComponents; j++)
	if ((*mat)[i][j]) {
	  MatCreateMPIAIJ(PETSC_COMM_WORLD,
			  nRankRows, nRankRows,
			  nOverallRows, nOverallRows,
			  30, PETSC_NULL, 30, PETSC_NULL,
47
48
49
50
			  &(nestMat[i * nComponents + j]));
	  setDofMatrix(nestMat[i * nComponents + j], (*mat)[i][j]);
	  MatAssemblyBegin(nestMat[i * nComponents + j], MAT_FINAL_ASSEMBLY);
	  MatAssemblyEnd(nestMat[i * nComponents + j], MAT_FINAL_ASSEMBLY);
51
	} else {
52
	  nestMat[i * nComponents + j] = PETSC_NULL;
53
54
55
56
	}

    MatCreateNest(PETSC_COMM_WORLD, 
		  nComponents, PETSC_NULL, nComponents, PETSC_NULL, 
57
		  &(nestMat[0]), &petscMatrix);
58
59
60
61
62
63
64
65
66
67
68
69

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

    MatAssemblyBegin(petscMatrix, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(petscMatrix, MAT_FINAL_ASSEMBLY);

    // === Init PETSc solver. ===
    KSPCreate(PETSC_COMM_WORLD, &solver);
    KSPSetOperators(solver, petscMatrix, petscMatrix, SAME_NONZERO_PATTERN); 
    KSPSetFromOptions(solver);
70
71
72

    KSPGetPC(solver, &pc);
    setBlockPreconditioner(pc);
73
74
75
76
77
78
79
80
81
82
83
84

    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();
85
86
87
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
    int nRankRows = meshDistributor->getNumberRankDofs(feSpace);
    int nOverallRows = meshDistributor->getNumberOverallDofs(feSpace);
88

89
    nestVec.resize(nComponents);
90

91
92
    for (int i = 0; i < nComponents; i++) {
      VecCreateMPI(PETSC_COMM_WORLD, nRankRows, nOverallRows, &(nestVec[i]));
93

94
95
96
97
98
99
100
101
      setDofVector(nestVec[i], vec->getDOFVector(i));
      
      VecAssemblyBegin(nestVec[i]);
      VecAssemblyEnd(nestVec[i]);
    }

    VecCreateNest(PETSC_COMM_WORLD, nComponents, PETSC_NULL, 
		  &(nestVec[0]), &petscRhsVec);
102
103
104
105
106
107
108
109
110
111
112

    VecAssemblyBegin(petscRhsVec);
    VecAssemblyEnd(petscRhsVec);
  }


  void PetscSolverGlobalBlockMatrix::solvePetscMatrix(SystemVector &vec, 
						      AdaptInfo *adaptInfo)
  {
    FUNCNAME("PetscSolverGlobalBlockMatrix::solvePetscMatrix()");

113
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
Thomas Witkowski's avatar
Thomas Witkowski committed
114
115
    VecDuplicate(petscRhsVec, &petscSolVec);

116
    // PETSc.
Thomas Witkowski's avatar
Thomas Witkowski committed
117
    KSPSolve(solver, petscRhsVec, petscSolVec);
118
119
120
121

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

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

126
      int nRankDofs = meshDistributor->getNumberRankDofs(feSpace);
127
128
129
      PetscScalar *vecPointer;
      VecGetArray(tmp, &vecPointer);

130
      for (int j = 0; j < nRankDofs; j++)
131
	dofvec[meshDistributor->mapLocalToDof(feSpace, j)] = vecPointer[j]; 
132

133
134
      VecRestoreArray(tmp, &vecPointer);
    }
135
136
137
138
139
140
141
142
143


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


    // === Destroy PETSc's variables. ===
    VecDestroy(&petscRhsVec);
144
145
    for (int i = 0; i < nComponents; i++)
      VecDestroy(&(nestVec[i]));
Thomas Witkowski's avatar
Thomas Witkowski committed
146
147

    VecDestroy(&petscSolVec);
148
149
150
151
152
153
154
  }


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

155
156
157
158
    for (unsigned int i = 0; i < nestMat.size(); i++)
      if (nestMat[i] != PETSC_NULL)
	MatDestroy(&(nestMat[i]));

159
160
161
162
163
164
165
166
167
168
169
170
    MatDestroy(&petscMatrix);
    KSPDestroy(&solver);
  }


  void PetscSolverGlobalBlockMatrix::setDofMatrix(Mat& petscMat, DOFMatrix* mat)
  {
    FUNCNAME("PetscSolverGlobalBlockMatrix::setDofMatrix()");

    TEST_EXIT(mat)("No DOFMatrix!\n");
    TEST_EXIT(petscMat)("No PETSc matrix!\n");

171
172
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
    using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
    namespace traits = mtl::traits;
    typedef DOFMatrix::base_matrix_type Matrix;

    traits::col<Matrix>::type col(mat->getBaseMatrix());
    traits::const_value<Matrix>::type value(mat->getBaseMatrix());

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

    vector<int> cols;
    vector<double> values;
    cols.reserve(300);
    values.reserve(300);
    
    // === Traverse all rows of the dof matrix and insert row wise the values ===
    // === to the PETSc matrix.                                               ===

    for (cursor_type cursor = begin<row>(mat->getBaseMatrix()), 
	   cend = end<row>(mat->getBaseMatrix()); cursor != cend; ++cursor) {

      // Global index of the current row DOF.
195
      int rowIndex = meshDistributor->mapDofToGlobal(feSpace, *cursor);
196
197
198
199
200
201
202

      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.
203
	int colIndex = meshDistributor->mapDofToGlobal(feSpace, col(*icursor));
204
205
206
207
208
209
210
211
212
	
	// 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));
      }
213
214
215

      MatSetValues(petscMat, 1, &rowIndex, cols.size(), 
 		   &(cols[0]), &(values[0]), ADD_VALUES);	
216
217
218
219
220
    }
  }
  
  
  void PetscSolverGlobalBlockMatrix::setDofVector(Vec& petscVec, 
221
						  DOFVector<double>* vec)
222
223
224
  {
    FUNCNAME("PetscSolverGlobalBlockMatrix::setDofVector()");

225
226
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

227
228
229
    // Traverse all used DOFs in the dof vector.
    DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
    for (dofIt.reset(); !dofIt.end(); ++dofIt) {
230
      int index = meshDistributor->mapDofToGlobal(feSpace, dofIt.getDOFIndex());
231
232
233
234
235
236
237
      double value = *dofIt;

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

}