PetscSolverFeti.cc 48.5 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


13
#include "AMDiS.h"
14
#include "parallel/PetscSolverFeti.h"
15
#include "parallel/PetscSolverFetiStructs.h"
16
17
#include "parallel/StdMpi.h"
#include "parallel/MpiHelper.h"
18
#include "parallel/SubDomainSolver.h"
19
#include "io/VtkWriter.h"
20
21
22
23
24

namespace AMDiS {

  using namespace std;

25
26
27
28
29
30
31
  // y = mat * x
  int petscMultMatSchurPrimal(Mat mat, Vec x, Vec y)
  {
    // S_PiPi = K_PiPi - K_PiB inv(K_BB) K_BPi

    void *ctx;
    MatShellGetContext(mat, &ctx);
32
    SchurPrimalData* data = static_cast<SchurPrimalData*>(ctx);
33

34
35
36
37
    MatMult(data->subSolver->getMatIntCoarse(), x, data->tmp_vec_b);
    data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
    MatMult(data->subSolver->getMatCoarseInt(), data->tmp_vec_b, data->tmp_vec_primal);
    MatMult(data->subSolver->getMatCoarseCoarse(), x, y);
38
39
40
41
42
43
44
45
46
    VecAXPBY(y, -1.0, 1.0, data->tmp_vec_primal);

    return 0;
  }


  // y = mat * x
  int petscMultMatFeti(Mat mat, Vec x, Vec y)
  {
47
48
    //    F = L inv(K_BB) trans(L) + L inv(K_BB) K_BPi inv(S_PiPi) K_PiB inv(K_BB) trans(L)
    // => F = L [I + inv(K_BB) K_BPi inv(S_PiPi) K_PiB] inv(K_BB) trans(L)
49
50
51

    void *ctx;
    MatShellGetContext(mat, &ctx);
52
    FetiData* data = static_cast<FetiData*>(ctx);
53

Thomas Witkowski's avatar
Thomas Witkowski committed
54
    MatMultTranspose(*(data->mat_lagrange), x, data->tmp_vec_b);
55
    data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
Thomas Witkowski's avatar
Thomas Witkowski committed
56
    MatMult(*(data->mat_lagrange), data->tmp_vec_b, data->tmp_vec_lagrange);
Thomas Witkowski's avatar
Thomas Witkowski committed
57

58
    MatMult(data->subSolver->getMatCoarseInt(), data->tmp_vec_b, data->tmp_vec_primal);
59
    KSPSolve(*(data->ksp_schur_primal), data->tmp_vec_primal, data->tmp_vec_primal);
60
   MatMult(data->subSolver->getMatIntCoarse(), data->tmp_vec_primal, data->tmp_vec_b);
61
    data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
Thomas Witkowski's avatar
Thomas Witkowski committed
62
    MatMult(*(data->mat_lagrange), data->tmp_vec_b, y);
63

Thomas Witkowski's avatar
Thomas Witkowski committed
64
    VecAXPBY(y, 1.0, 1.0, data->tmp_vec_lagrange);
65
66
67
68
69

    return 0;
  }


70
  // y = PC * x
71
  PetscErrorCode petscApplyFetiDirichletPrecon(PC pc, Vec x, Vec y)
72
  {
73
    // Get data for the preconditioner
74
75
    void *ctx;
    PCShellGetContext(pc, &ctx);
76
    FetiDirichletPreconData* data = static_cast<FetiDirichletPreconData*>(ctx);
77

78
    // Multiply with scaled Lagrange constraint matrix.
79
80
81
    MatMultTranspose(*(data->mat_lagrange_scaled), x, data->tmp_vec_b);


82
    // === Restriction of the B nodes to the boundary nodes. ===
83

84
85
86
87
    int nLocalB;
    int nLocalDuals;
    VecGetLocalSize(data->tmp_vec_b, &nLocalB);
    VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals);
88

89
90
91
92
    PetscScalar *local_b, *local_duals;
    VecGetArray(data->tmp_vec_b, &local_b);
    VecGetArray(data->tmp_vec_duals0, &local_duals);

93
94
95
    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_duals[it->second] = local_b[it->first];
96
97

