PetscSolverFeti.cc 76.6 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;

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53
  PetscErrorCode KSPMonitorFeti(KSP ksp, PetscInt n, PetscReal rnorm, void *dummy)
  {    
    Vec            Br,v,w;
    Mat            A;

    KSPGetOperators(ksp, &A, PETSC_NULL, PETSC_NULL);
    MatGetVecs(A, &w, &v);
    KSPBuildResidual(ksp, v, w, &Br);

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

    PetscScalar norm0, norm1;
    VecNorm(nest0, NORM_2, &norm0);
    VecNorm(nest1, NORM_2, &norm1);

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

    PetscPrintf(PETSC_COMM_WORLD, "%3D KSP Component residual norm [%1.12e %1.12e] and relative [%1.12e]\n",
		norm0, norm1, rnorm);

    return 0;
  }

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

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

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

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

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


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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
107

Thomas Witkowski's avatar
Thomas Witkowski committed
108 109 110 111 112 113 114 115
    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
116 117
    }
			   
Thomas Witkowski's avatar
Thomas Witkowski committed
118 119
    if (subdomain == NULL) {
      subdomain = new PetscSolverGlobalMatrix();
120

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

132 133 134 135
    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
136

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

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

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


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

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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174
    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;
      }
    }
  }


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

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

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

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

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

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

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

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

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

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

220 221 222
      createDuals(feSpace);

      createInterfaceNodes(feSpace);
223

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

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

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

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

241 242 243
    lagrangeMap.update();


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

    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();

256
      mpi::getDofNumbering(mpiCommGlobal, groupRowsInterior,
257 258 259 260 261 262 263 264 265
			   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);
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
271
      if (feSpace == pressureFeSpace) {
Thomas Witkowski's avatar
Thomas Witkowski committed
272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287
	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
288 289 290
// 	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
291
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
292
    }
293

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

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


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

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

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

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

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

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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
335 336 337 338 339 340 341
	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);
	}
      }
342
    }
343

344 345 346 347

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

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


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

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

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

    DofContainer allBoundaryDofs;
366
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
367
    
368
    for (DofContainer::iterator it = allBoundaryDofs.begin();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
369 370 371 372 373 374 375 376 377 378 379
	 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))
380
	  dualDofMap[feSpace].insertRankDof(**it);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
381 382
      }	  
    }
383 384 385
  }

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

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

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

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


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

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

416 417
    boundaryDofRanks[feSpace].clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
418 419 420
    // 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
421
    map<int, std::set<DegreeOfFreedom> > sdRankDofs;
422

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

      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
434
	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))
435
	    subdomainRankDofs.push_back(1);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
436
	  else
437 438 439 440 441 442 443 444 445 446 447 448 449 450 451
	    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
452 453 454 455 456 457
	   !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());
    }
458

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
459 460 461 462 463 464 465
    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)).         ===

466
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
467 468 469 470 471 472 473
    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);

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

481 482 483 484

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

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

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

    stdMpi.updateSendDataSize();

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

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

515 516
    stdMpi.startCommunication();

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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
531

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

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


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

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

    // === 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.                                 ===
559

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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
569 570 571 572 573 574 575 576 577 578 579 580 581 582 583
      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");	
584
      }
585
    }
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
586
    
587 588
    // === And finally, add the global indicies of all dual nodes. ===

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


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

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

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

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

619 620 621 622 623 624 625
    // === 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
626
    for (unsigned int k = 0; k < feSpaces.size(); k++) {
627
      DofMap &dualMap = dualDofMap[feSpaces[k]].getMap();
628

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

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

	TEST_EXIT_DBG(degree > 1)("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
643 644 645 646
	
	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
647 648
	      int colIndex = 
		localDofMap.getMatIndex(k, it->first) + rStartInterior;
649

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

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

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


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

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

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

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

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

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

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

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

716

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

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

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

	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]));

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

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

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

 	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
749

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

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

	VecDestroy(&tmpVec);
      }

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

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

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

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

779
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
780
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal,
781
		      SAME_NONZERO_PATTERN);
782 783 784 785 786 787
      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);
788
      KSPSetFromOptions(ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
789

Thomas Witkowski's avatar
Thomas Witkowski committed
790
      if (printTimings) {
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
791
	MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
792 793 794 795 796 797 798 799 800 801
	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
802
	MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
803 804 805
	MSG("FETI-DP timing 07: %.5f seconds (factorization of primal schur matrix).\n",
	    MPI::Wtime() - wtime);
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
806
    }
807 808 809 810 811 812 813
  }


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

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

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


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

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