PetscSolverFeti.cc 74.5 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/PetscHelper.h"
16
#include "parallel/PetscSolverFeti.h"
17 18
#include "parallel/PetscSolverFetiDebug.h"
#include "parallel/PetscSolverFetiMonitor.h"
19
#include "parallel/PetscSolverFetiStructs.h"
20 21
#include "parallel/PetscSolverFetiOperators.h"
#include "parallel/PetscSolverFetiTimings.h"
22 23
#include "parallel/StdMpi.h"
#include "parallel/MpiHelper.h"
24
#include "parallel/PetscSolverGlobalMatrix.h"
25
#include "io/VtkWriter.h"
26 27 28 29 30

namespace AMDiS {

  using namespace std;

31 32
  PetscSolverFeti::PetscSolverFeti()
    : PetscSolver(),
33
      schurPrimalSolver(0),
34
      multiLevelTest(false),
Thomas Witkowski's avatar
Thomas Witkowski committed
35
      subdomain(NULL),
36
      massMatrixSolver(NULL),
37 38
      meshLevel(0),
      rStartInterior(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
39
      nGlobalOverallInterior(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
40
      printTimings(false),
Thomas Witkowski's avatar
Thomas Witkowski committed
41
      stokesMode(false),
42
      augmentedLagrange(false),
Thomas Witkowski's avatar
Thomas Witkowski committed
43
      pressureComponent(-1)
44 45 46 47 48 49
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

    string preconditionerName = "";
    Parameters::get("parallel->solver->precon", preconditionerName);
    if (preconditionerName == "" || preconditionerName == "none") {
50
      MSG("Create FETI-DP solver with no preconditioner!\n");
51 52
      fetiPreconditioner = FETI_NONE;
    } else if (preconditionerName == "dirichlet") {
53
      MSG("Create FETI-DP solver with Dirichlet preconditioner!\n");
54 55
      fetiPreconditioner = FETI_DIRICHLET;
    } else if (preconditionerName == "lumped") {
56
      MSG("Create FETI-DP solver with lumped preconditioner!\n");
57 58
      fetiPreconditioner = FETI_LUMPED;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
59 60
      ERROR_EXIT("Preconditioner \"%s\" not available!\n", 
		 preconditionerName.c_str());
61
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
62 63 64 65 66

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

    Parameters::get("parallel->multi level test", multiLevelTest);
69 70
    if (multiLevelTest)
      meshLevel = 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
71 72

    Parameters::get("parallel->print timings", printTimings);
73 74

    Parameters::get("parallel->feti->augmented lagrange", augmentedLagrange);
75 76 77
  }


78
  void PetscSolverFeti::initialize(vector<const FiniteElemSpace*> feSpaces)
79
  {
80 81
    FUNCNAME("PetscSolverFeti::initialize()");

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

85
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();
86
    vector<const FiniteElemSpace*>& uniqueFe = meshDistributor->getFeSpaces();
87

Thomas Witkowski's avatar
Thomas Witkowski committed
88

Thomas Witkowski's avatar
Thomas Witkowski committed
89 90 91 92 93 94 95 96
    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
97 98
    }
			   
Thomas Witkowski's avatar
Thomas Witkowski committed
99 100
    if (subdomain == NULL) {
      subdomain = new PetscSolverGlobalMatrix();
101

102
      if (meshLevel == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
103
	subdomain->setMeshDistributor(meshDistributor, 
104
				      mpiCommGlobal, mpiCommLocal);
105
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
106
	subdomain->setMeshDistributor(meshDistributor, 
107 108
				      levelData.getMpiComm(meshLevel - 1),
				      levelData.getMpiComm(meshLevel));
Thomas Witkowski's avatar
Thomas Witkowski committed
109
	subdomain->setLevel(meshLevel);
110 111
      }
    }
112

113 114 115 116
    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
117

Thomas Witkowski's avatar
Thomas Witkowski committed
118
    if (stokesMode)
Thomas Witkowski's avatar
Thomas Witkowski committed
119
      interfaceDofMap.init(levelData, feSpaces, uniqueFe);
Thomas Witkowski's avatar
Thomas Witkowski committed
120

121
    if (fetiPreconditioner == FETI_DIRICHLET) {
122 123 124
      TEST_EXIT(meshLevel == 0)
	("Dirichlet preconditioner not yet implemented for multilevel FETI-DP\n");

125
      interiorDofMap.init(levelData, feSpaces, uniqueFe, false);
126
    }
127 128 129
  }


Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
130 131 132 133
  void PetscSolverFeti::createDirichletData(Matrix<DOFMatrix*> &mat)
  {
    FUNCNAME("PetscSolverFeti::createDirichletData()");

Thomas Witkowski's avatar
Thomas Witkowski committed
134 135 136 137 138
    bool removeDirichletRows = false;
    Parameters::get("parallel->feti->remove dirichlet", removeDirichletRows);
    if (!removeDirichletRows)
      return;

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
139 140 141 142 143 144 145 146 147
    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)) {
148 149 150 151 152 153 154 155 156
	// === Run test if Dirichlet rows are all the same in all components ===
	// === of one FE space.                                              ===
	TEST_EXIT(dirichletRows[feSpace].size() == dRows.size())
	  ("Wrong number of dirichlet rows!\n");
	
	for (std::set<DegreeOfFreedom>::iterator it0 = dirichletRows[feSpace].begin(),
	       it1 = dRows.begin(); it1 != dRows.end(); ++it0, ++it1) {
	  TEST_EXIT(*it0 == *it1)("Wrong DOFs %d %d!\n", *it0, *it1);
	}
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
157 158 159 160 161 162 163
      } else {
	dirichletRows[feSpace] = dRows;
      }
    }
  }