    VecRestoreArray(data->tmp_vec_b, &local_b);
98
    VecRestoreArray(data->tmp_vec_duals0, &local_duals);
99
100


101
    // === K_DD - K_DI inv(K_II) K_ID ===
102

103
    MatMult(*(data->mat_duals_duals), data->tmp_vec_duals0, data->tmp_vec_duals1);
104

105
    MatMult(*(data->mat_interior_duals), data->tmp_vec_duals0, data->tmp_vec_interior);
106
    KSPSolve(*(data->ksp_interior), data->tmp_vec_interior, data->tmp_vec_interior);
107
108
109
110
111
112
113
114
115
116
    MatMult(*(data->mat_duals_interior), data->tmp_vec_interior, data->tmp_vec_duals0);

    VecAXPBY(data->tmp_vec_duals0, 1.0, -1.0, data->tmp_vec_duals1);


    // === Prolongation from local dual nodes to B nodes.

    VecGetArray(data->tmp_vec_b, &local_b);
    VecGetArray(data->tmp_vec_duals0, &local_duals);

117
118
119
    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_b[it->first] = local_duals[it->second];
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141

    VecRestoreArray(data->tmp_vec_b, &local_b);
    VecRestoreArray(data->tmp_vec_duals0, &local_duals);


    // Multiply with scaled Lagrange constraint matrix.
    MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b, y);

    return 0;
  }


  // y = PC * x
  PetscErrorCode petscApplyFetiLumpedPrecon(PC pc, Vec x, Vec y)
  {
    // Get data for the preconditioner
    void *ctx;
    PCShellGetContext(pc, &ctx);
    FetiLumpedPreconData* data = static_cast<FetiLumpedPreconData*>(ctx);

    // Multiply with scaled Lagrange constraint matrix.
    MatMultTranspose(*(data->mat_lagrange_scaled), x, data->tmp_vec_b);
142
143


144
    // === Restriction of the B nodes to the boundary nodes. ===
145

146
147
148
149
    int nLocalB;
    int nLocalDuals;
    VecGetLocalSize(data->tmp_vec_b, &nLocalB);
    VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals);
150

151
    PetscScalar *local_b, *local_duals;
152
    VecGetArray(data->tmp_vec_b, &local_b);
153
    VecGetArray(data->tmp_vec_duals0, &local_duals);
154

Thomas Witkowski's avatar
Thomas Witkowski committed
155
156
157
    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_duals[it->second] = local_b[it->first];
158
159

    VecRestoreArray(data->tmp_vec_b, &local_b);
160
161
162
163
164
165
166
    VecRestoreArray(data->tmp_vec_duals0, &local_duals);


    // === K_DD ===

    MatMult(*(data->mat_duals_duals), data->tmp_vec_duals0, data->tmp_vec_duals1);

167

168
    // === Prolongation from local dual nodes to B nodes.
169

170
171
172
    VecGetArray(data->tmp_vec_b, &local_b);
    VecGetArray(data->tmp_vec_duals1, &local_duals);

Thomas Witkowski's avatar
Thomas Witkowski committed
173
174
175
    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_b[it->first] = local_duals[it->second];
176

177
178
    VecRestoreArray(data->tmp_vec_b, &local_b);
    VecRestoreArray(data->tmp_vec_duals0, &local_duals);
179

180
181

    // Multiply with scaled Lagrange constraint matrix.
182
183
184
185
186
187
    MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b, y);

    return 0;
  }


188
189
  PetscSolverFeti::PetscSolverFeti()
    : PetscSolver(),
190
      schurPrimalSolver(0),
191
      multiLevelTest(false),
192
      subDomainSolver(NULL),
