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

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

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

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

    Parameters::get("parallel->multi level test", multiLevelTest);
64 65
    if (multiLevelTest)
      meshLevel = 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
66 67

    Parameters::get("parallel->print timings", printTimings);
68 69 70
  }


71
  void PetscSolverFeti::initialize(vector<const FiniteElemSpace*> feSpaces)
72
  {
73 74
    FUNCNAME("PetscSolverFeti::initialize()");

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

78
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();
79
    vector<const FiniteElemSpace*>& uniqueFe = meshDistributor->getFeSpaces();
80

Thomas Witkowski's avatar
Thomas Witkowski committed
81

Thomas Witkowski's avatar
Thomas Witkowski committed
82 83 84 85 86 87 88
    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
89 90
      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
91
      pressureFeSpace = feSpaces[pressureComponent];
Thomas Witkowski's avatar
Thomas Witkowski committed
92 93
    }
			   
Thomas Witkowski's avatar
Thomas Witkowski committed
94 95
    if (subdomain == NULL) {
      subdomain = new PetscSolverGlobalMatrix();
96

97
      if (meshLevel == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
98
	subdomain->setMeshDistributor(meshDistributor, 
99
				      mpiCommGlobal, mpiCommLocal);
100
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
101
	subdomain->setMeshDistributor(meshDistributor, 
102 103
				      levelData.getMpiComm(meshLevel - 1),
				      levelData.getMpiComm(meshLevel));
Thomas Witkowski's avatar
Thomas Witkowski committed
104
	subdomain->setLevel(meshLevel);
105 106
      }
    }
107

108 109 110 111
    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
112

Thomas Witkowski's avatar
Thomas Witkowski committed
113
    if (stokesMode)
Thomas Witkowski's avatar
Thomas Witkowski committed
114
      interfaceDofMap.init(levelData, feSpaces, uniqueFe);
Thomas Witkowski's avatar
Thomas Witkowski committed
115

116
    if (fetiPreconditioner == FETI_DIRICHLET) {
117 118 119
      TEST_EXIT(meshLevel == 0)
	("Dirichlet preconditioner not yet implemented for multilevel FETI-DP\n");

120
      interiorDofMap.init(levelData, feSpaces, uniqueFe, false);
121
    }
122 123 124 125 126 127 128
  }


  void PetscSolverFeti::createFetiData()
  {
    FUNCNAME("PetscSolverFeti::createFetiData()");

Thomas Witkowski's avatar
Thomas Witkowski committed
129 130
    double timeCounter = MPI::Wtime();

131 132 133
    TEST_EXIT(meshDistributor)("No mesh distributor object defined!\n");
    TEST_EXIT(meshDistributor->getFeSpaces().size() > 0)
      ("No FE space defined in mesh distributor!\n");
134

135 136
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

137 138 139 140
    primalDofMap.clear();
    dualDofMap.clear();
    lagrangeMap.clear();
    localDofMap.clear();
141
    if (fetiPreconditioner == FETI_DIRICHLET)
142 143
      interiorDofMap.clear();

144 145
    primalDofMap.setDofComm(meshDistributor->getDofComm());
    lagrangeMap.setDofComm(meshDistributor->getDofComm());
146

147 148
    primalDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    dualDofMap.setMpiComm(levelData.getMpiComm(0), 0);
149
    lagrangeMap.setMpiComm(levelData.getMpiComm(0), 0);
150
    localDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
151
    if (fetiPreconditioner == FETI_DIRICHLET)
152
      interiorDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
153

154 155 156
    if (meshLevel == 0)
      localDofMap.setDofComm(meshDistributor->getDofComm());
    else
157 158
      localDofMap.setDofComm(meshDistributor->getDofCommSd());

Thomas Witkowski's avatar
Thomas Witkowski committed
159
    if (stokesMode) {
160 161 162 163 164
      interfaceDofMap.clear();
      interfaceDofMap.setDofComm(meshDistributor->getDofComm());
      interfaceDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    }

165 166
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
167
    
168 169
      createPrimals(feSpace);  

170 171 172
      createDuals(feSpace);

      createInterfaceNodes(feSpace);
173

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
174
      createIndexB(feSpace);     
175
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
176 177 178 179

    primalDofMap.update();
    dualDofMap.update();
    localDofMap.update();
180
    if (fetiPreconditioner == FETI_DIRICHLET)
181
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
182

Thomas Witkowski's avatar
Thomas Witkowski committed
183
    if (stokesMode)
184 185
      interfaceDofMap.update();

186 187 188 189
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
      createLagrange(feSpace);
    }
