PetscSolverFeti.cc 77.7 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12
//
// 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.


13
#include "AMDiS.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
14
#include "MatrixVector.h"
15
#include "parallel/PetscSolverFeti.h"
16
#include "parallel/PetscSolverFetiStructs.h"
17 18
#include "parallel/PetscSolverFetiOperators.h"
#include "parallel/PetscSolverFetiTimings.h"
19 20
#include "parallel/StdMpi.h"
#include "parallel/MpiHelper.h"
21
#include "parallel/PetscSolverGlobalMatrix.h"
22
#include "io/VtkWriter.h"
23 24 25 26 27

namespace AMDiS {

  using namespace std;

28
  PetscErrorCode KSPMonitorFeti(KSP ksp, PetscInt n, PetscReal rnorm, void *data)
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
29
  {    
30 31 32
    Vec Br,v,w;
    VecDuplicate(static_cast<FetiKspData*>(data)->draft, &v);
    VecDuplicate(static_cast<FetiKspData*>(data)->draft, &w);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
33 34 35 36 37 38 39

    KSPBuildResidual(ksp, v, w, &Br);

    Vec nest0, nest1;
    VecNestGetSubVec(Br, 0, &nest0);
    VecNestGetSubVec(Br, 1, &nest1);

40 41
    PetscScalar norm, norm0, norm1;
    VecNorm(Br, NORM_2, &norm);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
42 43 44 45 46 47
    VecNorm(nest0, NORM_2, &norm0);
    VecNorm(nest1, NORM_2, &norm1);

    VecDestroy(&v);
    VecDestroy(&w);

48 49
    PetscPrintf(PETSC_COMM_WORLD, "%3D KSP residual norm %1.12e [%1.12e %1.12e] and preconditioned norm [%1.12e]\n",
		n, norm, norm0, norm1, rnorm);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
50 51 52 53

    return 0;
  }

54 55
  PetscSolverFeti::PetscSolverFeti()
    : PetscSolver(),
56
      schurPrimalSolver(0),
57
      multiLevelTest(false),
Thomas Witkowski's avatar
Thomas Witkowski committed
58
      subdomain(NULL),
59
      massMatrixSolver(NULL),
60 61
      meshLevel(0),
      rStartInterior(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
62
      nGlobalOverallInterior(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
63
      printTimings(false),
Thomas Witkowski's avatar
Thomas Witkowski committed
64 65
      stokesMode(false),
      pressureComponent(-1)
66 67 68 69 70 71
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

    string preconditionerName = "";
    Parameters::get("parallel->solver->precon", preconditionerName);
    if (preconditionerName == "" || preconditionerName == "none") {
72
      MSG("Create FETI-DP solver with no preconditioner!\n");
73 74
      fetiPreconditioner = FETI_NONE;
    } else if (preconditionerName == "dirichlet") {
75
      MSG("Create FETI-DP solver with Dirichlet preconditioner!\n");
76 77
      fetiPreconditioner = FETI_DIRICHLET;
    } else if (preconditionerName == "lumped") {
78
      MSG("Create FETI-DP solver with lumped preconditioner!\n");
79 80
      fetiPreconditioner = FETI_LUMPED;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
81 82
      ERROR_EXIT("Preconditioner \"%s\" not available!\n", 
		 preconditionerName.c_str());
83
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
84 85 86 87 88

    Parameters::get("parallel->feti->schur primal solver", schurPrimalSolver);
    TEST_EXIT(schurPrimalSolver == 0 || schurPrimalSolver == 1)
      ("Wrong solver \"%d\"for the Schur primal complement!\n", 
       schurPrimalSolver);
89 90

    Parameters::get("parallel->multi level test", multiLevelTest);
91 92
    if (multiLevelTest)
      meshLevel = 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
93 94

    Parameters::get("parallel->print timings", printTimings);
95 96 97
  }


98
  void PetscSolverFeti::initialize(vector<const FiniteElemSpace*> feSpaces)
99
  {
100 101
    FUNCNAME("PetscSolverFeti::initialize()");

102 103 104
    TEST_EXIT_DBG(meshLevel + 1 == meshDistributor->getMeshLevelData().getLevelNumber())
      ("Mesh hierarchy does not contain %d levels!\n", meshLevel + 1);

105
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();
106
    vector<const FiniteElemSpace*>& uniqueFe = meshDistributor->getFeSpaces();
107

Thomas Witkowski's avatar
Thomas Witkowski committed
108

Thomas Witkowski's avatar
Thomas Witkowski committed
109 110 111 112 113 114 115 116
    stokesMode = false;
    Parameters::get("parallel->feti->stokes mode", stokesMode);
    if (stokesMode) {
      Parameters::get("parallel->feti->pressure component", pressureComponent);
      TEST_EXIT(pressureComponent >= 0)
	("FETI-DP in Stokes mode, no pressure component defined!\n");

      pressureFeSpace = feSpaces[pressureComponent];
Thomas Witkowski's avatar
Thomas Witkowski committed
117 118
    }
			   
Thomas Witkowski's avatar
Thomas Witkowski committed
119 120
    if (subdomain == NULL) {
      subdomain = new PetscSolverGlobalMatrix();
121

122
      if (meshLevel == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
123
	subdomain->setMeshDistributor(meshDistributor, 
124
				      mpiCommGlobal, mpiCommLocal);
125
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
126
	subdomain->setMeshDistributor(meshDistributor, 
127 128
				      levelData.getMpiComm(meshLevel - 1),
				      levelData.getMpiComm(meshLevel));
Thomas Witkowski's avatar
Thomas Witkowski committed
129
	subdomain->setLevel(meshLevel);
130 131
      }
    }
132

133 134 135 136
    primalDofMap.init(levelData, feSpaces, uniqueFe);   
    dualDofMap.init(levelData, feSpaces, uniqueFe, false);
    localDofMap.init(levelData, feSpaces, uniqueFe, meshLevel != 0);
    lagrangeMap.init(levelData, feSpaces, uniqueFe);
Thomas Witkowski's avatar
Thomas Witkowski committed
137

Thomas Witkowski's avatar
Thomas Witkowski committed
138
    if (stokesMode)
Thomas Witkowski's avatar
Thomas Witkowski committed
139
      interfaceDofMap.init(levelData, feSpaces, uniqueFe);
Thomas Witkowski's avatar
Thomas Witkowski committed
140

141
    if (fetiPreconditioner == FETI_DIRICHLET) {
142 143 144
      TEST_EXIT(meshLevel == 0)
	("Dirichlet preconditioner not yet implemented for multilevel FETI-DP\n");

145
      interiorDofMap.init(levelData, feSpaces, uniqueFe, false);
146
    }
147 148 149
  }


Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
150 151 152 153
  void PetscSolverFeti::createDirichletData(Matrix<DOFMatrix*> &mat)
  {
    FUNCNAME("PetscSolverFeti::createDirichletData()");

Thomas Witkowski's avatar
Thomas Witkowski committed
154 155 156 157 158
    bool removeDirichletRows = false;
    Parameters::get("parallel->feti->remove dirichlet", removeDirichletRows);
    if (!removeDirichletRows)
      return;

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175
    int nComponents = mat.getSize();
    for (int i = 0; i < nComponents; i++) {
      DOFMatrix* dofMat = mat[i][i];
      if (!dofMat)
	continue;
      
      const FiniteElemSpace *feSpace = dofMat->getRowFeSpace();
      std::set<DegreeOfFreedom>& dRows = dofMat->getDirichletRows();
      if (dirichletRows.count(feSpace)) {
	WARNING("Implement test that for all components dirichlet rows are equal!\n");
      } else {
	dirichletRows[feSpace] = dRows;
      }
    }
  }


176 177 178 179
  void PetscSolverFeti::createFetiData()
  {
    FUNCNAME("PetscSolverFeti::createFetiData()");

Thomas Witkowski's avatar
Thomas Witkowski committed
180 181
    double timeCounter = MPI::Wtime();

182 183 184
    TEST_EXIT(meshDistributor)("No mesh distributor object defined!\n");
    TEST_EXIT(meshDistributor->getFeSpaces().size() > 0)
      ("No FE space defined in mesh distributor!\n");
185

186 187
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

188 189 190 191
    primalDofMap.clear();
    dualDofMap.clear();
    lagrangeMap.clear();
    localDofMap.clear();
192
    if (fetiPreconditioner == FETI_DIRICHLET)
193 194
      interiorDofMap.clear();

195 196
    primalDofMap.setDofComm(meshDistributor->getDofComm());
    lagrangeMap.setDofComm(meshDistributor->getDofComm());
197

198 199
    primalDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    dualDofMap.setMpiComm(levelData.getMpiComm(0), 0);
200
    lagrangeMap.setMpiComm(levelData.getMpiComm(0), 0);
201
    localDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
202
    if (fetiPreconditioner == FETI_DIRICHLET)
203
      interiorDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
204

205 206 207
    if (meshLevel == 0)
      localDofMap.setDofComm(meshDistributor->getDofComm());
    else
208 209
      localDofMap.setDofComm(meshDistributor->getDofCommSd());

Thomas Witkowski's avatar
Thomas Witkowski committed
210
    if (stokesMode) {
211 212 213 214 215
      interfaceDofMap.clear();
      interfaceDofMap.setDofComm(meshDistributor->getDofComm());
      interfaceDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    }

216 217
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
218
    
219 220
      createPrimals(feSpace);  

221 222 223
      createDuals(feSpace);

      createInterfaceNodes(feSpace);
224

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
225
      createIndexB(feSpace);     
226
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
227 228 229 230

    primalDofMap.update();
    dualDofMap.update();
    localDofMap.update();
231
    if (fetiPreconditioner == FETI_DIRICHLET)
232
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
233

Thomas Witkowski's avatar
Thomas Witkowski committed
234
    if (stokesMode)
235 236
      interfaceDofMap.update();

237 238 239 240
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
      createLagrange(feSpace);
    }
241

242 243 244
    lagrangeMap.update();


245 246 247 248 249 250 251 252 253 254 255 256
    // === ===

    if (meshLevel == 0) {
      rStartInterior = 0;
      nGlobalOverallInterior = localDofMap.getOverallDofs();
    } else {
      MeshLevelData& levelData = meshDistributor->getMeshLevelData();

      int groupRowsInterior = 0;
      if (levelData.getMpiComm(1).Get_rank() == 0)
	groupRowsInterior = localDofMap.getOverallDofs();

257
      mpi::getDofNumbering(mpiCommGlobal, groupRowsInterior,
258 259 260 261 262 263 264 265 266
			   rStartInterior, nGlobalOverallInterior);

      int tmp = 0;
      if (levelData.getMpiComm(1).Get_rank() == 0)
	tmp = rStartInterior;

      levelData.getMpiComm(1).Allreduce(&tmp, &rStartInterior, 1, MPI_INT, MPI_SUM);
    }

267
    MSG("FETI-DP data created on mesh level %d\n", meshLevel);
Thomas Witkowski's avatar
Thomas Witkowski committed
268 269
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
Thomas Witkowski's avatar
Thomas Witkowski committed
270
      MSG("FETI-DP data for %d-ith FE space: %p\n", i, feSpace);
Thomas Witkowski's avatar
Thomas Witkowski committed
271

Thomas Witkowski's avatar
Thomas Witkowski committed
272
      if (feSpace == pressureFeSpace) {
Thomas Witkowski's avatar
Thomas Witkowski committed
273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288
	MSG("  nRankInterface = %d  nOverallInterface = %d\n",
	    interfaceDofMap[feSpace].nRankDofs, 
	    interfaceDofMap[feSpace].nOverallDofs);
      } else {
	MSG("  nRankPrimals = %d   nOverallPrimals = %d\n", 
	    primalDofMap[feSpace].nRankDofs, 
	    primalDofMap[feSpace].nOverallDofs);
	
	MSG("  nRankDuals = %d  nOverallDuals = %d\n",
	    dualDofMap[feSpace].nRankDofs, 
	    dualDofMap[feSpace].nOverallDofs);
	
	MSG("  nRankLagrange = %d  nOverallLagrange = %d\n",
	    lagrangeMap[feSpace].nRankDofs, 
	    lagrangeMap[feSpace].nOverallDofs);

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
289 290 291
// 	TEST_EXIT_DBG(localDofMap[feSpace].size() + primalDofMap[feSpace].size() == 
// 		      static_cast<unsigned int>(feSpace->getAdmin()->getUsedDofs()))
// 	  ("Should not happen!\n");	
Thomas Witkowski's avatar
Thomas Witkowski committed
292
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
293
    }
294

Thomas Witkowski's avatar
Thomas Witkowski committed
295
    subdomain->setDofMapping(&localDofMap);
Thomas Witkowski's avatar
Thomas Witkowski committed
296 297 298
    subdomain->setCoarseSpaceDofMapping(&primalDofMap); 
    if (stokesMode)
      subdomain->setCoarseSpaceDofMapping(&interfaceDofMap, pressureComponent);
Thomas Witkowski's avatar
Thomas Witkowski committed
299 300

    if (printTimings) {
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
301
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
302
      timeCounter = MPI::Wtime() - timeCounter;
303
      MSG("FETI-DP timing 01: %.5f seconds (creation of basic data structures)\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
304 305
	  timeCounter);
    }
306 307 308
  }


309
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
310
  {
311
    FUNCNAME("PetscSolverFeti::createPrimals()");  
312

Thomas Witkowski's avatar
Thomas Witkowski committed
313
    if (feSpace == pressureFeSpace)
314 315
      return;

316 317 318
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

Thomas Witkowski's avatar
Thomas Witkowski committed
319
    // Set of DOF indices that are considered to be primal variables.
320
    DofContainerSet& vertices = 
321
      meshDistributor->getBoundaryDofInfo(feSpace, meshLevel).geoDofs[VERTEX];
322 323

    DofIndexSet primals;
324
    for (DofContainerSet::iterator it = vertices.begin(); 
325
	 it != vertices.end(); ++it) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
326 327
      if (dirichletRows[feSpace].count(**it))
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
328

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
329 330
      if (meshLevel == 0) {
	primals.insert(**it);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
331
      } else {
332
	double e = 1e-8;
Thomas Witkowski's avatar
Thomas Witkowski committed
333 334 335
	WorldVector<double> c;
	feSpace->getMesh()->getDofIndexCoords(*it, feSpace, c);

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
336 337 338 339 340 341 342
	if (fabs(c[0]) < e || fabs(c[1]) < e ||
	    fabs(c[0] - 25.0) < e || fabs(c[1] - 25.0) < e ||
	    (fabs(c[0] - 12.5) < e && fabs(c[1] - 12.5) < e)) {
	  MSG("PRIMAL COORD %f %f\n", c[0], c[1]);
	  primals.insert(**it);
	}
      }
343
    }
344

345 346 347 348

    // === Calculate the number of primals that are owned by the rank and ===
    // === create local indices of the primals starting at zero.          ===

349
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
350 351
      if (meshDistributor->getDofMap()[feSpace].isRankDof(*it))
	primalDofMap[feSpace].insertRankDof(*it);
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
352
      else
353
  	primalDofMap[feSpace].insertNonRankDof(*it);
354 355 356
  }


357
  void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace)