164 165 166 167
  void PetscSolverFeti::createFetiData()
  {
    FUNCNAME("PetscSolverFeti::createFetiData()");

Thomas Witkowski's avatar
Thomas Witkowski committed
168 169
    double timeCounter = MPI::Wtime();

170 171 172
    TEST_EXIT(meshDistributor)("No mesh distributor object defined!\n");
    TEST_EXIT(meshDistributor->getFeSpaces().size() > 0)
      ("No FE space defined in mesh distributor!\n");
173

174 175
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

176 177 178 179
    primalDofMap.clear();
    dualDofMap.clear();
    lagrangeMap.clear();
    localDofMap.clear();
180
    if (fetiPreconditioner == FETI_DIRICHLET)
181 182
      interiorDofMap.clear();

183 184
    primalDofMap.setDofComm(meshDistributor->getDofComm());
    lagrangeMap.setDofComm(meshDistributor->getDofComm());
185

186 187
    primalDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    dualDofMap.setMpiComm(levelData.getMpiComm(0), 0);
188
    lagrangeMap.setMpiComm(levelData.getMpiComm(0), 0);
189
    localDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
190
    if (fetiPreconditioner == FETI_DIRICHLET)
191
      interiorDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
192

193 194 195
    if (meshLevel == 0)
      localDofMap.setDofComm(meshDistributor->getDofComm());
    else
196 197
      localDofMap.setDofComm(meshDistributor->getDofCommSd());

Thomas Witkowski's avatar
Thomas Witkowski committed
198
    if (stokesMode) {
199 200 201 202 203
      interfaceDofMap.clear();
      interfaceDofMap.setDofComm(meshDistributor->getDofComm());
      interfaceDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    }

204 205
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
206
    
207 208
      createPrimals(feSpace);  

209 210 211
      createDuals(feSpace);

      createInterfaceNodes(feSpace);
212

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
213
      createIndexB(feSpace);     
214
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
215 216 217 218

    primalDofMap.update();
    dualDofMap.update();
    localDofMap.update();
219
    if (fetiPreconditioner == FETI_DIRICHLET)
220
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
221

Thomas Witkowski's avatar
Thomas Witkowski committed
222
    if (stokesMode)
223 224
      interfaceDofMap.update();

225 226 227
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
      createLagrange(feSpace);
228
      createAugmentedLagrange(feSpace);
229
    }
230

231 232 233
    lagrangeMap.update();


234 235 236 237 238 239 240 241 242 243 244 245
    // === ===

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