193
194
195
      meshLevel(0),
      rStartInterior(0),
      nGlobalOverallInterior(0)
196
197
198
199
200
201
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

    string preconditionerName = "";
    Parameters::get("parallel->solver->precon", preconditionerName);
    if (preconditionerName == "" || preconditionerName == "none") {
202
      MSG("Create FETI-DP solver with no preconditioner!\n");
203
204
      fetiPreconditioner = FETI_NONE;
    } else if (preconditionerName == "dirichlet") {
205
      MSG("Create FETI-DP solver with Dirichlet preconditioner!\n");
206
207
      fetiPreconditioner = FETI_DIRICHLET;
    } else if (preconditionerName == "lumped") {
208
      MSG("Create FETI-DP solver with lumped preconditioner!\n");
209
210
      fetiPreconditioner = FETI_LUMPED;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
211
212
      ERROR_EXIT("Preconditioner \"%s\" not available!\n", 
		 preconditionerName.c_str());
213
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
214
215
216
217
218

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

    Parameters::get("parallel->multi level test", multiLevelTest);
221
222
    if (multiLevelTest)
      meshLevel = 1;
223
224
225
  }


226
  void PetscSolverFeti::initialize(vector<const FiniteElemSpace*> feSpaces)
227
  {
228
229
    FUNCNAME("PetscSolverFeti::initialize()");

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

233
234
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

235
236
237
238
239
240
241
242
243
    if (subDomainSolver == NULL) {
      if (meshLevel == 0) {
	subDomainSolver = 
	  new SubDomainSolver(meshDistributor, mpiComm, mpiSelfComm);
      } else {
	subDomainSolver = 
	  new SubDomainSolver(meshDistributor, 
			      levelData.getMpiComm(meshLevel - 1),
			      levelData.getMpiComm(meshLevel));
244
	subDomainSolver->setLevel(meshLevel);
245
246
      }
    }
247

248
249
    primalDofMap.init(levelData,
		      feSpaces, meshDistributor->getFeSpaces(), 
250
		      true, true);
251
252
    primalDofMap.setComputeMatIndex(true);
    
253
254
    dualDofMap.init(levelData, 
		    feSpaces, meshDistributor->getFeSpaces(),
255
		    false, false);
256
257
    dualDofMap.setComputeMatIndex(true);

258
259
    lagrangeMap.init(levelData, 
		     feSpaces, meshDistributor->getFeSpaces(),
260
		     true, true);
261
    lagrangeMap.setComputeMatIndex(true);
262
263
264
265
266
267
268
269
270
271
272
273

    if (meshLevel == 0) {
      localDofMap.init(levelData, 
		       feSpaces, meshDistributor->getFeSpaces(),
		       false, false);
      localDofMap.setComputeMatIndex(true);
    } else {
      localDofMap.init(levelData, 
		       feSpaces, meshDistributor->getFeSpaces(),
		       true, true);
      localDofMap.setComputeMatIndex(true);
    }
274
275

    if (fetiPreconditioner == FETI_DIRICHLET) {
276
277
278
      TEST_EXIT(meshLevel == 0)
	("Dirichlet preconditioner not yet implemented for multilevel FETI-DP\n");

279
280
      interiorDofMap.init(levelData, 
			  feSpaces,  meshDistributor->getFeSpaces(),
281
			  false, false);
282
283
      interiorDofMap.setComputeMatIndex(true);
    }
284
285
286
287
288
289
290
291
292
293
  }


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

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

295
296
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

297
298
299
300
    primalDofMap.clear();
    dualDofMap.clear();
    lagrangeMap.clear();
    localDofMap.clear();
301
302
303
    if (fetiPreconditioner == FETI_DIRICHLET)
      interiorDofMap.clear();

304
305
    primalDofMap.setDofComm(meshDistributor->getDofComm());
    lagrangeMap.setDofComm(meshDistributor->getDofComm());
306