358 359
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
360

Thomas Witkowski's avatar
Thomas Witkowski committed
361
    if (feSpace == pressureFeSpace)
362 363
      return;

364 365 366
    // === Create global index of the dual nodes on each rank. ===

    DofContainer allBoundaryDofs;
367
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
368
    
369
    for (DofContainer::iterator it = allBoundaryDofs.begin();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
370 371 372 373 374 375 376 377 378 379 380
	 it != allBoundaryDofs.end(); ++it) {
      if (dirichletRows[feSpace].count(**it))
	continue;

      if (isPrimal(feSpace, **it))
	continue;

      if (meshLevel == 0) {
	dualDofMap[feSpace].insertRankDof(**it);
      } else {
	if (meshDistributor->getDofMapSd()[feSpace].isRankDof(**it))
381
	  dualDofMap[feSpace].insertRankDof(**it);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
382 383
      }	  
    }
384 385 386
  }

  
387 388 389 390
  void PetscSolverFeti::createInterfaceNodes(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createInterfaceNodes()");

Thomas Witkowski's avatar
Thomas Witkowski committed
391
    if (feSpace != pressureFeSpace)
392 393 394 395 396 397
      return;

    DofContainer allBoundaryDofs;
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);

    for (DofContainer::iterator it = allBoundaryDofs.begin();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
398 399 400 401
	 it != allBoundaryDofs.end(); ++it) {
      if (dirichletRows[feSpace].count(**it))
	continue;      
      
402
      if (meshDistributor->getDofMap()[feSpace].isRankDof(**it))
Thomas Witkowski's avatar
Thomas Witkowski committed
403
	interfaceDofMap[feSpace].insertRankDof(**it);
404
      else
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
405 406
	interfaceDofMap[feSpace].insertNonRankDof(**it);      
    }
