PetscSolverGlobalMatrix.cc 31 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: 
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
 * 
 ******************************************************************************/
Thomas Witkowski's avatar
Thomas Witkowski committed
20

21
#include <mpi.h>
22
// #include "DirichletBC.h"
23
#include "DOFVector.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
24 25 26
#include "parallel/PetscSolverGlobalMatrix.h"
#include "parallel/StdMpi.h"
#include "parallel/MpiHelper.h"
27
#include "solver/PetscTypes.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
28

29 30
using namespace std;

31
namespace AMDiS { namespace Parallel {
32

33
  PetscSolverGlobalMatrix::PetscSolverGlobalMatrix(string name, bool setOptions)
34 35 36
    : PetscSolver(name),
      zeroStartVector(false),
      printMatInfo(false)
37 38
  { FUNCNAME_DBG("PetscSolverGlobalMatrix()");
  
39
    bool matSolverPackage = false;
40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68
    if (setOptions) {
      PetscParameters params;
      
      // set the solver
      std::string solverName = "petsc";
      Parameters::get(name, solverName);
      if (solverName == "petsc") 
	Parameters::get(name + "->ksp_type", solverName);
      
      std::string kspSolver = params.solverMap[solverName];
      
      if (params.matSolverPackage.find(kspSolver) != params.matSolverPackage.end()) {
	// direct solvers
	PetscOptionsInsertString(("-" + kspPrefix + "ksp_type preonly").c_str());
	PetscOptionsInsertString(("-" + kspPrefix + "pc_type lu").c_str());
	PetscOptionsInsertString(("-" + kspPrefix + "pc_factor_mat_solver_package " + kspSolver).c_str());
	setMaxIterations(1);
	zeroStartVector = true;
	matSolverPackage = true;
      } else if (params.emptyParam.find(kspSolver) == params.emptyParam.end() && solverName != "petsc") {    
	// other solvers
	PetscOptionsInsertString(("-" + kspPrefix + "ksp_type " + kspSolver).c_str());
      }
      
      // set the preconditioner
      string precon = "";
      Parameters::get(name + "->pc_type", precon);
      if (!precon.size())
	Parameters::get(name + "->left precon", precon);
69 70
      if (!precon.size())
	Parameters::get(name + "->right precon", precon);
71 72 73 74 75 76 77 78 79 80 81 82 83
      if (!matSolverPackage && params.emptyParam.find(precon) == params.emptyParam.end()) {
	precon = (params.preconMap.find(precon) != params.preconMap.end() ? params.preconMap[precon] : precon);
	PetscOptionsInsertString(("-" + kspPrefix + "pc_type " + precon).c_str());
      }
      
      PetscOptionsInsertString(("-" + kspPrefix + "ksp_max_it " + boost::lexical_cast<std::string>(getMaxIterations())).c_str());
      PetscOptionsInsertString(("-" + kspPrefix + "ksp_rtol " + boost::lexical_cast<std::string>(getRelative())).c_str());    
      PetscOptionsInsertString(("-" + kspPrefix + "ksp_atol " + boost::lexical_cast<std::string>(getTolerance())).c_str());   
      
      if (getInfo() >= 20)
	PetscOptionsInsertString(("-" + kspPrefix + "ksp_monitor_true_residual").c_str());
      else if (getInfo() >= 10)
	PetscOptionsInsertString(("-" + kspPrefix + "ksp_monitor").c_str());
84 85
    }
    if (!matSolverPackage) {
86
      Parameters::get(name + "->use zero start vector", zeroStartVector);
87
    }
88
    Parameters::get("parallel->print matrix info", printMatInfo);
89 90 91 92 93 94 95 96 97 98 99 100 101
    
#if DEBUG != 0
    bool printOptionsInfo = false;
    Parameters::get("parallel->debug->print options info", printOptionsInfo);
    if (printOptionsInfo) {
    MSG("PetscOptionsView:\n");
    PetscViewer viewer;
    PetscViewerCreate(PETSC_COMM_WORLD, &viewer);
    PetscViewerSetType(viewer, PETSCVIEWERASCII);
    PetscOptionsView(viewer);
    PetscViewerDestroy(&viewer);
    }
#endif
102 103 104
  }


Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
105
  void PetscSolverGlobalMatrix::fillPetscMatrix(Matrix<DOFMatrix*> *seqMat)
Thomas Witkowski's avatar
Thomas Witkowski committed
106 107 108
  {
    FUNCNAME("PetscSolverGlobalMatrix::fillPetscMatrix()");

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
109 110 111 112
    TEST_EXIT_DBG(meshDistributor)("No mesh distributor object defined!\n");
    TEST_EXIT_DBG(interiorMap)("No parallel mapping object defined!\n");
    TEST_EXIT_DBG(seqMat)("No DOF matrix defined!\n");
    
113
#if (DEBUG != 0)
114
    Timer t;
115
#endif
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
116

117
    createMatVec(*seqMat);
118

Thomas Witkowski's avatar
Thomas Witkowski committed
119
    if (coarseSpaceMap.size()) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
120
      fillPetscMatrixWithCoarseSpace(seqMat);
121 122
      return;
    }
123
    
124 125
    // === Create PETSc vector (solution and a temporary vector). ===

Thomas Witkowski's avatar
Thomas Witkowski committed
126
#if (DEBUG != 0)
127 128
    MSG("Fill petsc matrix 1 needed %.5f seconds\n", t.elapsed());
    t.reset();
Thomas Witkowski's avatar
Thomas Witkowski committed
129 130 131 132
#endif

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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
133
    int nComponents = seqMat->getNumRows();
Thomas Witkowski's avatar
Thomas Witkowski committed
134 135
    for (int i = 0; i < nComponents; i++)
      for (int j = 0; j < nComponents; j++)
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
136 137
	if ((*seqMat)[i][j])
	  setDofMatrix((*seqMat)[i][j], i, j);
Thomas Witkowski's avatar
Thomas Witkowski committed
138 139

#if (DEBUG != 0)
140 141
    MSG("Fill petsc matrix 2 needed %.5f seconds\n", t.elapsed());
    t.reset();
Thomas Witkowski's avatar
Thomas Witkowski committed
142 143
#endif

144
    matAssembly();
145

146 147
    if (printMatInfo) {
      MatInfo matInfo;
148
      MatGetInfo(getMatInterior(), MAT_GLOBAL_SUM, &matInfo);
149 150 151 152 153 154 155
      MSG("Matrix info:\n");
      MSG("  memory usage: %e MB\n", matInfo.memory / (1024.0 * 1024.0));
      MSG("  mallocs: %d\n", static_cast<int>(matInfo.mallocs));
      MSG("  nz allocated: %d\n", static_cast<int>(matInfo.nz_allocated));
      MSG("  nz used: %d\n", static_cast<int>(matInfo.nz_used));
      MSG("  nz unneeded: %d\n", static_cast<int>(matInfo.nz_unneeded));
    }
156 157


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

160 161
    initSolver(kspInterior);
    KSPGetPC(kspInterior, &pcInterior);
Praetorius, Simon's avatar
Praetorius, Simon committed
162
    initPreconditioner(*seqMat, mat[0][0]);
163

164

Thomas Witkowski's avatar
Thomas Witkowski committed
165
#if (DEBUG != 0)
166
    MSG("Fill petsc matrix 3 needed %.5f seconds\n", t.elapsed());
Thomas Witkowski's avatar
Thomas Witkowski committed
167
#endif
168 169


170 171 172 173 174 175 176 177 178 179 180
    // === For debugging allow to write the matrix to a file. ===