190

191 192 193
    lagrangeMap.update();


194 195 196 197 198 199 200 201 202 203 204 205
    // === ===

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

206
      mpi::getDofNumbering(mpiCommGlobal, groupRowsInterior,
207 208 209 210 211 212 213 214 215
			   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);
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
221
      if (feSpace == pressureFeSpace) {
Thomas Witkowski's avatar
Thomas Witkowski committed
222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241
	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);

	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
242
    }
243

Thomas Witkowski's avatar
Thomas Witkowski committed
244
    subdomain->setDofMapping(&localDofMap);
Thomas Witkowski's avatar
Thomas Witkowski committed
245 246 247
    subdomain->setCoarseSpaceDofMapping(&primalDofMap); 
    if (stokesMode)
      subdomain->setCoarseSpaceDofMapping(&interfaceDofMap, pressureComponent);
Thomas Witkowski's avatar
Thomas Witkowski committed
248 249

    if (printTimings) {
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
250
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
251
      timeCounter = MPI::Wtime() - timeCounter;
252
      MSG("FETI-DP timing 01: %.5f seconds (creation of basic data structures)\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
253 254
	  timeCounter);
    }
255 256 257
  }


258
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
259
  {
260
    FUNCNAME("PetscSolverFeti::createPrimals()");  
261

Thomas Witkowski's avatar
Thomas Witkowski committed
262
    if (feSpace == pressureFeSpace)
263 264
      return;

265 266 267
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

Thomas Witkowski's avatar
Thomas Witkowski committed
268
    // Set of DOF indices that are considered to be primal variables.
269
    DofContainerSet& vertices = 
270
      meshDistributor->getBoundaryDofInfo(feSpace, meshLevel).geoDofs[VERTEX];
271 272

    DofIndexSet primals;
273
    for (DofContainerSet::iterator it = vertices.begin(); 
274
	 it != vertices.end(); ++it) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
275
      if (meshLevel == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
276 277 278 279 280 281 282 283
	MSG("WARNING: Modified primal detection algorithm!\n");

	double e = 1e-3;
	WorldVector<double> c;
	feSpace->getMesh()->getDofIndexCoords(*it, feSpace, c);

	if (fabs(c[0]) > e && fabs(c[1]) > e && fabs(c[0] - 1.0) > e && fabs(c[1] - 1.0) > e)
	  primals.insert(**it);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
284
      } else {
285
	double e = 1e-8;
Thomas Witkowski's avatar
Thomas Witkowski committed
286 287 288
	WorldVector<double> c;
	feSpace->getMesh()->getDofIndexCoords(*it, feSpace, c);

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
289 290 291 292 293 294 295
	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);
	}
      }
296
    }
297

298 299 300 301

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

302
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
303 304
      if (meshDistributor->getDofMap()[feSpace].isRankDof(*it))
	primalDofMap[feSpace].insertRankDof(*it);
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
305
      else
306
  	primalDofMap[feSpace].insertNonRankDof(*it);
307 308 309
  }


310
  void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace)
311 312
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
313

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

317 318 319
    // === Create global index of the dual nodes on each rank. ===

    DofContainer allBoundaryDofs;
320
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
321 322 323 324

    for (DofContainer::iterator it = allBoundaryDofs.begin();
	 it != allBoundaryDofs.end(); ++it)
      if (!isPrimal(feSpace, **it))
325 326
	if (meshLevel == 0) {
	  dualDofMap[feSpace].insertRankDof(**it);
327 328 329 330
 	} else {
 	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(**it))
 	    dualDofMap[feSpace].insertRankDof(**it);
 	}	  
331 332 333
  }

  
334 335 336 337
  void PetscSolverFeti::createInterfaceNodes(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createInterfaceNodes()");

Thomas Witkowski's avatar
Thomas Witkowski committed
338
    if (feSpace != pressureFeSpace)
339 340 341 342 343 344 345
      return;

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

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


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

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

360 361
    boundaryDofRanks[feSpace].clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
362 363 364
    // 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
365
    map<int, std::set<DegreeOfFreedom> > sdRankDofs;
366

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
367
    if (meshLevel > 0) {
368
      StdMpi<vector<int> > stdMpi(mpiCommGlobal);
369 370 371 372 373 374 375 376 377

      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
378
	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))
379
	    subdomainRankDofs.push_back(1);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
380
	  else
381 382 383 384 385 386 387 388 389 390 391 392 393 394 395
	    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
396 397 398 399 400 401
	   !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());
    }