307
308
309
310
    primalDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    dualDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    lagrangeMap.setMpiComm(levelData.getMpiComm(0), 0);    
    localDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
311
312
    if (fetiPreconditioner != FETI_NONE)
      interiorDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
313
314
315
316

    if (meshLevel > 0)
      localDofMap.setDofComm(meshDistributor->getDofCommSd());

317
318
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
319

320
321
      createPrimals(feSpace);  

322
      createDuals(feSpace);      
323

324
325
      createIndexB(feSpace);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
326
327
328
329

    primalDofMap.update();
    dualDofMap.update();
    localDofMap.update();
330
    if (fetiPreconditioner != FETI_NONE)
331
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
332

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
333

334
335
336
337
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
      createLagrange(feSpace);
    }
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
338
    
339
340
341
    lagrangeMap.update();


342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
    // === ===

    if (meshLevel == 0) {
      rStartInterior = 0;
      nGlobalOverallInterior = localDofMap.getOverallDofs();
    } else {
      MeshLevelData& levelData = meshDistributor->getMeshLevelData();

      MSG("RANK %d FROM %d\n",
	  levelData.getMpiComm(1).Get_rank(),
	  levelData.getMpiComm(1).Get_size());

      int groupRowsInterior = 0;
      if (levelData.getMpiComm(1).Get_rank() == 0)
	groupRowsInterior = localDofMap.getOverallDofs();

      mpi::getDofNumbering(mpiComm, groupRowsInterior,
			   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);

      MSG("COMM TEST FETI-DP: %d %d %d %d %d\n",
	  levelData.getMpiComm(1).Get_size(),
	  localDofMap.getRankDofs(),
	  localDofMap.getOverallDofs(),
	  nGlobalOverallInterior, rStartInterior);
    }

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

379
      MSG("  nRankPrimals = %d   nOverallPrimals = %d\n", 
380
381
	  primalDofMap[feSpace].nRankDofs, 
	  primalDofMap[feSpace].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
382
      
383
      MSG("  nRankDuals = %d  nOverallDuals = %d\n",
384
385
	  dualDofMap[feSpace].nRankDofs, 
	  dualDofMap[feSpace].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
386

387
      MSG("  nRankLagrange = %d  nOverallLagrange = %d\n",
388
389
	  lagrangeMap[feSpace].nRankDofs, 
	  lagrangeMap[feSpace].nOverallDofs);
390

391
      TEST_EXIT_DBG(localDofMap[feSpace].size() + primalDofMap[feSpace].size() == 
392
393
		    static_cast<unsigned int>(feSpace->getAdmin()->getUsedDofs()))
	("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
394
    }
395
396
397


    // If multi level test, inform sub domain solver about coarse space.
398
    subDomainSolver->setDofMapping(&primalDofMap, &localDofMap);
399
400
401
  }


402
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
403
  {
404
    FUNCNAME("PetscSolverFeti::createPrimals()");  
405

406
407
408
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

409
410
    /// Set of DOF indices that are considered to be primal variables.
    
411
    DofContainerSet& vertices = 
412
      meshDistributor->getBoundaryDofInfo(feSpace, meshLevel).geoDofs[VERTEX];
413
414

    DofIndexSet primals;
415
    for (DofContainerSet::iterator it = vertices.begin(); 
416
417
418
	 it != vertices.end(); ++it) {
      WorldVector<double> c;
      feSpace->getMesh()->getDofIndexCoords(*it, feSpace, c);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
419

420
      double e = 1e-8;
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
421
      if (meshLevel == 0) {
422
	primals.insert(**it);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
423
424
425
426
427
428
429
430
      } else {
	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);
	}
      }
431
    }
432

433
434
435
436

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

437
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
438
439
      if (meshDistributor->getDofMap()[feSpace].isRankDof(*it))
	primalDofMap[feSpace].insertRankDof(*it);
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
440
      else
441
  	primalDofMap[feSpace].insertNonRankDof(*it);
