PetscSolverFeti.cc 51.2 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
    MSG("START A\n");

767
768
    // === Create FETI-DP solver object. ===

769
770
771
772
773
774
775
    FetiData &data = (!enableStokesMode ? fetiData : fetiDataInterface);
    data.mat_lagrange = &mat_lagrange;
    data.subSolver = subdomain;
    data.ksp_schur_primal = &ksp_schur_primal;   
    localDofMap.createVec(data.tmp_vec_b, nGlobalOverallInterior);
    lagrangeMap.createVec(data.tmp_vec_lagrange);
    primalDofMap.createVec(data.tmp_vec_primal);
776

777
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);
    } else {
      MSG("TEST A0\n");
788

789
790
791
792
793
794
      fetiDataInterface.mat_interior_interface = 
	getMatInteriorCoarseByComponent(2);
      fetiDataInterface.mat_interface_interior = 
	getMatCoarseInteriorByComponent(2);
      fetiDataInterface.mat_primal_interface = getMatCoarse(0, 1);
      fetiDataInterface.mat_interface_primal = getMatCoarse(1, 0);
795

796
797
798
799
800
801
802
803
804
805
806
807
808
      MSG("TEST A1\n");

      MatCreateShell(mpiCommGlobal,
		     interfaceDofMap.getRankDofs() + lagrangeMap.getRankDofs(), 
		     interfaceDofMap.getRankDofs() + lagrangeMap.getRankDofs(),
		     interfaceDofMap.getOverallDofs() + lagrangeMap.getOverallDofs(), 
		     interfaceDofMap.getOverallDofs() + lagrangeMap.getOverallDofs(),
		     &fetiDataInterface, &mat_feti);
      MatShellSetOperation(mat_feti, MATOP_MULT, 
			   (void(*)(void))petscMultMatFetiInterface);

      MSG("TEST A2\n");
    }
809

810
    KSPCreate(mpiCommGlobal, &ksp_feti);
811
    KSPSetOperators(ksp_feti, mat_feti, mat_feti, SAME_NONZERO_PATTERN);
812
813
814
    KSPSetOptionsPrefix(ksp_feti, "feti_");
    KSPSetType(ksp_feti, KSPGMRES);
    KSPSetTolerances(ksp_feti, 0, 1e-8, 1e+3, 1000);
815
    KSPSetFromOptions(ksp_feti);
816
817


818
    // === Create FETI-DP preconditioner object. ===
819

820
821
822
823
    if (fetiPreconditioner != FETI_NONE) {
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatScale(mat_lagrange_scaled, 0.5);
    }
824

825
    switch (fetiPreconditioner) {
826
827
828
829
    case FETI_DIRICHLET:      
      TEST_EXIT_DBG(meshLevel == 0)
	("Check for localDofMap.getLocalMatIndex, which should not work for multilevel FETI-DP!\n");

830
      KSPCreate(PETSC_COMM_SELF, &ksp_interior);
Thomas Witkowski's avatar
Thomas Witkowski committed
831
832
      KSPSetOperators(ksp_interior, mat_interior_interior, mat_interior_interior, 
		      SAME_NONZERO_PATTERN);
833
834
835
836
837
838
      KSPSetOptionsPrefix(ksp_interior, "precon_interior_");
      KSPSetType(ksp_interior, KSPPREONLY);
      PC pc_interior;
      KSPGetPC(ksp_interior, &pc_interior);
      PCSetType(pc_interior, PCLU);
      PCFactorSetMatSolverPackage(pc_interior, MATSOLVERUMFPACK);
839
840
841
842
843
844
845
846
847
      KSPSetFromOptions(ksp_interior);
            
      fetiDirichletPreconData.mat_lagrange_scaled = &mat_lagrange_scaled;
      fetiDirichletPreconData.mat_interior_interior = &mat_interior_interior;
      fetiDirichletPreconData.mat_duals_duals = &mat_duals_duals;
      fetiDirichletPreconData.mat_interior_duals = &mat_interior_duals;
      fetiDirichletPreconData.mat_duals_interior = &mat_duals_interior;
      fetiDirichletPreconData.ksp_interior = &ksp_interior;
      
848
849
      localDofMap.createVec(fetiDirichletPreconData.tmp_vec_b, 
			    nGlobalOverallInterior);
850
851
852
853
854
855
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiDirichletPreconData.tmp_vec_duals0));
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiDirichletPreconData.tmp_vec_duals1));
      MatGetVecs(mat_interior_interior, PETSC_NULL, 
		 &(fetiDirichletPreconData.tmp_vec_interior));
856

857
858
859
      TEST_EXIT_DBG(meshLevel == 0)
	("Should not happen, check usage of localDofMap!\n");

860
      for (unsigned int i = 0; i < feSpaces.size(); i++) {
861
	DofMap &dualMap = dualDofMap[feSpaces[i]].getMap();
862
	for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
863
	  DegreeOfFreedom d = it->first;	  
864
865
	  int matIndexLocal = localDofMap.getLocalMatIndex(i, d);
	  int matIndexDual = dualDofMap.getLocalMatIndex(i, d);
866
867
868
869
	  fetiDirichletPreconData.localToDualMap[matIndexLocal] = matIndexDual;
	}
      }

870
871
872
873
      KSPGetPC(ksp_feti, &precon_feti);
      PCSetType(precon_feti, PCSHELL);
      PCShellSetContext(precon_feti, static_cast<void*>(&fetiDirichletPreconData));
      PCShellSetApply(precon_feti, petscApplyFetiDirichletPrecon);
Thomas Witkowski's avatar
Thomas Witkowski committed
874
875
876
877
878
879
880
881


      // For the case, that we want to print the timings, we force the LU 
      // factorization of the local problems to be done here explicitly.
      if (printTimings) {
	double wtime = MPI::Wtime();
	KSPSetUp(ksp_interior);
	KSPSetUpOnBlocks(ksp_interior);
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
882
	MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
883
884
885
	MSG("FETI-DP timing 08: %.5f seconds (factorization of Dirichlet preconditoner matrices)\n",
	    MPI::Wtime() - wtime);
      }
886
887
888
889
890
891
892
      
      break;

    case FETI_LUMPED:
      fetiLumpedPreconData.mat_lagrange_scaled = &mat_lagrange_scaled;
      fetiLumpedPreconData.mat_duals_duals = &mat_duals_duals;

Thomas Witkowski's avatar
Thomas Witkowski committed
893
      for (unsigned int i = 0; i < feSpaces.size(); i++) {
894
	DofMap &dualMap = dualDofMap[feSpaces[i]].getMap();
895
	for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
Thomas Witkowski's avatar
Thomas Witkowski committed
896
	  DegreeOfFreedom d = it->first;
897
898
	  int matIndexLocal = localDofMap.getLocalMatIndex(i, d);
	  int matIndexDual = dualDofMap.getLocalMatIndex(i, d);
Thomas Witkowski's avatar
Thomas Witkowski committed
899
900
901
902
	  fetiLumpedPreconData.localToDualMap[matIndexLocal] = matIndexDual;
	}
      }