246
      mpi::getDofNumbering(mpiCommGlobal, groupRowsInterior,
247 248 249 250 251 252 253 254 255
			   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);
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
261
      if (feSpace == pressureFeSpace) {
Thomas Witkowski's avatar
Thomas Witkowski committed
262 263 264 265
	MSG("  nRankInterface = %d  nOverallInterface = %d\n",
	    interfaceDofMap[feSpace].nRankDofs, 
	    interfaceDofMap[feSpace].nOverallDofs);
      } else {
266
	MSG("  nRankPrimals = %d   nLocalPrimals = %d  nOverallPrimals = %d\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
267
	    primalDofMap[feSpace].nRankDofs, 
268
	    primalDofMap[feSpace].nLocalDofs,
Thomas Witkowski's avatar
Thomas Witkowski committed
269 270 271 272 273 274 275 276 277 278
	    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
279 280 281
// 	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
282
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
283
    }
284

Thomas Witkowski's avatar
Thomas Witkowski committed
285
    subdomain->setDofMapping(&localDofMap);
Thomas Witkowski's avatar
Thomas Witkowski committed
286 287 288
    subdomain->setCoarseSpaceDofMapping(&primalDofMap); 
    if (stokesMode)
      subdomain->setCoarseSpaceDofMapping(&interfaceDofMap, pressureComponent);
Thomas Witkowski's avatar
Thomas Witkowski committed
289 290

    if (printTimings) {
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
291
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
292
      timeCounter = MPI::Wtime() - timeCounter;
293
      MSG("FETI-DP timing 01: %.5f seconds (creation of basic data structures)\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
294 295
	  timeCounter);
    }
296 297 298
  }


299
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
300
  {
301
    FUNCNAME("PetscSolverFeti::createPrimals()");  
302

Thomas Witkowski's avatar
Thomas Witkowski committed
303
    if (feSpace == pressureFeSpace)
304 305
      return;

306 307 308
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

Thomas Witkowski's avatar
Thomas Witkowski committed
309
    // Set of DOF indices that are considered to be primal variables.
310
    DofContainerSet& vertices = 
311
      meshDistributor->getBoundaryDofInfo(feSpace, meshLevel).geoDofs[VERTEX];
312 313

    DofIndexSet primals;
314
    for (DofContainerSet::iterator it = vertices.begin(); 
315
	 it != vertices.end(); ++it) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
316 317
      if (dirichletRows[feSpace].count(**it))
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
318

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
319 320
      if (meshLevel == 0) {
	primals.insert(**it);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
321
      } else {
322
	double e = 1e-8;
Thomas Witkowski's avatar
Thomas Witkowski committed
323 324 325
	WorldVector<double> c;
	feSpace->getMesh()->getDofIndexCoords(*it, feSpace, c);

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
326 327 328 329 330 331 332
	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);
	}
      }
333
    }
334

335 336 337 338

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

339
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
340 341
      if (meshDistributor->getDofMap()[feSpace].isRankDof(*it))
	primalDofMap[feSpace].insertRankDof(*it);
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
342
      else
343
  	primalDofMap[feSpace].insertNonRankDof(*it);
344 345 346
  }


347
  void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace)
348 349
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
350

Thomas Witkowski's avatar
Thomas Witkowski committed
351
    if (feSpace == pressureFeSpace)
352 353
      return;

354 355 356
    // === Create global index of the dual nodes on each rank. ===

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

  
377 378 379 380
  void PetscSolverFeti::createInterfaceNodes(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createInterfaceNodes()");

Thomas Witkowski's avatar
Thomas Witkowski committed
381
    if (feSpace != pressureFeSpace)
382 383 384 385 386 387
      return;

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

    for (DofContainer::iterator it = allBoundaryDofs.begin();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
388 389 390 391
	 it != allBoundaryDofs.end(); ++it) {
      if (dirichletRows[feSpace].count(**it))
	continue;      
      
392
      if (meshDistributor->getDofMap()[feSpace].isRankDof(**it))
Thomas Witkowski's avatar
Thomas Witkowski committed
393
	interfaceDofMap[feSpace].insertRankDof(**it);
394
      else
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
395 396
	interfaceDofMap[feSpace].insertNonRankDof(**it);      
    }
397 398 399
  }


400 401 402 403
  void PetscSolverFeti::createLagrange(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
404
    if (feSpace == pressureFeSpace)
405 406
      return;

407 408
    boundaryDofRanks[feSpace].clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
409 410 411
    // 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
412
    map<int, std::set<DegreeOfFreedom> > sdRankDofs;
413

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
414
    if (meshLevel > 0) {
415
      StdMpi<vector<int> > stdMpi(mpiCommGlobal);
416 417 418 419 420 421 422 423 424

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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
450 451 452 453 454 455 456
    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)).         ===

457
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
458 459 460 461 462 463 464
    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);

465
 	  if (meshLevel == 0 ||
Thomas Witkowski's avatar
Thomas Witkowski committed
466
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
467
	    boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());	  
468 469 470 471
	}
      }
    }

472 473 474 475

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

476
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
477

478
    for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), meshLevel, feSpace);