442
443
444
  }


445
  void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace)
446
447
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
448

449
450
451
    // === Create global index of the dual nodes on each rank. ===

    DofContainer allBoundaryDofs;
452
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
453
454
455
456

    for (DofContainer::iterator it = allBoundaryDofs.begin();
	 it != allBoundaryDofs.end(); ++it)
      if (!isPrimal(feSpace, **it))
457
458
459
460
461
462
	if (meshLevel == 0) {
	  dualDofMap[feSpace].insertRankDof(**it);
	} else {
	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(**it))
	    dualDofMap[feSpace].insertRankDof(**it);
	}	  
463
464
465
466
467
468
469
  }

  
  void PetscSolverFeti::createLagrange(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

470
471
    boundaryDofRanks[feSpace].clear();

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
472
473
474
475
    // Stores for all rank own communication DOFs, if the counterpart is
    // a rank owmed DOF in its subdomain. Thus, the following map stores to
    // each rank number all DOFS that fulfill this requirenment.
    map<int, std::set<DegreeOfFreedom> > sdRankDofs;
476

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
477
    if (meshLevel > 0) {
478
479
480
481
482
483
484
485
486
487
      StdMpi<vector<int> > stdMpi(mpiComm);

      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
488
	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))
489
	    subdomainRankDofs.push_back(1);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
490
	  else
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
	    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
506
507
508
509
510
511
	   !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());
    }
512
513


Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
    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)).         ===

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

	  if (meshLevel == 0 ||
	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
	    boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());
531
532
533
534
	}
      }
    }

535
536
537
538

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

539
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
540

541
    for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), meshLevel, feSpace);
542
543
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
544
	if (!isPrimal(feSpace, it.getDofIndex()))
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
545
546
	  if (meshLevel == 0 ||
	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
547
	  stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]);
548
549
550

    stdMpi.updateSendDataSize();

551
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
552
	 !it.end(); it.nextRank()) {
553
      bool recvFromRank = false;
554
      for (; !it.endDofIter(); it.nextDof()) {
555
	if (!isPrimal(feSpace, it.getDofIndex())) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
556
557
558
559
560
561
	  if (meshLevel == 0 ||
	      (meshLevel > 0 && 
	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))) {
	    recvFromRank = true;
	    break;
	  }
562
	}
563
      }
564
565

      if (recvFromRank)
566
	stdMpi.recv(it.getRank());
567
    }
568

569
570
    stdMpi.startCommunication();

571
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
572
	 !it.end(); it.nextRank()) {
573
      int i = 0;
574
      for (; !it.endDofIter(); it.nextDof())
575
	if (!isPrimal(feSpace, it.getDofIndex()))
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
576
577
578
579
580
	  if (meshLevel == 0 ||
	      (meshLevel > 0 && 
	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex())))	    
	    boundaryDofRanks[feSpace][it.getDofIndex()] = 
	      stdMpi.getRecvData(it.getRank())[i++];
581
582
    }

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
583

584
585
586
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

587
    int nRankLagrange = 0;
588
    DofMap& dualMap = dualDofMap[feSpace].getMap();
589
    for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
590
591
      if (meshDistributor->getDofMap()[feSpace].isRankDof(it->first)) {
	lagrangeMap[feSpace].insertRankDof(it->first, nRankLagrange);
592
	int degree = boundaryDofRanks[feSpace][it->first].size();
593
	nRankLagrange += (degree * (degree - 1)) / 2;
Thomas Witkowski's avatar
Thomas Witkowski committed
594
      } else {
595
	lagrangeMap[feSpace].insertNonRankDof(it->first);
596
597
      }
    }
598
    lagrangeMap[feSpace].nRankDofs = nRankLagrange;
599
600
601
  }


602
  void PetscSolverFeti::createIndexB(const FiniteElemSpace *feSpace)