402

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
403 404 405 406 407 408 409
    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)).         ===

410
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
411 412 413 414 415 416 417
    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);

418
 	  if (meshLevel == 0 ||
Thomas Witkowski's avatar
Thomas Witkowski committed
419
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
420
	    boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());	  
421 422 423 424
	}
      }
    }

425 426 427 428

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

429
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
430

431
    for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), meshLevel, feSpace);
432 433
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
434
	if (!isPrimal(feSpace, it.getDofIndex()))
435 436 437
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
	    stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]);
438 439 440

    stdMpi.updateSendDataSize();

441
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
442
	 !it.end(); it.nextRank()) {
443
      bool recvFromRank = false;
444
      for (; !it.endDofIter(); it.nextDof()) {
445
	if (!isPrimal(feSpace, it.getDofIndex())) {
446 447 448
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
 	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
449 450 451
	    recvFromRank = true;
	    break;
	  }
452
	}
453
      }
454 455

      if (recvFromRank)
456
	stdMpi.recv(it.getRank());
457
    }
458

459 460
    stdMpi.startCommunication();

461
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
462
	 !it.end(); it.nextRank()) {
463
      int i = 0;
464
      for (; !it.endDofIter(); it.nextDof())
465
	if (!isPrimal(feSpace, it.getDofIndex()))
466 467 468
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
 	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex())))	    
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
469 470
	    boundaryDofRanks[feSpace][it.getDofIndex()] = 
	      stdMpi.getRecvData(it.getRank())[i++];
471 472
	  else
	    lagrangeMap[feSpace].insertNonRankDof(it.getDofIndex());
473 474
    }

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
475

476 477 478
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

479
    int nRankLagrange = 0;
480
    DofMap& dualMap = dualDofMap[feSpace].getMap();
481
    for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
482 483
      if (meshDistributor->getDofMap()[feSpace].isRankDof(it->first)) {
	lagrangeMap[feSpace].insertRankDof(it->first, nRankLagrange);
484
	int degree = boundaryDofRanks[feSpace][it->first].size();
485
	nRankLagrange += (degree * (degree - 1)) / 2;
Thomas Witkowski's avatar
Thomas Witkowski committed
486
      } else {
487
	lagrangeMap[feSpace].insertNonRankDof(it->first);
488 489
      }
    }
490
    lagrangeMap[feSpace].nRankDofs = nRankLagrange;
491 492 493
  }


494
  void PetscSolverFeti::createIndexB(const FiniteElemSpace *feSpace)
495
  {
496
    FUNCNAME("PetscSolverFeti::createIndexB()");
497

498
    DOFAdmin* admin = feSpace->getAdmin();
499 500 501 502

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

504
    int nLocalInterior = 0;    
505
    for (int i = 0; i < admin->getUsedSize(); i++) {
506
      if (admin->isDofFree(i) == false && 
507
	  isPrimal(feSpace, i) == false &&
508 509
	  isDual(feSpace, i) == false &&
	  isInterface(feSpace, i) == false) {
510 511
	if (meshLevel == 0) {
	  localDofMap[feSpace].insertRankDof(i, nLocalInterior);
512

513
	  if (fetiPreconditioner == FETI_DIRICHLET)
514 515 516 517 518 519 520 521
	    interiorDofMap[feSpace].insertRankDof(i, nLocalInterior);
	  
	  nLocalInterior++;	  
	} else {
	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(i))
	    localDofMap[feSpace].insertRankDof(i);
	  else
	    localDofMap[feSpace].insertNonRankDof(i);
522 523 524

	  TEST_EXIT_DBG(fetiPreconditioner == FETI_NONE)
	    ("Not yet implemnted!\n");
525
	}
526
      }
527
    }
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
528
    
529 530
    // === And finally, add the global indicies of all dual nodes. ===

531 532 533 534 535 536 537 538 539 540
    for (DofMap::iterator it = dualDofMap[feSpace].getMap().begin();
	 it != dualDofMap[feSpace].getMap().end(); ++it)
      if (meshLevel == 0)
	localDofMap[feSpace].insertRankDof(it->first);
      else {
	if (meshDistributor->getDofMapSd()[feSpace].isRankDof(it->first))
	  localDofMap[feSpace].insertRankDof(it->first);
	else 
	  localDofMap[feSpace].insertNonRankDof(it->first);
      }
541 542 543
  }


544
  void PetscSolverFeti::createMatLagrange(vector<const FiniteElemSpace*> &feSpaces)