407 408 409
  }


410 411 412 413
  void PetscSolverFeti::createLagrange(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
414
    if (feSpace == pressureFeSpace)
415 416
      return;

417 418
    boundaryDofRanks[feSpace].clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
419 420 421
    // Stores for all rank owned communication DOFs, if the counterpart is
    // a rank owned DOF in its subdomain. Thus, the following map stores to
    // each rank number all DOFs that fulfill this requirenment.
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
422
    map<int, std::set<DegreeOfFreedom> > sdRankDofs;
423

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
424
    if (meshLevel > 0) {
425
      StdMpi<vector<int> > stdMpi(mpiCommGlobal);
426 427 428 429 430 431 432 433 434

      for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), 
				meshLevel, feSpace);
	   !it.end(); it.nextRank()) {

	vector<int> subdomainRankDofs;
	subdomainRankDofs.reserve(it.getDofs().size());

	for (; !it.endDofIter(); it.nextDof()) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
435
	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))
436
	    subdomainRankDofs.push_back(1);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
437
	  else
438 439 440 441 442 443 444 445 446 447 448 449 450 451 452
	    subdomainRankDofs.push_back(0);
	}

	stdMpi.send(it.getRank(), subdomainRankDofs);
      }	     

      for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), 
				meshLevel, feSpace);
	   !it.end(); it.nextRank())
	stdMpi.recv(it.getRank());

      stdMpi.startCommunication();

      for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), 
				meshLevel, feSpace); 
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
453 454 455 456 457 458
	   !it.end(); it.nextRank())
	for (; !it.endDofIter(); it.nextDof())
	  if (!isPrimal(feSpace, it.getDofIndex()))
	    if (stdMpi.getRecvData(it.getRank())[it.getDofCounter()] == 1)
	      sdRankDofs[it.getRank()].insert(it.getDofIndex());
    }
