PetscSolverGlobalBlockMatrix.cc 8.78 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
    createMatVec(*seqMat);
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 61 62 63 64 65
	MatCreateAIJ(mpiCommGlobal,
		     nRankRows * blockSize[i], nRankRows * blockSize[j],
		     nOverallRows * blockSize[i], nOverallRows * blockSize[j],
		     30 * blockSize[i], PETSC_NULL, 
		     30 * blockSize[j], PETSC_NULL,
		     &(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 87
    MatCreateNest(mpiCommGlobal, nBlocks, PETSC_NULL, nBlocks, PETSC_NULL,
		  &(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 97 98 99 100

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

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

102 103 104 105 106 107 108 109 110 111 112
    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();
113
    const FiniteElemSpace *feSpace = componentSpaces[0];
114 115
    int nRankRows = (*interiorMap)[feSpace].nRankDofs;
    int nOverallRows = (*interiorMap)[feSpace].nOverallDofs;
116

117
    nestVec.resize(nComponents);
118

119
    for (int i = 0; i < nComponents; i++) {
120
      VecCreateMPI(mpiCommGlobal, nRankRows, nOverallRows, &(nestVec[i]));
121

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

128
    VecCreateNest(mpiCommGlobal, nComponents, PETSC_NULL, 
129
		  &(nestVec[0]), &(getVecRhsInterior()));
130

131
    vecRhsAssembly();
132 133 134
  }


135 136 137 138 139 140 141 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
  void PetscSolverGlobalBlockMatrix::initSolver(KSP &ksp)
  {
    FUNCNAME("PetscSolverGlobalBlockMatrix::initSolver()");

    KSPCreate(mpiCommGlobal, &ksp);
    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()");
  }


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

174
    const FiniteElemSpace *feSpace = componentSpaces[0];
175
    VecDuplicate(getVecRhsInterior(), &petscSolVec);
Thomas Witkowski's avatar
Thomas Witkowski committed
176

177
    // PETSc.
178
    solve(getVecRhsInterior(), petscSolVec);
179 180 181 182

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

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

187
      int nRankDofs = (*interiorMap)[feSpace].nRankDofs;
188 189 190
      PetscScalar *vecPointer;
      VecGetArray(tmp, &vecPointer);

191
      DofMap& d = (*interiorMap)[feSpace].getMap();
192 193 194
      for (DofMap::iterator it = d.begin(); it != d.end(); ++it)
	if (it->second.local != -1)
	  dofvec[it->first] = vecPointer[it->second.local];
195

196 197
      VecRestoreArray(tmp, &vecPointer);
    }
198 199 200 201 202 203 204 205 206 207 208 209


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


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

210 211 212 213
    for (unsigned int i = 0; i < nestMat.size(); i++)
      if (nestMat[i] != PETSC_NULL)
	MatDestroy(&(nestMat[i]));

214
    matDestroy();
215

216 217 218
    exitPreconditioner(pcInterior);

    exitSolver(kspInterior);
219 220 221
  }


222 223 224 225
  void PetscSolverGlobalBlockMatrix::destroyVectorData()
  {
    FUNCNAME("PetscSolverGlobalBlockMatrix::destroyVectorData()");

226
    vecDestroy();
227 228 229 230 231

    VecDestroy(&petscSolVec);
  }


232
  void PetscSolverGlobalBlockMatrix::setDofMatrix(Mat& petscMat, 
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
233
						  DOFMatrix* seqMat,
234 235
						  int dispRowBlock, 
						  int dispColBlock)
236 237 238 239
  {
    FUNCNAME("PetscSolverGlobalBlockMatrix::setDofMatrix()");

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

242
    const FiniteElemSpace *feSpace = componentSpaces[0];
243

244 245 246 247
    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
248 249
    traits::col<Matrix>::type col(seqMat->getBaseMatrix());
    traits::const_value<Matrix>::type value(seqMat->getBaseMatrix());
250 251 252 253

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

254 255
    int dispRowIndex = (*interiorMap)[feSpace].nRankDofs * dispRowBlock;
    int dispColIndex = (*interiorMap)[feSpace].nRankDofs * dispColBlock;
256

257 258 259 260 261 262 263 264
    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.                                               ===

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
265 266
    for (cursor_type cursor = begin<row>(seqMat->getBaseMatrix()), 
	   cend = end<row>(seqMat->getBaseMatrix()); cursor != cend; ++cursor) {
267 268

      // Global index of the current row DOF.
269
      int rowIndex = (*interiorMap)[feSpace][*cursor].global + dispRowIndex;
270 271 272 273 274 275 276

      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.
277
	int colIndex = (*interiorMap)[feSpace][col(*icursor)].global + dispColIndex;
278 279 280 281 282 283 284 285 286
	
	// 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));
      }
287 288 289

      MatSetValues(petscMat, 1, &rowIndex, cols.size(), 
 		   &(cols[0]), &(values[0]), ADD_VALUES);	
290 291 292 293 294
    }
  }
  
  
  void PetscSolverGlobalBlockMatrix::setDofVector(Vec& petscVec, 
295
						  DOFVector<double>* vec)
296 297 298
  {
    FUNCNAME("PetscSolverGlobalBlockMatrix::setDofVector()");

299
    const FiniteElemSpace *feSpace = componentSpaces[0];
300

301 302 303
    // Traverse all used DOFs in the dof vector.
    DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
    for (dofIt.reset(); !dofIt.end(); ++dofIt) {
304
      int index = (*interiorMap)[feSpace][dofIt.getDOFIndex()].global;
305 306 307 308 309 310 311
      double value = *dofIt;

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

}