545 546 547
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
548
    double wtime = MPI::Wtime();
549
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Thomas Witkowski committed
550

551 552
    // === Create distributed matrix for Lagrange constraints. ===

553 554 555 556 557
    MatCreateAIJ(mpiCommGlobal,
		 lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		 lagrangeMap.getOverallDofs(), nGlobalOverallInterior,
		 2, PETSC_NULL, 2, PETSC_NULL,
		 &mat_lagrange);
558
    MatSetOption(mat_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
559

560 561 562 563 564 565 566
    // === 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
567
    for (unsigned int k = 0; k < feSpaces.size(); k++) {
568
      DofMap &dualMap = dualDofMap[feSpaces[k]].getMap();
569

570
      for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
571
	TEST_EXIT_DBG(boundaryDofRanks[feSpaces[k]].count(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
572 573 574
	  ("Should not happen!\n");
	
	// Global index of the first Lagrange constriant for this node.
575
	int index = lagrangeMap.getMatIndex(k, it->first);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
576

Thomas Witkowski's avatar
Thomas Witkowski committed
577
	// Copy set of all ranks that contain this dual node.
578 579
	vector<int> W(boundaryDofRanks[feSpaces[k]][it->first].begin(), 
		      boundaryDofRanks[feSpaces[k]][it->first].end());
Thomas Witkowski's avatar
Thomas Witkowski committed
580 581
	// Number of ranks that contain this dual node.
	int degree = W.size();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
582 583

	TEST_EXIT_DBG(degree > 1)("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
584 585 586 587
	
	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
588 589
	      int colIndex = 
		localDofMap.getMatIndex(k, it->first) + rStartInterior;
590

591
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
592
	      MatSetValue(mat_lagrange, index, colIndex, value, INSERT_VALUES);
593
	    }
Thomas Witkowski's avatar
Thomas Witkowski committed
594
	    index++;	      
595 596 597 598 599 600 601
	  }
	}
      }
    }

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

Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
603 604
    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
605 606
      MSG("FETI-DP timing 05: %.5f seconds (creation of lagrange constraint matrix)\n", 
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
607
    }
608 609 610
  }


611
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
612
  {
613
    FUNCNAME("PetscSolverFeti::createSchurPrimalKsp()");
614

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

Thomas Witkowski's avatar
Thomas Witkowski committed
618
      schurPrimalData.subSolver = subdomain;
Thomas Witkowski's avatar
Thomas Witkowski committed
619
      
620 621
      localDofMap.createVec(schurPrimalData.tmp_vec_b, nGlobalOverallInterior);
      primalDofMap.createVec(schurPrimalData.tmp_vec_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
622

623
      MatCreateShell(mpiCommGlobal,
624 625 626 627
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getOverallDofs(), 
		     primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
628 629 630 631 632
		     &schurPrimalData, 
		     &mat_schur_primal);
      MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			   (void(*)(void))petscMultMatSchurPrimal);
      
633
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
634
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
635 636
      KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
      KSPSetType(ksp_schur_primal, KSPGMRES);
Thomas Witkowski's avatar
Thomas Witkowski committed
637 638
      KSPSetFromOptions(ksp_schur_primal);
    } else {
639 640
      MSG("Create direct schur primal solver!\n");

641 642
      double wtime = MPI::Wtime();

643 644 645
      TEST_EXIT_DBG(meshLevel == 0)
	("Does not support for multilevel, check usage of localDofMap.\n");

646 647 648
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
      int nRowsRankB = localDofMap.getRankDofs();
649

Thomas Witkowski's avatar
Thomas Witkowski committed
650
      Mat matBPi;
651 652 653
      MatCreateAIJ(mpiCommGlobal,
		   nRowsRankB, nRowsRankPrimal, 
		   nGlobalOverallInterior, nRowsOverallPrimal,
654 655 656
		   150, PETSC_NULL, 150, PETSC_NULL, &matBPi);
      MatSetOption(matBPi, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

657

658 659 660
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
661

662 663
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

664
      for (int i = 0; i < nRowsRankB; i++) {
665
	PetscInt row = localDofMap.getStartDofs() + i;
666
	MatGetRow(subdomain->getMatInteriorCoarse(), row, &nCols, &cols, &values);
667 668 669 670 671

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

672
	MatRestoreRow(subdomain->getMatInteriorCoarse(), row, &nCols, &cols, &values);
673 674
      }

675
      TEST_EXIT(static_cast<int>(mat_b_primal_cols.size()) == 
676
		primalDofMap.getLocalDofs())
677
	("Should not happen!\n");
678

679 680 681
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
682
	VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);
683 684 685 686 687 688 689

 	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);