603
  {
604
    FUNCNAME("PetscSolverFeti::createIndexB()");
605

606
    DOFAdmin* admin = feSpace->getAdmin();
607
608
609
610

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

612
    int nLocalInterior = 0;    
613
    for (int i = 0; i < admin->getUsedSize(); i++) {
614
      if (admin->isDofFree(i) == false && 
615
	  isPrimal(feSpace, i) == false &&
616
	  isDual(feSpace, i) == false) {
617
618
	if (meshLevel == 0) {
	  localDofMap[feSpace].insertRankDof(i, nLocalInterior);
619

620
621
622
623
624
625
626
627
628
	  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);
629
630
631

	  TEST_EXIT_DBG(fetiPreconditioner == FETI_NONE)
	    ("Not yet implemnted!\n");
632
	}
633
      }
634
    }
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
635
    
636
637
    // === And finally, add the global indicies of all dual nodes. ===

638
639
640
641
642
643
644
645
646
647
    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);
      }
648
649
650
  }


651
  void PetscSolverFeti::createMatLagrange(vector<const FiniteElemSpace*> &feSpaces)
652
653
654
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

655
656
    // === Create distributed matrix for Lagrange constraints. ===

657
    MatCreateMPIAIJ(mpiComm,
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
658
659
		    lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		    lagrangeMap.getOverallDofs(), nGlobalOverallInterior,
660
661
662
		    2, PETSC_NULL, 2, PETSC_NULL,
		    &mat_lagrange);

663
664
665
666
667
668
669
    // === 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
670
    for (unsigned int k = 0; k < feSpaces.size(); k++) {
671
      DofMap &dualMap = dualDofMap[feSpaces[k]].getMap();
672

673
      for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
674
	TEST_EXIT_DBG(boundaryDofRanks[feSpaces[k]].count(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
675
676
677
	  ("Should not happen!\n");
	
	// Global index of the first Lagrange constriant for this node.
678
	int index = lagrangeMap.getMatIndex(k, it->first);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
679

Thomas Witkowski's avatar
Thomas Witkowski committed
680
	// Copy set of all ranks that contain this dual node.
681
682
	vector<int> W(boundaryDofRanks[feSpaces[k]][it->first].begin(), 
		      boundaryDofRanks[feSpaces[k]][it->first].end());
Thomas Witkowski's avatar
Thomas Witkowski committed
683
684
	// Number of ranks that contain this dual node.
	int degree = W.size();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
685
686

	TEST_EXIT_DBG(degree > 1)("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
687
688
689
690
	
	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
691
692
	      int colIndex = 
		localDofMap.getMatIndex(k, it->first) + rStartInterior;
693

694
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
Thomas Witkowski's avatar
Thomas Witkowski committed
695
	      MSG("SET VALUE: %f\n", value);
696
	      MatSetValue(mat_lagrange, index, colIndex, value, INSERT_VALUES);
697
	    }
Thomas Witkowski's avatar
Thomas Witkowski committed
698
	    index++;	      
699
700
701
702
703
704
705
706
707
708
	  }
	}
      }
    }

    MatAssemblyBegin(mat_lagrange, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(mat_lagrange, MAT_FINAL_ASSEMBLY);
  }


709
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
710
  {
711
    FUNCNAME("PetscSolverFeti::createSchurPrimalKsp()");
712

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

716
      schurPrimalData.subSolver = subDomainSolver;
Thomas Witkowski's avatar
Thomas Witkowski committed
717
      
718
      VecCreateMPI(mpiComm, 
719
		   localDofMap.getRankDofs(), 
720
		   nGlobalOverallInterior,
Thomas Witkowski's avatar
Thomas Witkowski committed
721
		   &(schurPrimalData.tmp_vec_b));
722
      VecCreateMPI(mpiComm, 
723
724
		   primalDofMap.getRankDofs(), 
		   primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
725
726
		   &(schurPrimalData.tmp_vec_primal));

727
      MatCreateShell(mpiComm,
728
729
730
731
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getOverallDofs(), 
		     primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
732
733
734
735
736
		     &schurPrimalData, 
		     &mat_schur_primal);
      MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			   (void(*)(void))petscMultMatSchurPrimal);
      
