PetscSolverFeti.cc 71.3 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
    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");

Thomas Witkowski's avatar
Thomas Witkowski committed
115 116
      MSG("========= !!!! SUPER WARNING !!! ========\n");
      MSG("Use make use of stokes mode which is work in progress! We guarantee for nothing and even less than this!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
117
      pressureFeSpace = feSpaces[pressureComponent];
Thomas Witkowski's avatar
Thomas Witkowski committed
118 119
    }
			   
Thomas Witkowski's avatar
Thomas Witkowski committed
120 121
    if (subdomain == NULL) {
      subdomain = new PetscSolverGlobalMatrix();
122

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

134 135 136 137
    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
138

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

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

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


Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171
  void PetscSolverFeti::createDirichletData(Matrix<DOFMatrix*> &mat)
  {
    FUNCNAME("PetscSolverFeti::createDirichletData()");

    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;
      }
    }
  }


172 173 174 175
  void PetscSolverFeti::createFetiData()
  {
    FUNCNAME("PetscSolverFeti::createFetiData()");

Thomas Witkowski's avatar
Thomas Witkowski committed
176 177
    double timeCounter = MPI::Wtime();

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

182 183
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

184 185 186 187
    primalDofMap.clear();
    dualDofMap.clear();
    lagrangeMap.clear();
    localDofMap.clear();
188
    if (fetiPreconditioner == FETI_DIRICHLET)
189 190
      interiorDofMap.clear();

191 192
    primalDofMap.setDofComm(meshDistributor->getDofComm());
    lagrangeMap.setDofComm(meshDistributor->getDofComm());
193

194 195
    primalDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    dualDofMap.setMpiComm(levelData.getMpiComm(0), 0);
196
    lagrangeMap.setMpiComm(levelData.getMpiComm(0), 0);
197
    localDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
198
    if (fetiPreconditioner == FETI_DIRICHLET)
199
      interiorDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
200

201 202 203
    if (meshLevel == 0)
      localDofMap.setDofComm(meshDistributor->getDofComm());
    else
204 205
      localDofMap.setDofComm(meshDistributor->getDofCommSd());

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

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

217 218 219
      createDuals(feSpace);

      createInterfaceNodes(feSpace);
220

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
221
      createIndexB(feSpace);     
222
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
223 224 225 226

    primalDofMap.update();
    dualDofMap.update();
    localDofMap.update();
227
    if (fetiPreconditioner == FETI_DIRICHLET)
228
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
229

Thomas Witkowski's avatar
Thomas Witkowski committed
230
    if (stokesMode)
231 232
      interfaceDofMap.update();

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

238 239 240
    lagrangeMap.update();


241 242 243 244 245 246 247 248 249 250 251 252
    // === ===

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

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

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

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

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

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


305
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
306
  {
307
    FUNCNAME("PetscSolverFeti::createPrimals()");  
308

Thomas Witkowski's avatar
Thomas Witkowski committed
309
    if (feSpace == pressureFeSpace)
310 311
      return;

312 313 314
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

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

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

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

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

341 342 343 344

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

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


353
  void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace)
354 355
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
356

Thomas Witkowski's avatar
Thomas Witkowski committed
357
    if (feSpace == pressureFeSpace)
358 359
      return;

360 361 362
    // === Create global index of the dual nodes on each rank. ===

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

  
383 384 385 386
  void PetscSolverFeti::createInterfaceNodes(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createInterfaceNodes()");

Thomas Witkowski's avatar
Thomas Witkowski committed
387
    if (feSpace != pressureFeSpace)
388 389 390 391 392 393
      return;

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

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


406 407 408 409
  void PetscSolverFeti::createLagrange(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
410
    if (feSpace == pressureFeSpace)
411 412
      return;

413 414
    boundaryDofRanks[feSpace].clear();

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

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

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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
456 457 458 459 460 461 462
    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)).         ===

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

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

478 479 480 481

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

482
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
483

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

    stdMpi.updateSendDataSize();

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

      if (recvFromRank)
509
	stdMpi.recv(it.getRank());
510
    }
511

512 513
    stdMpi.startCommunication();

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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
528

529 530 531
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

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


547
  void PetscSolverFeti::createIndexB(const FiniteElemSpace *feSpace)
