PetscSolverFeti.cc 51.9 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 86 87 88

    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");
      fullInterface = feSpaces[2];
    }
			   
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_NONE) {
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_NONE)
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 147
    if (fetiPreconditioner != FETI_NONE)
      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_NONE)
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 243 244 245 246 247 248

    if (enableStokesMode) {
      MSG("WARNING: MAKE THIS MORE GENERAL!!!!!\n");
      subdomain->setCoarseSpaceDofMapping(&primalDofMap, 0);
      subdomain->setCoarseSpaceDofMapping(&primalDofMap, 1);
      subdomain->setCoarseSpaceDofMapping(&interfaceDofMap, 2);
    } else {	  
      subdomain->setCoarseSpaceDofMapping(&primalDofMap); 
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
249 250

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


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

263 264 265
    if (feSpace == fullInterface)
      return;

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

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

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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
283 284 285 286 287 288 289
	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);
	}
      }
290
    }
291

292 293 294 295

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

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


304
  void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace)
305 306
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
307

308 309 310
    if (feSpace == fullInterface)
      return;

311 312 313
    // === Create global index of the dual nodes on each rank. ===

    DofContainer allBoundaryDofs;
314
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
315 316 317 318

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

  
328 329 330 331 332 333 334 335 336 337 338 339
  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)
340
      if (meshDistributor->getDofMap()[feSpace].isRankDof(**it))
Thomas Witkowski's avatar
Thomas Witkowski committed
341
	interfaceDofMap[feSpace].insertRankDof(**it);
342
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
343
	interfaceDofMap[feSpace].insertNonRankDof(**it);
344 345 346
  }


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

351 352 353
    if (feSpace == fullInterface)
      return;

354 355
    boundaryDofRanks[feSpace].clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
356 357 358
    // 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
359
    map<int, std::set<DegreeOfFreedom> > sdRankDofs;
360

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
361
    if (meshLevel > 0) {
362
      StdMpi<vector<int> > stdMpi(mpiCommGlobal);
363 364 365 366 367 368 369 370 371

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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
397 398 399 400 401 402 403
    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)).         ===

404
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
405 406 407 408 409 410 411
    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);

412
 	  if (meshLevel == 0 ||
Thomas Witkowski's avatar
Thomas Witkowski committed
413
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
414
	    boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());	  
415 416 417 418
	}
      }
    }

419 420 421 422

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

423
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
424

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

    stdMpi.updateSendDataSize();

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

      if (recvFromRank)
450
	stdMpi.recv(it.getRank());
451
    }
452

453 454
    stdMpi.startCommunication();

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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
469

470 471 472
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

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


488
  void PetscSolverFeti::createIndexB(const FiniteElemSpace *feSpace)
489
  {
490
    FUNCNAME("PetscSolverFeti::createIndexB()");
491

492
    DOFAdmin* admin = feSpace->getAdmin();
493 494 495 496

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

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

507 508 509 510 511 512 513 514 515
	  if (fetiPreconditioner != FETI_NONE)
	    interiorDofMap[feSpace].insertRankDof(i, nLocalInterior);
	  
	  nLocalInterior++;	  
	} else {
	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(i))
	    localDofMap[feSpace].insertRankDof(i);
	  else
	    localDofMap[feSpace].insertNonRankDof(i);
516 517 518

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

525 526 527 528 529 530 531 532 533 534
    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);
      }
535 536 537
  }


538
  void PetscSolverFeti::createMatLagrange(vector<const FiniteElemSpace*> &feSpaces)
539 540 541
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
542
    double wtime = MPI::Wtime();
543
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Thomas Witkowski committed
544

545 546
    // === Create distributed matrix for Lagrange constraints. ===

547 548 549 550 551
    MatCreateAIJ(mpiCommGlobal,
		 lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		 lagrangeMap.getOverallDofs(), nGlobalOverallInterior,
		 2, PETSC_NULL, 2, PETSC_NULL,
		 &mat_lagrange);