479 480
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
481
	if (!isPrimal(feSpace, it.getDofIndex()))
482 483 484
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
	    stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]);
485 486 487

    stdMpi.updateSendDataSize();

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

      if (recvFromRank)
503
	stdMpi.recv(it.getRank());
504
    }
505

506 507
    stdMpi.startCommunication();

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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
522

523 524 525
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

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


541 542 543 544 545 546 547 548 549
  void PetscSolverFeti::createAugmentedLagrange(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createAugmentedLagrange()");

    if (!augmentedLagrange)
      return;
  }


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
    Vec vec_scale_lagrange;
    lagrangeMap.createVec(vec_scale_lagrange);

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

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

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

	TEST_EXIT_DBG(degree > 1)("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
646
	
647
	int counter = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
648 649 650
	for (int i = 0; i < degree; i++) {
	  for (int j = i + 1; j < degree; j++) {
	    if (W[i] == mpiRank || W[j] == mpiRank) {
651 652 653 654 655
	      MatSetValue(mat_lagrange, 
			  index + counter, 
			  localDofMap.getMatIndex(k, it->first) + rStartInterior,
			  (W[i] == mpiRank ? 1.0 : -1.0),
			  INSERT_VALUES);
656
	    }
657 658 659 660 661 662 663 664 665 666 667 668 669 670
	    counter++;
	  }
	}

	// === Create scaling factors for scaling the lagrange matrix, which ===
	// === is required for FETI-DP preconditioners.                      ===
	
	if (meshDistributor->getDofMap()[feSpaces[k]].isRankDof(it->first)) {
	  int nConstraints = (degree * (degree - 1)) / 2;
	  for (int i = 0; i < nConstraints; i++) {
	    VecSetValue(vec_scale_lagrange,
			index + i,
			1.0 / static_cast<double>(degree),
			INSERT_VALUES);
671 672 673 674 675 676 677
	  }
	}
      }
    }

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

679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694

    // === If required, create \ref mat_lagrange_scaled ===

    VecAssemblyBegin(vec_scale_lagrange);
    VecAssemblyEnd(vec_scale_lagrange);

    if (fetiPreconditioner != FETI_NONE || stokesMode) {
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatDiagonalScale(mat_lagrange_scaled, vec_scale_lagrange, PETSC_NULL);
    }

    VecDestroy(&vec_scale_lagrange);


    // === Print final timings. ===

Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
695 696
    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
697 698
      MSG("FETI-DP timing 05: %.5f seconds (creation of lagrange constraint matrix)\n", 
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
699
    }
700 701 702
  }