548
  {
549
    FUNCNAME("PetscSolverFeti::createIndexB()");
550

551
    DOFAdmin* admin = feSpace->getAdmin();
552 553 554 555

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

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

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

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


600
  void PetscSolverFeti::createMatLagrange(vector<const FiniteElemSpace*> &feSpaces)
601 602 603
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
604
    double wtime = MPI::Wtime();
605
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Thomas Witkowski committed
606

607 608
    // === Create distributed matrix for Lagrange constraints. ===

609 610 611 612 613
    MatCreateAIJ(mpiCommGlobal,
		 lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		 lagrangeMap.getOverallDofs(), nGlobalOverallInterior,
		 2, PETSC_NULL, 2, PETSC_NULL,
		 &mat_lagrange);
614
    MatSetOption(mat_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
615

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

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

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

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

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

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

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


667
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
668
  {
669
    FUNCNAME("PetscSolverFeti::createSchurPrimalKsp()");
670

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

Thomas Witkowski's avatar
Thomas Witkowski committed
674
      schurPrimalData.subSolver = subdomain;
Thomas Witkowski's avatar
Thomas Witkowski committed
675
      
676 677
      localDofMap.createVec(schurPrimalData.tmp_vec_b, nGlobalOverallInterior);
      primalDofMap.createVec(schurPrimalData.tmp_vec_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
678

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

697 698
      double wtime = MPI::Wtime();

699 700 701
      TEST_EXIT_DBG(meshLevel == 0)
	("Does not support for multilevel, check usage of localDofMap.\n");

702 703 704
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
      int nRowsRankB = localDofMap.getRankDofs();
705

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

713

714 715 716
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
717

718 719
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

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

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

728
	MatRestoreRow(subdomain->getMatInteriorCoarse(), row, &nCols, &cols, &values);
729 730
      }

731
      TEST_EXIT(static_cast<int>(mat_b_primal_cols.size()) == 
732
		primalDofMap.getLocalDofs())
733
	("Should not happen!\n");
734

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

 	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
746

Thomas Witkowski's avatar
Thomas Witkowski committed
747
	subdomain->solve(tmpVec, tmpVec);
748

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

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
762 763
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
764

765 766
      MatDuplicate(subdomain->getMatCoarse(), MAT_COPY_VALUES, 
		   &mat_schur_primal);
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
767 768

      Mat matPrimal;
769 770
      MatMatMult(subdomain->getMatCoarseInterior(), matBPi, MAT_INITIAL_MATRIX, 
		 PETSC_DEFAULT, &matPrimal);
771
      MatAXPY(mat_schur_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);
772 773

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
774
      MatDestroy(&matBPi);
775

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

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


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

811
    if (schurPrimalSolver == 0) {
812
      schurPrimalData.subSolver = NULL;
813

814 815 816
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
    }
817 818 819
    
    MatDestroy(&mat_schur_primal);
    KSPDestroy(&ksp_schur_primal);
820 821 822
  }


823
  void PetscSolverFeti::createFetiKsp(vector<const FiniteElemSpace*> &feSpaces)
824 825 826
  {
    FUNCNAME("PetscSolverFeti::createFetiKsp()");

827 828
    // === Create FETI-DP solver object. ===

829 830 831 832
    fetiData.mat_lagrange = &mat_lagrange;
    fetiData.subSolver = subdomain;
    fetiData.ksp_schur_primal = &ksp_schur_primal;

Thomas Witkowski's avatar
Thomas Witkowski committed
833
    localDofMap.createVec(fetiData.tmp_vec_b0, nGlobalOverallInterior);
834
    lagrangeMap.createVec(fetiData.tmp_vec_lagrange);
Thomas Witkowski's avatar
Thomas Witkowski committed
835
    primalDofMap.createVec(fetiData.tmp_vec_primal0);
836

Thomas Witkowski's avatar
Thomas Witkowski committed
837
    if (stokesMode == false) {      
838 839 840 841 842 843 844 845
      MatCreateShell(mpiCommGlobal,
		     lagrangeMap.getRankDofs(), 
		     lagrangeMap.getRankDofs(),
		     lagrangeMap.getOverallDofs(), 
		     lagrangeMap.getOverallDofs(),
		     &fetiData, &mat_feti);
      MatShellSetOperation(mat_feti, MATOP_MULT, 
			   (void(*)(void))petscMultMatFeti);
846
    }  else {
Thomas Witkowski's avatar
Thomas Witkowski committed
847 848 849 850
      localDofMap.createVec(fetiData.tmp_vec_b1, nGlobalOverallInterior);
      primalDofMap.createVec(fetiData.tmp_vec_primal1);
      interfaceDofMap.createVec(fetiData.tmp_vec_interface);

851 852 853 854 855
      MatCreateShell(mpiCommGlobal,
		     interfaceDofMap.getRankDofs() + lagrangeMap.getRankDofs(), 
		     interfaceDofMap.getRankDofs() + lagrangeMap.getRankDofs(),
		     interfaceDofMap.getOverallDofs() + lagrangeMap.getOverallDofs(), 
		     interfaceDofMap.getOverallDofs() + lagrangeMap.getOverallDofs(),
856
		     &fetiData, &mat_feti);
857
      MatShellSetOperation(mat_feti, MATOP_MULT, 
Thomas Witkowski's avatar
Thomas Witkowski committed
858
			   (void(*)(void))petscMultMatFetiInterface);      
859
    }