PetscSolverGlobalBlockMatrix.cc 8.32 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 30 31
    petscData.init(interiorMap, coarseSpaceMap, 
		   subdomainLevel, mpiCommLocal, mpiCommGlobal,
		   meshDistributor);
32

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

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

42 43 44 45 46 47 48 49 50 51 52 53 54 55 56
    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);
57 58 59

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

60 61
    for (int i = 0; i < nBlocks; i++)
      for (int j = 0; j < nBlocks; j++)
62 63 64 65 66 67
	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]));
68
			
69 70
    for (int i = 0; i < nComponents; i++)
      for (int j = 0; j < nComponents; j++)
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
71
	if ((*seqMat)[i][j]) {
72
	  int idx = componentInBlock[i] * nBlocks + componentInBlock[j];
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
73
	  setDofMatrix(nestMat[idx], (*seqMat)[i][j], 
74
		       compNthInBlock[i], compNthInBlock[j]);
75 76
	}

77 78 79 80 81 82 83 84 85 86 87
    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);
	}
      }
    }	  
	

88
    MatCreateNest(mpiCommGlobal, 
89
		  nBlocks, PETSC_NULL, nBlocks, PETSC_NULL, 
90
		  &(nestMat[0]), &petscData.getInteriorMat());
91 92 93 94 95

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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
96
    petscData.matAssembly();
97 98

    // === Init PETSc solver. ===
99
    KSPCreate(mpiCommGlobal, &kspInterior);
100 101 102
    KSPSetOperators(kspInterior, 
		    petscData.getInteriorMat(), 
		    petscData.getInteriorMat(), SAME_NONZERO_PATTERN); 
103
    KSPSetFromOptions(kspInterior);
104

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

120
    nestVec.resize(nComponents);
121

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

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

131
    VecCreateNest(mpiCommGlobal, nComponents, PETSC_NULL, 
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
132
		  &(nestVec[0]), &(getRhsInterior()));
133

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
134
    petscData.vecRhsAssembly();
135 136 137 138 139 140 141 142
  }


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

143 144
    KSPGetPC(kspInterior, &pcInterior);
    setBlockPreconditioner(pcInterior);
145

146
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
147
    VecDuplicate(getRhsInterior(), &petscSolVec);
Thomas Witkowski's avatar
Thomas Witkowski committed
148

149
    // PETSc.
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
150
    solve(getRhsInterior(), petscSolVec);
151 152 153 154

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

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

159
      int nRankDofs = (*interiorMap)[feSpace].nRankDofs;
160 161 162
      PetscScalar *vecPointer;
      VecGetArray(tmp, &vecPointer);

163
      DofMap& d = (*interiorMap)[feSpace].getMap();
164 165 166
      for (DofMap::iterator it = d.begin(); it != d.end(); ++it)
	if (it->second.local != -1)
	  dofvec[it->first] = vecPointer[it->second.local];
167

168 169
      VecRestoreArray(tmp, &vecPointer);
    }
170 171 172 173 174 175 176 177 178 179 180 181


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


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

182 183 184 185
    for (unsigned int i = 0; i < nestMat.size(); i++)
      if (nestMat[i] != PETSC_NULL)
	MatDestroy(&(nestMat[i]));

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
186
    petscData.matDestroy();
187

188
    KSPDestroy(&kspInterior);
189 190 191
  }


192 193 194 195
  void PetscSolverGlobalBlockMatrix::destroyVectorData()
  {
    FUNCNAME("PetscSolverGlobalBlockMatrix::destroyVectorData()");

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
196
    petscData.vecDestroy();
197 198 199 200 201

    VecDestroy(&petscSolVec);
  }


202
  void PetscSolverGlobalBlockMatrix::setDofMatrix(Mat& petscMat, 
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
203
						  DOFMatrix* seqMat,
204 205
						  int dispRowBlock, 
						  int dispColBlock)
206 207 208 209
  {
    FUNCNAME("PetscSolverGlobalBlockMatrix::setDofMatrix()");

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

212 213
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

214 215 216 217
    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
218 219
    traits::col<Matrix>::type col(seqMat->getBaseMatrix());
    traits::const_value<Matrix>::type value(seqMat->getBaseMatrix());
220 221 222 223

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

224 225
    int dispRowIndex = (*interiorMap)[feSpace].nRankDofs * dispRowBlock;
    int dispColIndex = (*interiorMap)[feSpace].nRankDofs * dispColBlock;
226

227 228 229 230 231 232 233 234
    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
235 236
    for (cursor_type cursor = begin<row>(seqMat->getBaseMatrix()), 
	   cend = end<row>(seqMat->getBaseMatrix()); cursor != cend; ++cursor) {
237 238

      // Global index of the current row DOF.
239
      int rowIndex = (*interiorMap)[feSpace][*cursor].global + dispRowIndex;
240 241 242 243 244 245 246

      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.
247
	int colIndex = (*interiorMap)[feSpace][col(*icursor)].global + dispColIndex;
248 249 250 251 252 253 254 255 256
	
	// 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));
      }
257 258 259

      MatSetValues(petscMat, 1, &rowIndex, cols.size(), 
 		   &(cols[0]), &(values[0]), ADD_VALUES);	
260 261 262 263 264
    }
  }
  
  
  void PetscSolverGlobalBlockMatrix::setDofVector(Vec& petscVec, 
265
						  DOFVector<double>* vec)
266 267 268
  {
    FUNCNAME("PetscSolverGlobalBlockMatrix::setDofVector()");

269 270
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

271 272 273
    // Traverse all used DOFs in the dof vector.
    DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
    for (dofIt.reset(); !dofIt.end(); ++dofIt) {
274
      int index = (*interiorMap)[feSpace][dofIt.getDOFIndex()].global;
275 276 277 278 279 280 281
      double value = *dofIt;

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

}