    bool dbgWriteMatrix = false;
    Parameters::get("parallel->debug->write matrix", dbgWriteMatrix);
    if (dbgWriteMatrix) {
      PetscViewer matView;
      PetscViewerBinaryOpen(PETSC_COMM_WORLD, "mpi.mat",
			    FILE_MODE_WRITE, &matView);
      MatView(getMatInterior(), matView);
      PetscViewerDestroy(&matView);
    }
181 182 183
  }


Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
184
  void PetscSolverGlobalMatrix::fillPetscMatrixWithCoarseSpace(Matrix<DOFMatrix*> *seqMat)
185
  {
186
    FUNCNAME_DBG("PetscSolverGlobalMatrix::fillPetscMatrixWithCoarseSpace()");
187

188
    TEST_EXIT_DBG(interiorMap)("No interiorMap! Should not happen!\n");
189
    TEST_EXIT_DBG(coarseSpaceMap.size() == (size_t)seqMat->getSize())
190
      ("Wrong sizes %d %d\n", coarseSpaceMap.size(), seqMat->getSize());
191

192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207
    // === Prepare traverse of sequentially created matrices. ===

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

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

    vector<int> cols, colsOther;
    vector<double> values, valuesOther;
    cols.reserve(300);
    colsOther.reserve(300);
    values.reserve(300);
    valuesOther.reserve(300);

208 209 210
    bool localMatrix = 
      (meshDistributor->getMeshLevelData().getMpiComm(meshLevel) == MPI::COMM_SELF);

211 212 213
    // === Traverse all sequentially created matrices and add the values to ===
    // === the global PETSc matrices.                                       ===

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
214
    int nComponents = seqMat->getSize();
Thomas Witkowski's avatar
Thomas Witkowski committed
215 216 217 218 219
    for (int rowComponent = 0; rowComponent < nComponents; rowComponent++) {
      for (int colComponent = 0; colComponent < nComponents; colComponent++) {
	DOFMatrix* dofMat = (*seqMat)[rowComponent][colComponent];

	if (!dofMat)
220 221
	  continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
222 223
	ParallelDofMapping *rowCoarseSpace = coarseSpaceMap[rowComponent];
	ParallelDofMapping *colCoarseSpace = coarseSpaceMap[colComponent];
Thomas Witkowski's avatar
Thomas Witkowski committed
224

Thomas Witkowski's avatar
Thomas Witkowski committed
225 226
	std::set<DegreeOfFreedom> &dirichletRows = dofMat->getDirichletRows();

Thomas Witkowski's avatar
Thomas Witkowski committed
227 228
	traits::col<Matrix>::type col(dofMat->getBaseMatrix());
	traits::const_value<Matrix>::type value(dofMat->getBaseMatrix());
229 230
	
	// Traverse all rows.
Thomas Witkowski's avatar
Thomas Witkowski committed
231 232
	for (cursor_type cursor = begin<row>(dofMat->getBaseMatrix()), 
	       cend = end<row>(dofMat->getBaseMatrix()); cursor != cend; ++cursor) {
233

234
	  bool isRowCoarse = isCoarseSpace(rowComponent, cursor.value());
Thomas Witkowski's avatar
Thomas Witkowski committed
235 236 237

	  // For the case, this is a dirichlet row we have to check whether the 
	  // rank is also owner of this row DOF.
238 239 240
	  if (dirichletRows.count(cursor.value())) {
	    if ((!isRowCoarse && !(*interiorMap)[rowComponent].isRankDof(cursor.value())) ||
		(isRowCoarse && !(*rowCoarseSpace)[rowComponent].isRankDof(cursor.value())))
241
	      continue;    
Thomas Witkowski's avatar
Thomas Witkowski committed
242
	  }
243 244 245 246 247 248 249 250 251 252
  
	  cols.clear();
	  colsOther.clear();
	  values.clear();	  
	  valuesOther.clear();

	  // Traverse all columns.
	  for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); 
	       icursor != icend; ++icursor) {

Thomas Witkowski's avatar
Thomas Witkowski committed
253
	    bool isColCoarse = isCoarseSpace(colComponent, col(*icursor));
254

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
255
	    if (isColCoarse == false)
256
	      if ((*interiorMap)[colComponent].isSet(col(*icursor)) == false)
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
257 258
		continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
259 260 261
	    if (isColCoarse == isRowCoarse) {
	      cols.push_back(col(*icursor));
	      values.push_back(value(*icursor));
262
	    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
263 264
	      colsOther.push_back(col(*icursor));
	      valuesOther.push_back(value(*icursor));
265 266 267 268 269 270
	    }
	  }  // for each nnz in row


	  // === Set matrix values. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
271
	  if (isRowCoarse) {
Thomas Witkowski's avatar
Thomas Witkowski committed
272 273
	    for (unsigned int i = 0; i < cols.size(); i++)
	      cols[i] = colCoarseSpace->getMatIndex(colComponent, cols[i]);
274

275
	    int rowIndex = rowCoarseSpace->getMatIndex(rowComponent, cursor.value());
Thomas Witkowski's avatar
Thomas Witkowski committed
276
	    MatSetValues(getMatCoarseByComponent(rowComponent, colComponent),
277 278
	    		 1, &rowIndex, cols.size(),
	    		 &(cols[0]), &(values[0]), ADD_VALUES);
279 280

	    if (colsOther.size()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
281 282 283
	      for (unsigned int i = 0; i < colsOther.size(); i++)
		colsOther[i] = 
		  interiorMap->getMatIndex(colComponent, colsOther[i]) + 
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
284
		  rStartInterior;	      
285

Thomas Witkowski's avatar
Thomas Witkowski committed
286
	      MatSetValues(getMatCoarseInteriorByComponent(rowComponent), 
287 288
	       		   1, &rowIndex, colsOther.size(),
 	       		   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
289 290
	    }
	  } else {
291
	    if ((*interiorMap)[rowComponent].isSet(cursor.value()) == false)
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
292 293
	      continue;

294
	    int localRowIndex = 
295
	      (localMatrix ? 
296 297
	       interiorMap->getLocalMatIndex(rowComponent, cursor.value()) :
	       interiorMap->getMatIndex(rowComponent, cursor.value()));
298

Thomas Witkowski's avatar
Thomas Witkowski committed
299
	    for (unsigned int i = 0; i < cols.size(); i++) {
300
	      if (localMatrix)
Thomas Witkowski's avatar
Thomas Witkowski committed
301
		cols[i] = interiorMap->getLocalMatIndex(colComponent, cols[i]);
302
	      else
Thomas Witkowski's avatar
Thomas Witkowski committed
303
		cols[i] = interiorMap->getMatIndex(colComponent, cols[i]);
304 305
	    }
	    
306 307
	    MatSetValues(getMatInterior(), 1, &localRowIndex, cols.size(),
  	     		 &(cols[0]), &(values[0]), ADD_VALUES);
Thomas Witkowski's avatar
Thomas Witkowski committed
308
	    
309
	    if (colsOther.size()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
310
	      int globalRowIndex = 
311
		interiorMap->getMatIndex(rowComponent, cursor.value()) + rStartInterior;
Thomas Witkowski's avatar
Thomas Witkowski committed
312
      
Thomas Witkowski's avatar
Thomas Witkowski committed
313 314 315
	      for (unsigned int i = 0; i < colsOther.size(); i++)
		colsOther[i] = 
		  colCoarseSpace->getMatIndex(colComponent, colsOther[i]);
316

317 318 319
	      MatSetValues(getMatInteriorCoarseByComponent(colComponent), 
			   1, &globalRowIndex, colsOther.size(),
			   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
320 321 322 323 324 325
	    }
	  }
	} 
      }
    }

326
    matAssembly();
327

328 329
    // === Create solver for the non primal (thus local) variables. ===

330
    KSPCreate(domainComm, &kspInterior);
331 332 333 334 335
#if (PETSC_VERSION_MINOR >= 5)
    KSPSetOperators(kspInterior, getMatInterior(), getMatInterior());
#else
    KSPSetOperators(kspInterior, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN);
#endif
336 337 338
    KSPSetOptionsPrefix(kspInterior, "interior_");
    KSPSetType(kspInterior, KSPPREONLY);
    KSPGetPC(kspInterior, &pcInterior);
339 340
    if (isSymmetric) {
      PCSetType(pcInterior, PCCHOLESKY);
341
      PCFactorSetMatSolverPackage(pcInterior, MATSOLVERMUMPS);
342 343
    } else {
      PCSetType(pcInterior, PCLU);
344
      if (localMatrix)
345 346 347 348
	PCFactorSetMatSolverPackage(pcInterior, MATSOLVERUMFPACK);
      else
	PCFactorSetMatSolverPackage(pcInterior, MATSOLVERMUMPS);
    }
349
    KSPSetFromOptions(kspInterior);
350 351 352
  }


353 354
  void PetscSolverGlobalMatrix::fillPetscRhs(SystemVector *vec)
  {
355
    FUNCNAME_DBG("PetscSolverGlobalMatrix::fillPetscRhs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
356

357
    TEST_EXIT_DBG(vec)("No DOF vector defined!\n");
358
    TEST_EXIT_DBG(interiorMap)("No parallel DOF map defined!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
359
    
Thomas Witkowski's avatar
Thomas Witkowski committed
360
    // === Transfer values from DOF vector to the PETSc vector. === 
Thomas Witkowski's avatar
Thomas Witkowski committed
361
    if (coarseSpaceMap.size()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
362
      for (int i = 0; i < vec->getSize(); i++)
363
	setDofVector(getVecRhsInterior(), 
Thomas Witkowski's avatar
Thomas Witkowski committed
364
		     getVecRhsCoarseByComponent(i), vec->getDOFVector(i), i);
Thomas Witkowski's avatar
Thomas Witkowski committed
365 366
    } else {
      for (int i = 0; i < vec->getSize(); i++)
367
	setDofVector(getVecRhsInterior(), vec->getDOFVector(i), i);
Thomas Witkowski's avatar
Thomas Witkowski committed
368 369
    }

370
    vecRhsAssembly();
371 372 373 374 375 376 377 378 379 380 381 382

    // === For debugging allow to write the rhs vector to a file. ===

    bool dbgWriteRhs = false;
    Parameters::get("parallel->debug->write rhs", dbgWriteRhs);
    if (dbgWriteRhs) {
      PetscViewer vecView;
      PetscViewerBinaryOpen(PETSC_COMM_WORLD, "mpi.vec",
			    FILE_MODE_WRITE, &vecView);
      VecView(getVecRhsInterior(), vecView);
      PetscViewerDestroy(&vecView);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
383 384 385
  }


386 387 388 389 390
  /// 1.) set startsolution
  /// 2.) create null-space
  /// 3.) solve Ax=b
  /// 4.) destroy null-space
  /// 5.) transfer solution back to DOFVector
391 392
  void PetscSolverGlobalMatrix::solvePetscMatrix(SystemVector &vec, 
						 AdaptInfo *adaptInfo)
Thomas Witkowski's avatar
Thomas Witkowski committed
393 394 395 396 397 398 399
  {
    FUNCNAME("PetscSolverGlobalMatrix::solvePetscMatrix()");

    int nComponents = vec.getSize();

    // === Set old solution to be initiual guess for PETSc solver. ===
    if (!zeroStartVector) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
400 401
      TEST_EXIT(coarseSpaceMap.size() == 0)("Not yet supported!\n");

402
      VecSet(getVecSolInterior(), 0.0);
Thomas Witkowski's avatar
Thomas Witkowski committed
403 404
      
      for (int i = 0; i < nComponents; i++)
405
	setDofVector(getVecSolInterior(), vec.getDOFVector(i), i, true);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
406

407
      vecSolAssembly();
Thomas Witkowski's avatar
Thomas Witkowski committed
408 409
    }

410

Thomas Witkowski's avatar
Thomas Witkowski committed
411 412
    MatNullSpace matNullspace;
    Vec nullspaceBasis;
413 414 415
    if (nullspace.size() > 0 || 
	hasConstantNullspace ||
	constNullspaceComponent.size() > 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
416
      TEST_EXIT_DBG(nullspace.size() <= 1)("Not yet implemented!\n");
417

418 419 420 421 422 423 424 425 426 427
      if (constNullspaceComponent.size() > 0) {
	nullspace.clear();
	SystemVector *basisVec = new SystemVector(vec);
	basisVec->set(0.0);
	for (unsigned int i = 0; i < constNullspaceComponent.size(); i++)
	  basisVec->getDOFVector(constNullspaceComponent[i])->set(1.0);

	nullspace.push_back(basisVec);
      } 

Thomas Witkowski's avatar
Thomas Witkowski committed
428
      if (nullspace.size() > 0) {
429
	VecDuplicate(getVecSolInterior(), &nullspaceBasis);
430
	setDofVector(nullspaceBasis, *(nullspace[0]), true);
Thomas Witkowski's avatar
Thomas Witkowski committed
431 432 433
	
	VecAssemblyBegin(nullspaceBasis);
	VecAssemblyEnd(nullspaceBasis);
434 435

	VecNormalize(nullspaceBasis, PETSC_NULL);
Thomas Witkowski's avatar
Thomas Witkowski committed
436
	
437
	MatNullSpaceCreate(domainComm, (hasConstantNullspace ? PETSC_TRUE : PETSC_FALSE), 
Thomas Witkowski's avatar
Thomas Witkowski committed
438 439
			   1, &nullspaceBasis, &matNullspace);

440
	MatMult(getMatInterior(), nullspaceBasis, getVecSolInterior());
Thomas Witkowski's avatar
Thomas Witkowski committed
441
	PetscReal n;
442
	VecNorm(getVecSolInterior(), NORM_2, &n);
Thomas Witkowski's avatar
Thomas Witkowski committed
443 444
	MSG("NORM IS: %e\n", n);
      } else {
445
	MatNullSpaceCreate(domainComm, PETSC_TRUE, 0, PETSC_NULL, &matNullspace);
Thomas Witkowski's avatar
Thomas Witkowski committed
446
      }
447

Thomas Witkowski's avatar
Thomas Witkowski committed
448 449 450
      MSG("NULLSPACE IS NOT REMOVED!\n");
      // MatSetNullSpace(getMatInterior(), matNullspace);
      // KSPSetNullSpace(kspInterior, matNullspace);
451

Thomas Witkowski's avatar
Thomas Witkowski committed
452
      // === Remove null space, if requested. ===
453

Thomas Witkowski's avatar
Thomas Witkowski committed
454
      if (removeRhsNullspace) {
Thomas Witkowski's avatar
Thomas Witkowski committed
455
	TEST_EXIT_DBG(coarseSpaceMap.empty())("Not supported!\n");
456 457

	MSG("Remove nullspace from rhs vector.\n");
458 459 460
#if (PETSC_VERSION_MINOR >= 5)
	MatNullSpaceRemove(matNullspace, getVecRhsInterior());
#else
461
	MatNullSpaceRemove(matNullspace, getVecRhsInterior(), PETSC_NULL);
462
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
463 464 465 466
      }
    } else {
      TEST_EXIT(removeRhsNullspace == false)
	("No nullspace provided that should be removed from rhs!\n");
467 468
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
469
    // PETSc.
470
    solve(getVecRhsInterior(), getVecSolInterior());
Thomas Witkowski's avatar
Thomas Witkowski committed
471

472 473

    if (nullspace.size() > 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
474
      MatNullSpaceDestroy(&matNullspace);
475 476 477 478
      VecDestroy(&nullspaceBasis);
    }


Thomas Witkowski's avatar
Thomas Witkowski committed
479 480
    // === Transfere values from PETSc's solution vectors to the DOF vectors. ===
    PetscScalar *vecPointer;
481
    VecGetArray(getVecSolInterior(), &vecPointer);    
Thomas Witkowski's avatar
Thomas Witkowski committed
482

483
    int c = 0;
484 485
    for (int component = 0; component < nComponents; component++) {
      DOFVector<double> &dv = *(vec.getDOFVector(component));
486

487
      DofMap& d = (*interiorMap)[component].getMap();
488 489 490
      for (DofMap::iterator it = d.begin(); it != d.end(); ++it)
	if (it->second.local != -1)
	  dv[it->first] = vecPointer[c++];
Thomas Witkowski's avatar
Thomas Witkowski committed
491 492
    }

493
    VecRestoreArray(getVecSolInterior(), &vecPointer);
Thomas Witkowski's avatar
Thomas Witkowski committed
494 495 496 497

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


501 502 503
  void PetscSolverGlobalMatrix::solveGlobal(Vec &rhs, Vec &sol)
  {
    Vec tmp;
504
    if (domainComm.Get_size() == 1)
505
      createLocalVec(*interiorMap, tmp);
506
    else
507
      createVec(*interiorMap, tmp);
508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533

    PetscScalar *tmpValues, *rhsValues;
    VecGetArray(tmp, &tmpValues);
    VecGetArray(rhs, &rhsValues);

    for (int i = 0; i < interiorMap->getRankDofs(); i++)
      tmpValues[i] = rhsValues[i];

    VecRestoreArray(rhs, &rhsValues);
    VecRestoreArray(tmp, &tmpValues);

    KSPSolve(kspInterior, tmp, tmp);

    VecGetArray(tmp, &tmpValues);
    VecGetArray(sol, &rhsValues);

    for (int i = 0; i < interiorMap->getRankDofs(); i++) 
      rhsValues[i] = tmpValues[i];

    VecRestoreArray(sol, &rhsValues);
    VecRestoreArray(tmp, &tmpValues);

    VecDestroy(&tmp);
  }


534 535
  void PetscSolverGlobalMatrix::destroyMatrixData()
  {
536
    matDestroy();
Thomas Witkowski's avatar
Thomas Witkowski committed
537

538 539
    exitPreconditioner(pcInterior);
    exitSolver(kspInterior);
Thomas Witkowski's avatar
Thomas Witkowski committed
540 541 542
  }


543 544
  void PetscSolverGlobalMatrix::destroyVectorData()
  {
545
    vecDestroy();
546 547 548
  }


549 550 551 552 553 554 555 556 557 558 559 560 561 562 563
  void PetscSolverGlobalMatrix::createFieldSplit(PC pc)
  {
    FUNCNAME("PetscSolverGlobalMatrix::createFieldSplit()");

    vector<string> isNames;
    Parameters::get("parallel->solver->is blocks", isNames);

    int nBlocks = isNames.size();
    if (nBlocks == 0)
      return;

    for (int i = 0; i < nBlocks; i++) {
      MSG("Create for block %s\n", isNames[i].c_str());

      vector<int> blockComponents;
564
      Parameters::get("parallel->solver->is block " + boost::lexical_cast<string>(i),
565 566 567
		      blockComponents);
      int nComponents = static_cast<int>(blockComponents.size());

568
      TEST_EXIT(nComponents > 0)("No IS block for block %d defined!\n", i);
569 570 571 572 573 574 575 576

      // Check if blocks are continous
      for (int j = 0; j < nComponents; j++) {
	TEST_EXIT(blockComponents[j] == blockComponents[0] + j)
	  ("Does not yet support not continous IS blocks! Block %s\n", 
	   isNames[i].c_str());
      }

577
      createFieldSplit(pc, isNames[i].c_str(), blockComponents);
578 579 580 581
    }
  }


582
  void PetscSolverGlobalMatrix::createFieldSplit(PC pc, 
583
						 const char* splitName, 
584 585 586 587
						 vector<int> &components)
  {
    IS is;
    interiorMap->createIndexSet(is, components[0], components.size());
588
    PCFieldSplitSetIS(pc, splitName, is);
589 590
    ISDestroy(&is);
  }
591 592 593 594 595 596 597 598 599
  
  
  void PetscSolverGlobalMatrix::extractVectorComponent(Vec input, int i, Vec *output, int numberOfComponents)
  {
    IS is;
    interiorMap->createIndexSet(is, i, numberOfComponents);
    VecGetSubVector(input, is, output);
    ISDestroy(&is);
  }
600

601 602 603 604 605 606 607 608 609 610
  void PetscSolverGlobalMatrix::extractMatrixComponent(Mat input, int startRow, int numberOfRows, int startCol, int numberOfCols, Mat *output)
  {
    IS isrow, iscol;
    interiorMap->createIndexSet(isrow, startRow, numberOfRows);
    interiorMap->createIndexSet(iscol, startCol, numberOfCols);
    MatGetSubMatrix(input, isrow, iscol, MAT_INITIAL_MATRIX, output);
    ISDestroy(&iscol);
    ISDestroy(&isrow);
  }
  
611

612
  void PetscSolverGlobalMatrix::initSolver(KSP &ksp)
613
  {
614
    KSPCreate(domainComm, &ksp);
615 616 617 618 619
#if (PETSC_VERSION_MINOR >= 5)
    KSPSetOperators(ksp, getMatInterior(), getMatInterior());
#else
    KSPSetOperators(ksp, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN);
#endif
620 621 622
    KSPSetTolerances(ksp, 0.0, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT);
    KSPSetType(ksp, KSPBCGS);
    KSPSetOptionsPrefix(ksp, kspPrefix.c_str());
623
    MatSetOptionsPrefix(getMatInterior(), kspPrefix.c_str());
624 625 626 627 628 629 630 631
    KSPSetFromOptions(ksp);

    // Do not delete the solution vector, use it for the initial guess.
    if (!zeroStartVector)
      KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
  }  


632
  void PetscSolverGlobalMatrix::exitSolver(KSP &ksp)
633 634 635 636 637
  {
    KSPDestroy(&ksp);
  }  


638 639 640 641 642 643 644 645
  void PetscSolverGlobalMatrix::initPreconditioner(PC pc)
  {
    PCSetFromOptions(pc);
    createFieldSplit(pc);
  }

 
  void PetscSolverGlobalMatrix::exitPreconditioner(PC pc)
646
  { }
647 648


649
  /// called by \ref fillPetscMatrix
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
650
  void PetscSolverGlobalMatrix::setDofMatrix(DOFMatrix* seqMat,
Thomas Witkowski's avatar
Thomas Witkowski committed
651
					     int rowComp, int colComp)
Thomas Witkowski's avatar
Thomas Witkowski committed
652 653 654
  {
    FUNCNAME("PetscSolverGlobalMatrix::setDofMatrix()");

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
655
    TEST_EXIT(seqMat)("No DOFMatrix!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
656 657

    using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
658
    namespace traits = mtl::traits;
Thomas Witkowski's avatar
Thomas Witkowski committed
659 660
    typedef DOFMatrix::base_matrix_type Matrix;

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
661 662
    traits::col<Matrix>::type col(seqMat->getBaseMatrix());
    traits::const_value<Matrix>::type value(seqMat->getBaseMatrix());
Thomas Witkowski's avatar
Thomas Witkowski committed
663 664 665 666 667 668 669 670 671 672 673

    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);
    
    vector<int> globalCols;

674 675
    // Get periodic mapping object
    PeriodicMap &perMap = meshDistributor->getPeriodicMap();
Thomas Witkowski's avatar
Thomas Witkowski committed
676 677
    std::set<DegreeOfFreedom> &dirichletRows = seqMat->getDirichletRows();
      
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
678 679
    const FiniteElemSpace *rowFe = seqMat->getRowFeSpace();
    const FiniteElemSpace *colFe = seqMat->getColFeSpace();
Thomas Witkowski's avatar
Thomas Witkowski committed
680 681

    // === Traverse all rows of the DOF matrix and insert row wise the values ===
Thomas Witkowski's avatar
Thomas Witkowski committed
682 683
    // === to the PETSc matrix.                                               ===

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
684 685
    for (cursor_type cursor = begin<row>(seqMat->getBaseMatrix()), 
	   cend = end<row>(seqMat->getBaseMatrix()); cursor != cend; ++cursor) {
Thomas Witkowski's avatar
Thomas Witkowski committed
686
      // Global index of the current row DOF.
687
      MultiIndex rowMultiIndex;
688
      if ((*interiorMap)[rowComp].find(cursor.value(), rowMultiIndex) == false)
689 690 691
	continue;

      int globalRowDof = rowMultiIndex.global;
692

Thomas Witkowski's avatar
Thomas Witkowski committed
693
      // Test if the current row DOF is a periodic DOF.
694
      bool periodicRow = perMap.isPeriodic(rowFe, globalRowDof);
695

Thomas Witkowski's avatar
Thomas Witkowski committed
696
      // Dirichlet rows can be set only be the owner ranks.
697
      if (dirichletRows.count(cursor.value()) && !((*interiorMap)[rowComp].isRankDof(cursor.value())))
Thomas Witkowski's avatar
Thomas Witkowski committed
698 699
	continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
700 701 702
      if (!periodicRow) {
	// === Row DOF index is not periodic. ===

703
	// Get PETSc's mat row index.
Thomas Witkowski's avatar
Thomas Witkowski committed
704
	int rowIndex = interiorMap->getMatIndex(rowComp, globalRowDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
705 706 707 708 709 710 711 712

	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.
713
	  MultiIndex colMultiIndex;
Thomas Witkowski's avatar
Thomas Witkowski committed
714
	  if ((*interiorMap)[colComp].find(col(*icursor), colMultiIndex) == false)
715 716 717
	    continue;

	  int globalColDof = colMultiIndex.global;
Thomas Witkowski's avatar
Thomas Witkowski committed
718
	  // Test if the current col dof is a periodic dof.
719
	  bool periodicCol = perMap.isPeriodic(colFe, globalColDof);
720
	  // Get PETSc's mat col index.
Thomas Witkowski's avatar
Thomas Witkowski committed
721
	  int colIndex = interiorMap->getMatIndex(colComp, globalColDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
722 723 724 725 726 727 728

	  // Ignore all zero entries, expect it is a diagonal entry.
 	  if (value(*icursor) == 0.0 && rowIndex != colIndex)
 	    continue;

	  if (!periodicCol) {
	    // Calculate the exact position of the column index in the PETSc matrix.
Thomas Witkowski's avatar
Thomas Witkowski committed
729 730
 	    cols.push_back(colIndex);
 	    values.push_back(value(*icursor));
Thomas Witkowski's avatar
Thomas Witkowski committed
731 732 733 734 735
	  } else {
	    // === Row index is not periodic, but column index is. ===

	    // Create set of all periodic associations of the column index.
	    std::set<int> perAsc;
736 737 738
	    perMap.fillAssociations(colFe, globalColDof, 
				    meshDistributor->getElementObjectDb(), perAsc);

Thomas Witkowski's avatar
Thomas Witkowski committed
739 740 741 742 743 744 745 746 747
	    // Scale value to the number of periodic associations of the column index.
	    double scaledValue = 
	      value(*icursor) * pow(0.5, static_cast<double>(perAsc.size()));

	    
	    // === Create set of all matrix column indices due to the periodic ===
	    // === associations of the column DOF index.                       ===

	    vector<int> newCols;
748
	    perMap.mapDof(colFe, globalColDof, perAsc, newCols);
Thomas Witkowski's avatar
Thomas Witkowski committed
749
	    for (unsigned int i = 0; i < newCols.size(); i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
750
	      cols.push_back(interiorMap->getMatIndex(colComp, newCols[i]));
Thomas Witkowski's avatar
Thomas Witkowski committed
751 752 753 754 755
	      values.push_back(scaledValue);	      
	    }
	  }
	}

756
  	MatSetValues(getMatInterior(), 1, &rowIndex, cols.size(), 
Thomas Witkowski's avatar
Thomas Witkowski committed
757
  		     &(cols[0]), &(values[0]), ADD_VALUES);	
Thomas Witkowski's avatar
Thomas Witkowski committed
758 759 760 761 762 763 764 765 766 767 768 769 770 771
      } else {
	// === Row DOF index is periodic. ===

	// Because this row is periodic, we will have to add the entries of this 
	// matrix row to multiple rows. The following maps store to each row an
	// array of column indices and values of the entries that must be added to
	// the PETSc matrix.
	map<int, vector<int> > colsMap;
	map<int, vector<double> > valsMap;

	// Traverse all column entries.
	for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); 
	     icursor != icend; ++icursor) {
	  // Global index of the current column index.
Thomas Witkowski's avatar
Thomas Witkowski committed
772
	  int globalColDof = (*interiorMap)[colComp][col(*icursor)].global;
Thomas Witkowski's avatar
Thomas Witkowski committed
773 774 775 776 777 778 779 780 781

	  // Ignore all zero entries, expect it is a diagonal entry.
 	  if (value(*icursor) == 0.0 && globalRowDof != globalColDof)
 	    continue;

	  // === Add all periodic associations of both, the row and the column ===
	  // === indices to the set perAsc.                                    ===

	  std::set<int> perAsc;
782 783 784 785
	  perMap.fillAssociations(colFe, globalColDof, 
				  meshDistributor->getElementObjectDb(), perAsc);
	  perMap.fillAssociations(rowFe, globalRowDof, 
				  meshDistributor->getElementObjectDb(), perAsc);
Thomas Witkowski's avatar
Thomas Witkowski committed
786 787 788 789 790 791 792 793 794 795

	  // Scale the value with respect to the number of periodic associations.
	  double scaledValue = 
	    value(*icursor) * pow(0.5, static_cast<double>(perAsc.size()));


	  // === Create all matrix entries with respect to the periodic  ===
	  // === associations of the row and column indices.             ===

	  vector<pair<int, int> > entry;
Thomas Witkowski's avatar
Blbu  
Thomas Witkowski committed
796 797
	  perMap.mapDof(rowFe, colFe, make_pair(globalRowDof, globalColDof),
			perAsc, entry);
Thomas Witkowski's avatar
Thomas Witkowski committed
798 799 800

	  // === Translate the matrix entries to PETSc's matrix.

801
	  for (unsigned int i = 0; i < entry.size(); i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
802 803
	    int rowIdx = interiorMap->getMatIndex(rowComp, entry[i].first);
	    int colIdx = interiorMap->getMatIndex(colComp, entry[i].second);
Thomas Witkowski's avatar
Thomas Witkowski committed
804

805 806
	    colsMap[rowIdx].push_back(colIdx);
	    valsMap[rowIdx].push_back(scaledValue);
Thomas Witkowski's avatar
Thomas Witkowski committed
807 808 809 810 811 812 813 814 815 816 817 818
	  }
	}


	// === Finally, add all periodic rows to the PETSc matrix. ===

	for (map<int, vector<int> >::iterator rowIt = colsMap.begin();
	     rowIt != colsMap.end(); ++rowIt) {
	  TEST_EXIT_DBG(rowIt->second.size() == valsMap[rowIt->first].size())
	    ("Should not happen!\n");

	  int rowIndex = rowIt->first;
819

820
	  MatSetValues(getMatInterior(), 1, &rowIndex, rowIt->second.size(),
Thomas Witkowski's avatar
Thomas Witkowski committed
821 822 823 824 825 826 827
		       &(rowIt->second[0]), &(valsMap[rowIt->first][0]), ADD_VALUES);
	}
      }
    }
  }


828
  /// called by \ref fillPetscRhs
Thomas Witkowski's avatar
Thomas Witkowski committed
829 830
  void PetscSolverGlobalMatrix::setDofVector(Vec vecInterior, 
					     Vec vecCoarse,
831
					     DOFVector<double>* vec, 
Thomas Witkowski's avatar
Thomas Witkowski committed
832
					     int rowComp, 
833
					     bool rankOnly)
Thomas Witkowski's avatar
Thomas Witkowski committed
834
  {
835
    FUNCNAME_DBG("PetscSolverGlobalMatrix::setDofVector()");
Thomas Witkowski's avatar
Thomas Witkowski committed
836

837
    const FiniteElemSpace *feSpace = vec->getFeSpace();
838
    PeriodicMap &perMap = meshDistributor->getPeriodicMap();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
839

840
    ParallelDofMapping *rowCoarseSpace = 
841
      (coarseSpaceMap.size() ? coarseSpaceMap[rowComp] : NULL);
Thomas Witkowski's avatar
Thomas Witkowski committed
842 843

    map<DegreeOfFreedom, double> &dirichletValues = vec->getDirichletValues();
844

Thomas Witkowski's avatar
Thomas Witkowski committed
845 846 847
    // Traverse all used DOFs in the dof vector.
    DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
    for (dofIt.reset(); !dofIt.end(); ++dofIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
848 849
      
      DegreeOfFreedom dof = dofIt.getDOFIndex();
850

Thomas Witkowski's avatar
Thomas Witkowski committed
851 852
      if (rankOnly && !(*interiorMap)[rowComp].isRankDof(dof))
	continue;
853

854 855
      bool isCoarseDof = isCoarseSpace(rowComp, dof);

Thomas Witkowski's avatar
Thomas Witkowski committed
856
      // Dirichlet rows can be set only be the owner ranks.
857 858 859 860 861
      if (dirichletValues.count(dof)) {
	if ((!isCoarseDof && !((*interiorMap)[rowComp].isRankDof(dof))) ||
	    (isCoarseDof && !((*rowCoarseSpace)[rowComp].isRankDof(dof))))
	  continue;
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
862

863
      if (isCoarseDof) {
864
	TEST_EXIT_DBG(vecCoarse != PETSC_NULL)("vecCoarse not set! Should not happen!\n");
865

Thomas Witkowski's avatar
Thomas Witkowski committed
866
	int index = rowCoarseSpace->getMatIndex(rowComp, dof);
867 868
	VecSetValue(vecCoarse, index, *dofIt, ADD_VALUES);
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
869
	if ((*interiorMap)[rowComp].isSet(dof) == false)
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
870 871
	  continue;

872
	// Calculate global row index of the DOF.
Thomas Witkowski's avatar
Thomas Witkowski committed
873
	DegreeOfFreedom globalRowDof = (*interiorMap)[rowComp][dof].global;
874 875 876 877
	
	// Get PETSc's mat index of the row DOF.
	int index = 0;
	if (interiorMap->isMatIndexFromGlobal())
878
	  index = interiorMap->getMatIndex(rowComp, globalRowDof) + rStartInterior;
879 880
	else
	  index =
Thomas Witkowski's avatar
Thomas Witkowski committed
881
	    interiorMap->getMatIndex(rowComp, dof) + rStartInterior;
882

883
	if (perMap.isPeriodic(feSpace, globalRowDof)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
884 885 886
	  std::set<int>& perAsc = perMap.getAssociations(feSpace, globalRowDof);
	  double value = *dofIt / (perAsc.size() + 1.0);
	  VecSetValue(vecInterior, index, value, ADD_VALUES);
887

Thomas Witkowski's avatar
Thomas Witkowski committed
888 889 890
	  for (std::set<int>::iterator perIt = perAsc.begin(); 
	       perIt != perAsc.end(); ++perIt) {
	    int mappedDof = perMap.map(feSpace, *perIt, globalRowDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
891
	    int mappedIndex = interiorMap->getMatIndex(rowComp, mappedDof);
892

Thomas Witkowski's avatar
Thomas Witkowski committed
893
	    VecSetValue(vecInterior, mappedIndex, value, ADD_VALUES);
894 895 896
	  }	  
	} else {	  
	  // The DOF index is not periodic.
Thomas Witkowski's avatar
Thomas Witkowski committed
897 898
	  VecSetValue(vecInterior, index, *dofIt, ADD_VALUES);
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
899 900 901 902
      }
    }
  }

903 904 905 906 907 908 909 910

  PetscSolver* PetscSolverGlobalMatrix::createSubSolver(int component,
							string kspPrefix)
  {
    vector<const FiniteElemSpace*> fe;
    fe.push_back(componentSpaces[component]);

    PetscSolver* subSolver = new PetscSolverGlobalMatrix("");
911
    subSolver->setKspPrefix(kspPrefix);
912
    subSolver->setMeshDistributor(meshDistributor, 0);
913 914 915 916 917 918 919 920 921
    subSolver->init(fe, fe);

    ParallelDofMapping &subDofMap = subSolver->getDofMap();
    subDofMap[0] = dofMap[component];
    subDofMap.update();

    return subSolver;
  }

922 923 924 925

  void PetscSolverGlobalMatrix::setConstantNullSpace(KSP ksp,
						     int constFeSpace,
						     bool test)
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942
  {
    FUNCNAME("PetscSolverGlobalMatrix::setConstantNullSpace()");

    Vec nullSpaceBasis;
    VecDuplicate(getVecSolInterior(), &nullSpaceBasis);
    
    SystemVector basisVec("tmp", componentSpaces, componentSpaces.size(), true);
    basisVec.set(0.0);
    basisVec.getDOFVector(constFeSpace)->set(1.0);
    
    setDofVector(nullSpaceBasis, basisVec, true);
    VecAssemblyBegin(nullSpaceBasis);
    VecAssemblyEnd(nullSpaceBasis);
    VecNormalize(nullSpaceBasis, PETSC_NULL);
    
    if (test) {
      Vec tmp;
943 944 945
#if (PETSC_VERSION_MINOR >= 6)
      MatCreateVecs(getMatInterior(), &tmp, PETSC_NULL);
#else
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
946
      MatGetVecs(getMatInterior(), &tmp, PETSC_NULL);
947
#endif
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
948 949 950 951 952
      MatMult(getMatInterior(), nullSpaceBasis, tmp);
      PetscReal n;
      VecNorm(tmp, NORM_2, &n);
      MSG("NORM IS: %e\n", n);
      VecDestroy(&tmp);