459

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
460 461 462 463 464 465 466
    if (dualDofMap[feSpace].nLocalDofs == 0)
      return;


    // === Create for each dual node that is owned by the rank, the set ===
    // === of ranks that contain this node (denoted by W(x_j)).         ===

467
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
468 469 470 471 472 473 474
    for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), 
			      meshLevel, feSpace); 
	 !it.end(); it.nextRank()) {
      for (; !it.endDofIter(); it.nextDof()) {
	if (!isPrimal(feSpace, it.getDofIndex())) {
	  boundaryDofRanks[feSpace][it.getDofIndex()].insert(mpiRank);

475
 	  if (meshLevel == 0 ||
Thomas Witkowski's avatar
Thomas Witkowski committed
476
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
477
	    boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());	  
478 479 480 481
	}
      }
    }

482 483 484 485

    // === Communicate these sets for all rank owned dual nodes to other ===
    // === ranks that also have this node.                               ===

486
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
487

488
    for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), meshLevel, feSpace);
489 490
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
491
	if (!isPrimal(feSpace, it.getDofIndex()))
492 493 494
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
	    stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]);
495 496 497

    stdMpi.updateSendDataSize();

498
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
499
	 !it.end(); it.nextRank()) {
500
      bool recvFromRank = false;
501
      for (; !it.endDofIter(); it.nextDof()) {
502
	if (!isPrimal(feSpace, it.getDofIndex())) {
503 504 505
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
 	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
506 507 508
	    recvFromRank = true;
	    break;
	  }
509
	}
510
      }
