PetscSolverFeti.cc 59.4 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 37
      printTimings(false),
      enableStokesMode(false)
38 39 40 41 42 43
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

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

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

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

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


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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
80 81 82 83 84 85

    enableStokesMode = false;
    Parameters::get("parallel->stokes mode", enableStokesMode);
    if (enableStokesMode) {
      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");
86
      fullInterface = feSpaces[0];
Thomas Witkowski's avatar
Thomas Witkowski committed
87 88
    }
			   
Thomas Witkowski's avatar
Thomas Witkowski committed
89 90
    if (subdomain == NULL) {
      subdomain = new PetscSolverGlobalMatrix();
91

92
      if (meshLevel == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
93
	subdomain->setMeshDistributor(meshDistributor, 
94
				      mpiCommGlobal, mpiCommLocal);
95
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
96
	subdomain->setMeshDistributor(meshDistributor, 
97 98
				      levelData.getMpiComm(meshLevel - 1),
				      levelData.getMpiComm(meshLevel));
Thomas Witkowski's avatar
Thomas Witkowski committed
99
	subdomain->setLevel(meshLevel);
100 101
      }
    }
102

103 104 105 106
    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
107

108
    if (fullInterface != NULL)
Thomas Witkowski's avatar
Thomas Witkowski committed
109
      interfaceDofMap.init(levelData, feSpaces, uniqueFe);
Thomas Witkowski's avatar
Thomas Witkowski committed
110

111
    if (fetiPreconditioner == FETI_DIRICHLET) {
112 113 114
      TEST_EXIT(meshLevel == 0)
	("Dirichlet preconditioner not yet implemented for multilevel FETI-DP\n");

115
      interiorDofMap.init(levelData, feSpaces, uniqueFe, false);
116
    }
117 118 119 120 121 122 123
  }


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

Thomas Witkowski's avatar
Thomas Witkowski committed
124 125
    double timeCounter = MPI::Wtime();

126 127 128
    TEST_EXIT(meshDistributor)("No mesh distributor object defined!\n");
    TEST_EXIT(meshDistributor->getFeSpaces().size() > 0)
      ("No FE space defined in mesh distributor!\n");
129

130 131
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

132 133 134 135
    primalDofMap.clear();
    dualDofMap.clear();
    lagrangeMap.clear();
    localDofMap.clear();
136
    if (fetiPreconditioner == FETI_DIRICHLET)
137 138
      interiorDofMap.clear();

139 140
    primalDofMap.setDofComm(meshDistributor->getDofComm());
    lagrangeMap.setDofComm(meshDistributor->getDofComm());
141

142 143
    primalDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    dualDofMap.setMpiComm(levelData.getMpiComm(0), 0);
144
    lagrangeMap.setMpiComm(levelData.getMpiComm(0), 0);
145
    localDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
146
    if (fetiPreconditioner == FETI_DIRICHLET)
147
      interiorDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
148

149 150 151
    if (meshLevel == 0)
      localDofMap.setDofComm(meshDistributor->getDofComm());
    else
152 153
      localDofMap.setDofComm(meshDistributor->getDofCommSd());

154 155 156 157 158 159
    if (fullInterface != NULL) {
      interfaceDofMap.clear();
      interfaceDofMap.setDofComm(meshDistributor->getDofComm());
      interfaceDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    }

160 161
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
162
    
163 164
      createPrimals(feSpace);  

165 166 167
      createDuals(feSpace);

      createInterfaceNodes(feSpace);
168

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
169
      createIndexB(feSpace);     
170
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
171 172 173 174

    primalDofMap.update();
    dualDofMap.update();
    localDofMap.update();
175
    if (fetiPreconditioner == FETI_DIRICHLET)
176
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
177

178 179 180
    if (fullInterface != NULL)
      interfaceDofMap.update();

181 182 183 184
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
      createLagrange(feSpace);
    }
185

186 187 188
    lagrangeMap.update();


189 190 191 192 193 194 195 196 197 198 199 200
    // === ===

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

201
      mpi::getDofNumbering(mpiCommGlobal, groupRowsInterior,
202 203 204 205 206 207 208 209 210
			   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);
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236
      if (feSpace == fullInterface) {
	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
237
    }
238

Thomas Witkowski's avatar
Thomas Witkowski committed
239
    subdomain->setDofMapping(&localDofMap);