552
    MatSetOption(mat_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
553

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

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

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

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

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

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

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


605
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
606
  {
607
    FUNCNAME("PetscSolverFeti::createSchurPrimalKsp()");
608

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

Thomas Witkowski's avatar
Thomas Witkowski committed
612
      schurPrimalData.subSolver = subdomain;
Thomas Witkowski's avatar
Thomas Witkowski committed
613
      
614 615
      localDofMap.createVec(schurPrimalData.tmp_vec_b, nGlobalOverallInterior);
      primalDofMap.createVec(schurPrimalData.tmp_vec_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
616

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

635 636
      double wtime = MPI::Wtime();

637 638 639
      TEST_EXIT_DBG(meshLevel == 0)
	("Does not support for multilevel, check usage of localDofMap.\n");

640 641 642
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
      int nRowsRankB = localDofMap.getRankDofs();
643

Thomas Witkowski's avatar
Thomas Witkowski committed
644
      Mat matBPi;
645 646 647
      MatCreateAIJ(mpiCommGlobal,
		   nRowsRankB, nRowsRankPrimal, 
		   nGlobalOverallInterior, nRowsOverallPrimal,
648 649 650
		   150, PETSC_NULL, 150, PETSC_NULL, &matBPi);
      MatSetOption(matBPi, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

651

652 653 654
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
655

656 657
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

658
      for (int i = 0; i < nRowsRankB; i++) {
659
	PetscInt row = localDofMap.getStartDofs() + i;
660
	MatGetRow(subdomain->getMatInteriorCoarse(), row, &nCols, &cols, &values);
661 662 663 664 665

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

666
	MatRestoreRow(subdomain->getMatInteriorCoarse(), row, &nCols, &cols, &values);
667 668
      }

669
      TEST_EXIT(static_cast<int>(mat_b_primal_cols.size()) == 
670
		primalDofMap.getLocalDofs())
671
	("Should not happen!\n");
672

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

 	for (unsigned int i = 0; i < it->second.size(); i++) 
 	  VecSetValue(tmpVec, 
 		      it->second[i].first, it->second[i].second, INSERT_VALUES);

	VecAssemblyBegin(tmpVec);
	VecAssemblyEnd(tmpVec);

Thomas Witkowski's avatar
Thomas Witkowski committed
685
	subdomain->solve(tmpVec, tmpVec);
686

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

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
700 701
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
702

703 704
      MatDuplicate(subdomain->getMatCoarse(), MAT_COPY_VALUES, 
		   &mat_schur_primal);
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
705 706

      Mat matPrimal;
707 708
      MatMatMult(subdomain->getMatCoarseInterior(), matBPi, MAT_INITIAL_MATRIX, 
		 PETSC_DEFAULT, &matPrimal);
709
      MatAXPY(mat_schur_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);
710 711

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
712
      MatDestroy(&matBPi);
713

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

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


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

749
    if (schurPrimalSolver == 0) {
750
      schurPrimalData.subSolver = NULL;
751

752 753 754
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
    }
755 756 757
    
    MatDestroy(&mat_schur_primal);
    KSPDestroy(&ksp_schur_primal);
758 759 760
  }


761
  void PetscSolverFeti::createFetiKsp(vector<const FiniteElemSpace*> &feSpaces)
762 763 764
  {
    FUNCNAME("PetscSolverFeti::createFetiKsp()");

765 766
    // === Create FETI-DP solver object. ===

767 768 769 770
    fetiData.mat_lagrange = &mat_lagrange;
    fetiData.subSolver = subdomain;
    fetiData.ksp_schur_primal = &ksp_schur_primal;

Thomas Witkowski's avatar
Thomas Witkowski committed
771
    localDofMap.createVec(fetiData.tmp_vec_b0, nGlobalOverallInterior);
772
    lagrangeMap.createVec(fetiData.tmp_vec_lagrange);
Thomas Witkowski's avatar
Thomas Witkowski committed
773
    primalDofMap.createVec(fetiData.tmp_vec_primal0);
774

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

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

799
    KSPCreate(mpiCommGlobal, &ksp_feti);
800
    KSPSetOperators(ksp_feti, mat_feti, mat_feti, SAME_NONZERO_PATTERN);
801 802 803
    KSPSetOptionsPrefix(ksp_feti, "feti_");
    KSPSetType(ksp_feti, KSPGMRES);
    KSPSetTolerances(ksp_feti, 0, 1e-8, 1e+3, 1000);
804
    KSPSetFromOptions(ksp_feti);
805

Thomas Witkowski's avatar
Thomas Witkowski committed
806 807 808 809 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);
      MatSetNullSpace(mat_feti, matNullspace);
      KSPSetNullSpace(ksp_feti, matNullspace);
    }

827

828
    // === Create FETI-DP preconditioner object. ===