511 512

      if (recvFromRank)
513
	stdMpi.recv(it.getRank());
514
    }
515

516 517
    stdMpi.startCommunication();

518
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
519
	 !it.end(); it.nextRank()) {
520
      int i = 0;
521
      for (; !it.endDofIter(); it.nextDof())
522
	if (!isPrimal(feSpace, it.getDofIndex()))
523 524 525
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
 	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex())))	    
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
526 527
	    boundaryDofRanks[feSpace][it.getDofIndex()] = 
	      stdMpi.getRecvData(it.getRank())[i++];
528 529
	  else
	    lagrangeMap[feSpace].insertNonRankDof(it.getDofIndex());
530 531
    }

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
532

533 534 535
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

536
    int nRankLagrange = 0;
537
    DofMap& dualMap = dualDofMap[feSpace].getMap();
538
    for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
539 540
      if (meshDistributor->getDofMap()[feSpace].isRankDof(it->first)) {
	lagrangeMap[feSpace].insertRankDof(it->first, nRankLagrange);
541
	int degree = boundaryDofRanks[feSpace][it->first].size();
542
	nRankLagrange += (degree * (degree - 1)) / 2;
Thomas Witkowski's avatar
Thomas Witkowski committed
543
      } else {
544
	lagrangeMap[feSpace].insertNonRankDof(it->first);
545 546
      }
    }