Thomas Witkowski's avatar
Thomas Witkowski committed
240 241 242

    if (enableStokesMode) {
      MSG("WARNING: MAKE THIS MORE GENERAL!!!!!\n");
Thomas Witkowski's avatar
So  
Thomas Witkowski committed
243 244
      subdomain->setCoarseSpaceDofMapping(&interfaceDofMap, 0);
      //      subdomain->setCoarseSpaceDofMapping(&primalDofMap, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
245
      subdomain->setCoarseSpaceDofMapping(&primalDofMap, 1);
Thomas Witkowski's avatar
So  
Thomas Witkowski committed
246 247 248 249
      subdomain->setCoarseSpaceDofMapping(&primalDofMap, 2);
      subdomain->setCoarseSpaceDofMapping(&primalDofMap, 3);
      subdomain->setCoarseSpaceDofMapping(&primalDofMap, 4);

Thomas Witkowski's avatar
Thomas Witkowski committed
250 251 252
    } else {	  
      subdomain->setCoarseSpaceDofMapping(&primalDofMap); 
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
253 254

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


263
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
264
  {
265
    FUNCNAME("PetscSolverFeti::createPrimals()");  
266

267 268 269
    if (feSpace == fullInterface)
      return;

270 271 272
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

Thomas Witkowski's avatar
Thomas Witkowski committed
273
    // Set of DOF indices that are considered to be primal variables.
274
    DofContainerSet& vertices = 
275
      meshDistributor->getBoundaryDofInfo(feSpace, meshLevel).geoDofs[VERTEX];
276 277

    DofIndexSet primals;
278
    for (DofContainerSet::iterator it = vertices.begin(); 
279
	 it != vertices.end(); ++it) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
280
      if (meshLevel == 0) {
281
	primals.insert(**it);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
282
      } else {
283
	double e = 1e-8;
Thomas Witkowski's avatar
Thomas Witkowski committed
284 285 286
	WorldVector<double> c;
	feSpace->getMesh()->getDofIndexCoords(*it, feSpace, c);

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

296 297 298 299

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

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


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

312 313 314
    if (feSpace == fullInterface)
      return;

315 316 317
    // === Create global index of the dual nodes on each rank. ===

    DofContainer allBoundaryDofs;
318
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
319 320 321 322

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

  
332 333 334 335 336 337 338 339 340 341 342 343
  void PetscSolverFeti::createInterfaceNodes(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createInterfaceNodes()");

    if (feSpace != fullInterface)
      return;

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

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


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

355 356 357
    if (feSpace == fullInterface)
      return;

358 359
    boundaryDofRanks[feSpace].clear();

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

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

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

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

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

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

423 424 425 426

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

427
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
428

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

    stdMpi.updateSendDataSize();

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

      if (recvFromRank)
454
	stdMpi.recv(it.getRank());
455
    }
456

457 458
    stdMpi.startCommunication();

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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
473

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

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


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

496
    DOFAdmin* admin = feSpace->getAdmin();
497 498 499 500

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

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

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

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

529 530 531 532 533 534 535 536 537 538
    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);
      }
539 540 541
  }


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

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

549 550
    // === Create distributed matrix for Lagrange constraints. ===

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

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

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

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

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

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

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

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


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

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

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

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

639 640
      double wtime = MPI::Wtime();

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

644 645 646
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
      int nRowsRankB = localDofMap.getRankDofs();
647

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

655

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

660 661
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

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

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

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

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

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

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

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

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
704 705
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
706

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

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

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
716
      MatDestroy(&matBPi);
717

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

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


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

753
    if (schurPrimalSolver == 0) {
754
      schurPrimalData.subSolver = NULL;
755

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


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

769 770
    // === Create FETI-DP solver object. ===

771 772 773 774
    fetiData.mat_lagrange = &mat_lagrange;
    fetiData.subSolver = subdomain;
    fetiData.ksp_schur_primal = &ksp_schur_primal;

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

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

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

803
    KSPCreate(mpiCommGlobal, &ksp_feti);
804
    KSPSetOperators(ksp_feti, mat_feti, mat_feti, SAME_NONZERO_PATTERN);
805 806 807
    KSPSetOptionsPrefix(ksp_feti, "feti_");
    KSPSetType(ksp_feti, KSPGMRES);
    KSPSetTolerances(ksp_feti, 0, 1e-8, 1e+3, 1000);
808
    KSPSetFromOptions(ksp_feti);
809

Thomas Witkowski's avatar
Thomas Witkowski committed
810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826
    if (enableStokesMode) {
      Vec nullSpaceInterface;
      Vec nullSpaceLagrange;
      Vec nullSpaceBasis;

      interfaceDofMap.createVec(nullSpaceInterface);
      VecSet(nullSpaceInterface, 1.0);

      lagrangeMap.createVec(nullSpaceLagrange);
      VecSet(nullSpaceLagrange, 0.0);

      Vec vecArray[2] = {nullSpaceInterface, nullSpaceLagrange};
      VecCreateNest(mpiCommGlobal, 2, PETSC_NULL, vecArray, &nullSpaceBasis);

      
      MatNullSpace matNullspace;
      MatNullSpaceCreate(mpiCommGlobal, PETSC_FALSE, 1, &nullSpaceBasis, &matNullspace);
827 828
      MatSetNullSpace(mat_feti, matNullspace);
      KSPSetNullSpace(ksp_feti, matNullspace);
Thomas Witkowski's avatar
Thomas Witkowski committed
829 830
    }

831

832
    // === Create FETI-DP preconditioner object. ===
833

834 835 836 837
    if (fetiPreconditioner != FETI_NONE) {
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatScale(mat_lagrange_scaled, 0.5);
    }
838

839

840
    switch (fetiPreconditioner) {
841 842
    case FETI_DIRICHLET:
      TEST_EXIT(meshLevel == 0)
843 844
	("Check for localDofMap.getLocalMatIndex, which should not work for multilevel FETI-DP!\n");