737
      KSPCreate(mpiComm, &ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
738
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
739
740
      KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
      KSPSetType(ksp_schur_primal, KSPGMRES);
Thomas Witkowski's avatar
Thomas Witkowski committed
741
742
      KSPSetFromOptions(ksp_schur_primal);
    } else {
743
744
      MSG("Create direct schur primal solver!\n");

745
746
      double wtime = MPI::Wtime();

747
748
749
      TEST_EXIT_DBG(meshLevel == 0)
	("Does not support for multilevel, check usage of localDofMap.\n");

750
751
752
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
      int nRowsRankB = localDofMap.getRankDofs();
753

Thomas Witkowski's avatar
Thomas Witkowski committed
754
      Mat matBPi;
755
      MatCreateMPIAIJ(mpiComm,
Thomas Witkowski's avatar
Thomas Witkowski committed
756
		      nRowsRankB, nRowsRankPrimal, 
757
		      nGlobalOverallInterior, nRowsOverallPrimal,
Thomas Witkowski's avatar
Thomas Witkowski committed
758
759
		      30, PETSC_NULL, 30, PETSC_NULL, &matBPi);
      Mat matPrimal;
760

761
762
763
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
764

765
766
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

767
      for (int i = 0; i < nRowsRankB; i++) {
768
	PetscInt row = localDofMap.getStartDofs() + i;
769
	MatGetRow(subDomainSolver->getMatIntCoarse(), row, &nCols, &cols, &values);
770
771
772
773
774

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

775
	MatRestoreRow(subDomainSolver->getMatIntCoarse(), row, &nCols, &cols, &values);
776
777
      }

778
      TEST_EXIT(static_cast<int>(mat_b_primal_cols.size()) == 
779
		primalDofMap.getLocalDofs())
780
	("Should not happen!\n");
781

782
783
784
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
785
	VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);
786
787
788
789
790
791
792
793

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

794
	subDomainSolver->solve(tmpVec, tmpVec);
795

Thomas Witkowski's avatar
Thomas Witkowski committed
796
	PetscScalar *tmpValues;
797
	VecGetArray(tmpVec, &tmpValues);
798
	for (int i  = 0; i < nRowsRankB; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
799
	  MatSetValue(matBPi, 
800
		      localDofMap.getStartDofs() + i,
Thomas Witkowski's avatar
Thomas Witkowski committed
801
802
803
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
804
805
806
807
808
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
809
810
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
811
812

      MatDuplicate(subDomainSolver->getMatCoarseCoarse(), MAT_COPY_VALUES, &mat_schur_primal);
813
      MatMatMult(subDomainSolver->getMatCoarseInt(), matBPi, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &matPrimal);
814
      MatAXPY(mat_schur_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);
815
816

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
817
      MatDestroy(&matBPi);
818
819

      MatInfo minfo;
820
      MatGetInfo(mat_schur_primal, MAT_GLOBAL_SUM, &minfo);
821
822
      MSG("Schur primal matrix nnz = %f\n", minfo.nz_used);

823
      KSPCreate(mpiComm, &ksp_schur_primal);
824
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal,
825
		      SAME_NONZERO_PATTERN);
826
827
828
829
830
831
      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);
832
      KSPSetFromOptions(ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
833

834
835
      MSG("Creating Schur primal matrix needed %.5f seconds.\n",
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
Thomas Witkowski committed
836
    }
837
838
839
840
841
842
843
  }


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

844
    if (schurPrimalSolver == 0) {
845
      schurPrimalData.subSolver = NULL;
846

847
848
849
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
    }
850
851
852
    
    MatDestroy(&mat_schur_primal);
    KSPDestroy(&ksp_schur_primal);
853
854
855
  }