547
    lagrangeMap[feSpace].nRankDofs = nRankLagrange;
548 549 550
  }


551
  void PetscSolverFeti::createIndexB(const FiniteElemSpace *feSpace)
552
  {
553
    FUNCNAME("PetscSolverFeti::createIndexB()");
554

555
    DOFAdmin* admin = feSpace->getAdmin();
556 557 558 559

    // === To ensure that all interior node on each rank are listen first in ===
    // === the global index of all B nodes, insert all interior nodes first, ===
    // === without defining a correct index.                                 ===
560

561
    int nLocalInterior = 0;    
562
    for (int i = 0; i < admin->getUsedSize(); i++) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
563 564 565 566 567 568
      if (admin->isDofFree(i) ||
	  isPrimal(feSpace, i) ||
	  isDual(feSpace, i) ||
	  isInterface(feSpace, i) ||
	  dirichletRows[feSpace].count(i))
	continue;      
569

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
570 571 572 573 574 575 576 577 578 579 580 581 582 583 584
      if (meshLevel == 0) {
	localDofMap[feSpace].insertRankDof(i, nLocalInterior);
	
	if (fetiPreconditioner == FETI_DIRICHLET)
	  interiorDofMap[feSpace].insertRankDof(i, nLocalInterior);
	
	nLocalInterior++;	
      } else {
	if (meshDistributor->getDofMapSd()[feSpace].isRankDof(i))
	  localDofMap[feSpace].insertRankDof(i);
	else
	  localDofMap[feSpace].insertNonRankDof(i);
	
	TEST_EXIT_DBG(fetiPreconditioner == FETI_NONE)
	  ("Not yet implemnted!\n");	
585
      }
586
    }
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
587
    
588 589
    // === And finally, add the global indicies of all dual nodes. ===

590
    for (DofMap::iterator it = dualDofMap[feSpace].getMap().begin();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
591 592
	 it != dualDofMap[feSpace].getMap().end(); ++it) {
      if (meshLevel == 0) {
593
	localDofMap[feSpace].insertRankDof(it->first);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
594
      } else {
595 596 597 598 599
	if (meshDistributor->getDofMapSd()[feSpace].isRankDof(it->first))
	  localDofMap[feSpace].insertRankDof(it->first);
	else 
	  localDofMap[feSpace].insertNonRankDof(it->first);
      }
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
600
    }
601 602 603
  }


604
  void PetscSolverFeti::createMatLagrange(vector<const FiniteElemSpace*> &feSpaces)
605 606 607
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
608
    double wtime = MPI::Wtime();
609
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Thomas Witkowski committed
610

611 612
    // === Create distributed matrix for Lagrange constraints. ===

613 614 615 616 617
    MatCreateAIJ(mpiCommGlobal,
		 lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		 lagrangeMap.getOverallDofs(), nGlobalOverallInterior,
		 2, PETSC_NULL, 2, PETSC_NULL,
		 &mat_lagrange);