703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728
  void PetscSolverFeti::createMatAugmentedLagrange(vector<const FiniteElemSpace*> &feSpaces)
  {
    FUNCNAME("PetscSolverFeti::createMatAugmentedLagrange()");

    if (!augmentedLagrange)
      return;

    double wtime = MPI::Wtime();

    nRankEdges = 0;
    nOverallEdges = 0;
    InteriorBoundary &intBound = meshDistributor->getIntBoundary();
    for (InteriorBoundary::iterator it(intBound.getOwn()); !it.end(); ++it)
      if (it->rankObj.subObj == EDGE)
	nRankEdges++;
    
    int rStartEdges = 0;
    mpi::getDofNumbering(mpiCommGlobal, nRankEdges, rStartEdges, nOverallEdges);

    MSG("nRankEdges = %d, nOverallEdges = %d\n", nRankEdges, nOverallEdges);

    nRankEdges *= feSpaces.size();
    rStartEdges *= feSpaces.size();
    nOverallEdges *= feSpaces.size();

    MatCreateAIJ(mpiCommGlobal,
729 730
		 nRankEdges, lagrangeMap.getRankDofs(),
		 nOverallEdges, lagrangeMap.getOverallDofs(),
731 732 733 734
		 1, PETSC_NULL, 1, PETSC_NULL, 
		 &mat_augmented_lagrange);
    MatSetOption(mat_augmented_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

735
    int rowCounter = rStartEdges;
736 737 738 739 740 741 742 743 744 745 746 747 748
    for (InteriorBoundary::iterator it(intBound.getOwn()); !it.end(); ++it) {
      if (it->rankObj.subObj == EDGE) {
	for (int i = 0; i < feSpaces.size(); i++) {
	  DofContainer edgeDofs;
	  it->rankObj.el->getAllDofs(feSpaces[i], it->rankObj, edgeDofs);

	  MSG("SIZE = %d\n", edgeDofs.size());

	  for (DofContainer::iterator it = edgeDofs.begin();
	       it != edgeDofs.end(); it++) {
	    TEST_EXIT_DBG(isPrimal(feSpaces[i], **it) == false)
	      ("Should not be primal!\n");

749 750 751
	    if (dirichletRows[feSpaces[i]].count(**it))
	      continue;

752
	    int col = lagrangeMap.getMatIndex(i, **it);
753
	    double value = 1.0;
754
	    MatSetValue(mat_augmented_lagrange, rowCounter, col, value, INSERT_VALUES);
755 756
	  }

757
	  rowCounter++;
758 759 760 761 762 763 764 765 766 767 768 769 770 771 772
	}
      }
    }

    MatAssemblyBegin(mat_augmented_lagrange, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(mat_augmented_lagrange, MAT_FINAL_ASSEMBLY);

    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
      MSG("FETI-DP timing 05a: %.5f seconds (creation of augmented lagrange constraint matrix)\n", 
	  MPI::Wtime() - wtime);
    }
  }


773
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
774
  {
775
    FUNCNAME("PetscSolverFeti::createSchurPrimalKsp()");
776

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

780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796
      if (augmentedLagrange == false) {
	schurPrimalData.subSolver = subdomain;
	
	localDofMap.createVec(schurPrimalData.tmp_vec_b, nGlobalOverallInterior);
	primalDofMap.createVec(schurPrimalData.tmp_vec_primal);
	
	MatCreateShell(mpiCommGlobal,
		       primalDofMap.getRankDofs(), 
		       primalDofMap.getRankDofs(), 
		       primalDofMap.getOverallDofs(), 
		       primalDofMap.getOverallDofs(),
		       &schurPrimalData, 
		       &mat_schur_primal);
	MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			     (void(*)(void))petscMultMatSchurPrimal);	
      } else {
	schurPrimalAugmentedData.subSolver = subdomain;
797
	schurPrimalAugmentedData.nestedVec = true;
798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816

	localDofMap.createVec(schurPrimalAugmentedData.tmp_vec_b0, nGlobalOverallInterior);
	localDofMap.createVec(schurPrimalAugmentedData.tmp_vec_b1, nGlobalOverallInterior);
	primalDofMap.createVec(schurPrimalAugmentedData.tmp_vec_primal);
	lagrangeMap.createVec(schurPrimalAugmentedData.tmp_vec_lagrange);

	schurPrimalAugmentedData.mat_lagrange = &mat_lagrange;
	schurPrimalAugmentedData.mat_augmented_lagrange = &mat_augmented_lagrange;

	MatCreateShell(mpiCommGlobal,
		       primalDofMap.getRankDofs() + nRankEdges, 
		       primalDofMap.getRankDofs() + nRankEdges, 
		       primalDofMap.getOverallDofs() + nOverallEdges, 
		       primalDofMap.getOverallDofs() + nOverallEdges,
		       &schurPrimalAugmentedData, 
		       &mat_schur_primal);
	MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			     (void(*)(void))petscMultMatSchurPrimalAugmented);
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
817

818
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
819
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
820 821
      KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
      KSPSetType(ksp_schur_primal, KSPGMRES);
Thomas Witkowski's avatar
Thomas Witkowski committed
822 823
      KSPSetFromOptions(ksp_schur_primal);
    } else {
824 825
      MSG("Create direct schur primal solver!\n");

826 827
      double wtime = MPI::Wtime();

828 829 830
      TEST_EXIT_DBG(meshLevel == 0)
	("Does not support for multilevel, check usage of localDofMap.\n");

831

832
      // === Create explicit matrix representation of the Schur primal system. ===
833

834 835 836 837
      if (!augmentedLagrange)
	createMatExplicitSchurPrimal();
      else
	createMatExplicitAugmentedSchurPrimal();
838 839 840


      // === Create KSP solver object and set appropriate solver options. ====
841

842
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
843
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal,
844
		      SAME_NONZERO_PATTERN);
845 846 847 848 849 850
      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);
851
      KSPSetFromOptions(ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
852

853 854 855

      // === And finally print timings, if required. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
856
      if (printTimings) {
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
857
	MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
858 859 860 861 862 863 864 865 866 867
	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
868
	MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
869 870 871
	MSG("FETI-DP timing 07: %.5f seconds (factorization of primal schur matrix).\n",
	    MPI::Wtime() - wtime);
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
872
    }