856
  void PetscSolverFeti::createFetiKsp(vector<const FiniteElemSpace*> &feSpaces)
857
858
859
  {
    FUNCNAME("PetscSolverFeti::createFetiKsp()");

860
861
    // === Create FETI-DP solver object. ===

862
    fetiData.mat_lagrange = &mat_lagrange;
863
    fetiData.subSolver = subDomainSolver;
864
    fetiData.ksp_schur_primal = &ksp_schur_primal;
865

866
    VecCreateMPI(mpiComm, 
867
		 localDofMap.getRankDofs(), 
868
		 nGlobalOverallInterior,
Thomas Witkowski's avatar
Thomas Witkowski committed
869
		 &(fetiData.tmp_vec_b));
870
    VecCreateMPI(mpiComm,
871
872
		 lagrangeMap.getRankDofs(), 
		 lagrangeMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
873
		 &(fetiData.tmp_vec_lagrange));
874
    VecCreateMPI(mpiComm, 
875
876
		 primalDofMap.getRankDofs(), 
		 primalDofMap.getOverallDofs(),
877
		 &(fetiData.tmp_vec_primal));
878

879
    MatCreateShell(mpiComm,
880
881
882
883
		   lagrangeMap.getRankDofs(), 
		   lagrangeMap.getRankDofs(),
		   lagrangeMap.getOverallDofs(), 
		   lagrangeMap.getOverallDofs(),
884
		   &fetiData, &mat_feti);
885
886
887
    MatShellSetOperation(mat_feti, MATOP_MULT, (void(*)(void))petscMultMatFeti);


888
    KSPCreate(mpiComm, &ksp_feti);
889
    KSPSetOperators(ksp_feti, mat_feti, mat_feti, SAME_NONZERO_PATTERN);
890
891
892
    KSPSetOptionsPrefix(ksp_feti, "feti_");
    KSPSetType(ksp_feti, KSPGMRES);
    KSPSetTolerances(ksp_feti, 0, 1e-8, 1e+3, 1000);
893
    KSPSetFromOptions(ksp_feti);
894
895


896
    // === Create FETI-DP preconditioner object. ===
897

898
899
900
901
    if (fetiPreconditioner != FETI_NONE) {
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatScale(mat_lagrange_scaled, 0.5);
    }
902

903
    switch (fetiPreconditioner) {
904
905
906
907
    case FETI_DIRICHLET:      
      TEST_EXIT_DBG(meshLevel == 0)
	("Check for localDofMap.getLocalMatIndex, which should not work for multilevel FETI-DP!\n");

908
      KSPCreate(PETSC_COMM_SELF, &ksp_interior);
Thomas Witkowski's avatar
Thomas Witkowski committed
909
910
      KSPSetOperators(ksp_interior, mat_interior_interior, mat_interior_interior, 
		      SAME_NONZERO_PATTERN);
911
912
913
914
915
916
      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);
917
918
919
920
921
922
923
924
925
      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;
      
926
      VecCreateMPI(mpiComm, 
927
		   localDofMap.getRankDofs(),
928
		   nGlobalOverallInterior,
929
		   &(fetiDirichletPreconData.tmp_vec_b));      
930
931
932
933
934
935
      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));
936

937
938
939
      TEST_EXIT_DBG(meshLevel == 0)
	("Should not happen, check usage of localDofMap!\n");

940
      for (unsigned int i = 0; i < feSpaces.size(); i++) {
941
	DofMap &dualMap = dualDofMap[feSpaces[i]].getMap();
942
	for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
943
	  DegreeOfFreedom d = it->first;	  
944
945
	  int matIndexLocal = localDofMap.getLocalMatIndex(i, d);
	  int matIndexDual = dualDofMap.getLocalMatIndex(i, d);
946
947
948
949
	  fetiDirichletPreconData.localToDualMap[matIndexLocal] = matIndexDual;
	}
      }

Thomas Witkowski's avatar