618
    MatSetOption(mat_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
619

620 621 622 623 624 625 626
    // === Create for all duals the corresponding Lagrange constraints. On ===
    // === each rank we traverse all pairs (n, m) of ranks, with n < m,    ===
    // === that contain this node. If the current rank number is r, and    ===
    // === n == r, the rank sets 1.0 for the corresponding constraint, if  ===
    // === m == r, than the rank sets -1.0 for the corresponding           ===
    // === constraint.                                                     ===

Thomas Witkowski's avatar
Thomas Witkowski committed
627
    for (unsigned int k = 0; k < feSpaces.size(); k++) {
628
      DofMap &dualMap = dualDofMap[feSpaces[k]].getMap();
629

630
      for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
631
	TEST_EXIT_DBG(boundaryDofRanks[feSpaces[k]].count(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
632 633 634
	  ("Should not happen!\n");
	
	// Global index of the first Lagrange constriant for this node.
635
	int index = lagrangeMap.getMatIndex(k, it->first);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
636

Thomas Witkowski's avatar
Thomas Witkowski committed
637
	// Copy set of all ranks that contain this dual node.
638 639
	vector<int> W(boundaryDofRanks[feSpaces[k]][it->first].begin(), 
		      boundaryDofRanks[feSpaces[k]][it->first].end());
Thomas Witkowski's avatar
Thomas Witkowski committed
640 641
	// Number of ranks that contain this dual node.
	int degree = W.size();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
642 643

	TEST_EXIT_DBG(degree > 1)("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
644 645 646 647
	
	for (int i = 0; i < degree; i++) {
	  for (int j = i + 1; j < degree; j++) {
	    if (W[i] == mpiRank || W[j] == mpiRank) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
648 649
	      int colIndex = 
		localDofMap.getMatIndex(k, it->first) + rStartInterior;
650

651
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
652
	      MatSetValue(mat_lagrange, index, colIndex, value, INSERT_VALUES);
653
	    }
Thomas Witkowski's avatar
Thomas Witkowski committed
654
	    index++;	      
655 656 657 658 659 660 661
	  }
	}
      }
    }

    MatAssemblyBegin(mat_lagrange, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(mat_lagrange, MAT_FINAL_ASSEMBLY);
Thomas Witkowski's avatar
Thomas Witkowski committed
662

Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
663 664
    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
665 666
      MSG("FETI-DP timing 05: %.5f seconds (creation of lagrange constraint matrix)\n", 
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
667
    }
668 669 670
  }


671
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
672
  {
673
    FUNCNAME("PetscSolverFeti::createSchurPrimalKsp()");
674

Thomas Witkowski's avatar
Thomas Witkowski committed
675
    if (schurPrimalSolver == 0) {
676 677
      MSG("Create iterative schur primal solver!\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
678
      schurPrimalData.subSolver = subdomain;
Thomas Witkowski's avatar
Thomas Witkowski committed
679
      
680 681
      localDofMap.createVec(schurPrimalData.tmp_vec_b, nGlobalOverallInterior);
      primalDofMap.createVec(schurPrimalData.tmp_vec_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
682

683
      MatCreateShell(mpiCommGlobal,
684 685 686 687
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getOverallDofs(), 
		     primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
688 689 690 691 692
		     &schurPrimalData, 
		     &mat_schur_primal);
      MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			   (void(*)(void))petscMultMatSchurPrimal);
      
693
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
694
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
695 696
      KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
      KSPSetType(ksp_schur_primal, KSPGMRES);
Thomas Witkowski's avatar
Thomas Witkowski committed
697 698
      KSPSetFromOptions(ksp_schur_primal);
    } else {
699 700
      MSG("Create direct schur primal solver!\n");

701 702
      double wtime = MPI::Wtime();

703 704 705
      TEST_EXIT_DBG(meshLevel == 0)
	("Does not support for multilevel, check usage of localDofMap.\n");

706 707 708
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
      int nRowsRankB = localDofMap.getRankDofs();
709

Thomas Witkowski's avatar
Thomas Witkowski committed
710
      Mat matBPi;
711 712 713
      MatCreateAIJ(mpiCommGlobal,
		   nRowsRankB, nRowsRankPrimal, 
		   nGlobalOverallInterior, nRowsOverallPrimal,
714 715 716
		   150, PETSC_NULL, 150, PETSC_NULL, &matBPi);
      MatSetOption(matBPi, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

717

718 719 720
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
721

722 723
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

724
      for (int i = 0; i < nRowsRankB; i++) {
725
	PetscInt row = localDofMap.getStartDofs() + i;
726
	MatGetRow(subdomain->getMatInteriorCoarse(), row, &nCols, &cols, &values);
727 728 729 730 731

	for (int j = 0; j < nCols; j++)
	  if (values[j] != 0.0)
	    mat_b_primal_cols[cols[j]].push_back(make_pair(i, values[j]));

732
	MatRestoreRow(subdomain->getMatInteriorCoarse(), row, &nCols, &cols, &values);
733 734
      }

735
      TEST_EXIT(static_cast<int>(mat_b_primal_cols.size()) == 
736
		primalDofMap.getLocalDofs())
737
	("Should not happen!\n");
738

739 740 741
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
742
	VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);
743 744 745 746 747 748 749

 	for (unsigned int i = 0; i < it->second.size(); i++) 
 	  VecSetValue(tmpVec, 
 		      it->second[i].first, it->second[i].second, INSERT_VALUES);

	VecAssemblyBegin(tmpVec);
	VecAssemblyEnd(tmpVec);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
750

Thomas Witkowski's avatar
Thomas Witkowski committed
751
	subdomain->solve(tmpVec, tmpVec);
752

Thomas Witkowski's avatar
Thomas Witkowski committed
753
	PetscScalar *tmpValues;
754
	VecGetArray(tmpVec, &tmpValues);
755
	for (int i  = 0; i < nRowsRankB; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
756
	  MatSetValue(matBPi, 
757
		      localDofMap.getStartDofs() + i,
Thomas Witkowski's avatar
Thomas Witkowski committed
758 759 760
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
761 762 763 764 765
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
766 767
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
768

769 770
      MatDuplicate(subdomain->getMatCoarse(), MAT_COPY_VALUES, 
		   &mat_schur_primal);
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
771 772

      Mat matPrimal;
773 774
      MatMatMult(subdomain->getMatCoarseInterior(), matBPi, MAT_INITIAL_MATRIX, 
		 PETSC_DEFAULT, &matPrimal);
775
      MatAXPY(mat_schur_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);
776 777

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
778
      MatDestroy(&matBPi);
779

780
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
781
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal,
782
		      SAME_NONZERO_PATTERN);
783 784 785 786 787 788
      KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
      KSPSetType(ksp_schur_primal, KSPPREONLY);
      PC pc_schur_primal;      
      KSPGetPC(ksp_schur_primal, &pc_schur_primal);
      PCSetType(pc_schur_primal, PCLU);
      PCFactorSetMatSolverPackage(pc_schur_primal, MATSOLVERMUMPS);
789
      KSPSetFromOptions(ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
790

Thomas Witkowski's avatar
Thomas Witkowski committed
791
      if (printTimings) {
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
792
	MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
793 794 795 796 797 798 799 800 801 802
	MatInfo minfo;
	MatGetInfo(mat_schur_primal, MAT_GLOBAL_SUM, &minfo);
	MSG("Schur primal matrix nnz = %f\n", minfo.nz_used);
	
	MSG("FETI-DP timing 06: %.5f seconds (creation of schur primal matrix)\n",
	    MPI::Wtime() - wtime);

	wtime = MPI::Wtime();
	KSPSetUp(ksp_schur_primal);
	KSPSetUpOnBlocks(ksp_schur_primal);
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
803
	MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
804 805 806
	MSG("FETI-DP timing 07: %.5f seconds (factorization of primal schur matrix).\n",
	    MPI::Wtime() - wtime);
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
807
    }
808 809 810 811 812 813 814
  }


  void PetscSolverFeti::destroySchurPrimalKsp()
  {
    FUNCNAME("PetscSolverFeti::destroySchurPrimal()");

815
    if (schurPrimalSolver == 0) {
816
      schurPrimalData.subSolver = NULL;
817

818 819 820
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
    }
821 822 823
    
    MatDestroy(&mat_schur_primal);
    KSPDestroy(&ksp_schur_primal);
824 825 826
  }


827
  void PetscSolverFeti::createFetiKsp(vector<const FiniteElemSpace*> &feSpaces)
828 829 830
  {
    FUNCNAME("PetscSolverFeti::createFetiKsp()");

831 832
    // === Create FETI-DP solver object. ===

833 834 835 836
    fetiData.mat_lagrange = &mat_lagrange;
    fetiData.subSolver = subdomain;
    fetiData.ksp_schur_primal = &ksp_schur_primal;

Thomas Witkowski's avatar
Thomas Witkowski committed
837
    localDofMap.createVec(fetiData.tmp_vec_b0, nGlobalOverallInterior);
838
    lagrangeMap.createVec(fetiData.tmp_vec_lagrange);
Thomas Witkowski's avatar
Thomas Witkowski committed
839
    primalDofMap.createVec(fetiData.tmp_vec_primal0);
840

Thomas Witkowski's avatar
Thomas Witkowski committed
841
    if (stokesMode == false) {      
842 843 844 845 846 847 848 849
      MatCreateShell(mpiCommGlobal,
		     lagrangeMap.getRankDofs(), 
		     lagrangeMap.getRankDofs(),
		     lagrangeMap.getOverallDofs(), 
		     lagrangeMap.getOverallDofs(),
		     &fetiData, &mat_feti);
      MatShellSetOperation(mat_feti, MATOP_MULT, 
			   (void(*)(void))petscMultMatFeti);
850
    }  else {
Thomas Witkowski's avatar
Thomas Witkowski committed
851 852 853 854