690
	
Thomas Witkowski's avatar
Thomas Witkowski committed
691
	subdomain->solve(tmpVec, tmpVec);
692

Thomas Witkowski's avatar
Thomas Witkowski committed
693
	PetscScalar *tmpValues;
694
	VecGetArray(tmpVec, &tmpValues);
695
	for (int i  = 0; i < nRowsRankB; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
696
	  MatSetValue(matBPi, 
697
		      localDofMap.getStartDofs() + i,
Thomas Witkowski's avatar
Thomas Witkowski committed
698 699 700
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
701 702 703 704 705
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
706 707
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
708

709 710
      MatDuplicate(subdomain->getMatCoarse(), MAT_COPY_VALUES, 
		   &mat_schur_primal);
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
711 712

      Mat matPrimal;
713 714
      MatMatMult(subdomain->getMatCoarseInterior(), matBPi, MAT_INITIAL_MATRIX, 
		 PETSC_DEFAULT, &matPrimal);
715
      MatAXPY(mat_schur_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);
716 717

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
718
      MatDestroy(&matBPi);
719

720
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
721
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal,
722
		      SAME_NONZERO_PATTERN);
723 724 725 726 727 728
      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);
729
      KSPSetFromOptions(ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
730

Thomas Witkowski's avatar
Thomas Witkowski committed
731
      if (printTimings) {
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
732
	MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
733 734 735 736 737 738 739 740 741 742
	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
743
	MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
744 745 746
	MSG("FETI-DP timing 07: %.5f seconds (factorization of primal schur matrix).\n",
	    MPI::Wtime() - wtime);
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
747
    }
748 749 750 751 752 753 754
  }


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

755
    if (schurPrimalSolver == 0) {
756
      schurPrimalData.subSolver = NULL;
757

758 759 760
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
    }
761 762 763
    
    MatDestroy(&mat_schur_primal);
    KSPDestroy(&ksp_schur_primal);
764 765 766
  }


767
  void PetscSolverFeti::createFetiKsp(vector<const FiniteElemSpace*> &feSpaces)
768 769 770
  {
    FUNCNAME("PetscSolverFeti::createFetiKsp()");

771 772
    // === Create FETI-DP solver object. ===

773 774 775 776
    fetiData.mat_lagrange = &mat_lagrange;
    fetiData.subSolver = subdomain;
    fetiData.ksp_schur_primal = &ksp_schur_primal;

Thomas Witkowski's avatar
Thomas Witkowski committed
777
    localDofMap.createVec(fetiData.tmp_vec_b0, nGlobalOverallInterior);
778
    lagrangeMap.createVec(fetiData.tmp_vec_lagrange);
Thomas Witkowski's avatar
Thomas Witkowski committed
779
    primalDofMap.createVec(fetiData.tmp_vec_primal0);
780

Thomas Witkowski's avatar
Thomas Witkowski committed
781
    if (stokesMode == false) {      
782 783 784 785 786 787 788 789
      MatCreateShell(mpiCommGlobal,
		     lagrangeMap.getRankDofs(), 
		     lagrangeMap.getRankDofs(),
		     lagrangeMap.getOverallDofs(), 
		     lagrangeMap.getOverallDofs(),
		     &fetiData, &mat_feti);
      MatShellSetOperation(mat_feti, MATOP_MULT, 
			   (void(*)(void))petscMultMatFeti);
790
    }  else {
Thomas Witkowski's avatar
Thomas Witkowski committed
791 792 793 794
      localDofMap.createVec(fetiData.tmp_vec_b1, nGlobalOverallInterior);
      primalDofMap.createVec(fetiData.tmp_vec_primal1);
      interfaceDofMap.createVec(fetiData.tmp_vec_interface);

795 796 797 798 799
      MatCreateShell(mpiCommGlobal,
		     interfaceDofMap.getRankDofs() + lagrangeMap.getRankDofs(), 
		     interfaceDofMap.getRankDofs() + lagrangeMap.getRankDofs(),
		     interfaceDofMap.getOverallDofs() + lagrangeMap.getOverallDofs(), 
		     interfaceDofMap.getOverallDofs() + lagrangeMap.getOverallDofs(),
800
		     &fetiData, &mat_feti);
801
      MatShellSetOperation(mat_feti, MATOP_MULT, 
Thomas Witkowski's avatar
Thomas Witkowski committed
802
			   (void(*)(void))petscMultMatFetiInterface);