PetscSolverFeti.cc 55 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/PetscSolverGlobalMatrix.h"
19
#include "io/VtkWriter.h"
20
21
22
23
24

namespace AMDiS {

  using namespace std;

25
26
27
28
29
30
  double FetiTimings::fetiSolve = 0.0;
  double FetiTimings::fetiSolve01 = 0.0;
  double FetiTimings::fetiSolve02 = 0.0;
  double FetiTimings::fetiPreconditioner = 0.0;


31
32
33
34
35
36
37
  // 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);
38
    SchurPrimalData* data = static_cast<SchurPrimalData*>(ctx);
39

40
41
42
43
    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);
44
45
46
47
48
49
50
51
52
    VecAXPBY(y, -1.0, 1.0, data->tmp_vec_primal);

    return 0;
  }


  // y = mat * x
  int petscMultMatFeti(Mat mat, Vec x, Vec y)
  {
53
54
    FUNCNAME("petscMultMatFeti()");

55
56
    //    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)
57

58
59
    double wtime = MPI::Wtime();

60
61
    void *ctx;
    MatShellGetContext(mat, &ctx);
62
    FetiData* data = static_cast<FetiData*>(ctx);
63

Thomas Witkowski's avatar
Thomas Witkowski committed
64
    MatMultTranspose(*(data->mat_lagrange), x, data->tmp_vec_b);
65
66

    double wtime01 = MPI::Wtime();
67
    data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
68
69
70

    FetiTimings::fetiSolve01 += (MPI::Wtime() - wtime01);

Thomas Witkowski's avatar
Thomas Witkowski committed
71
    MatMult(*(data->mat_lagrange), data->tmp_vec_b, data->tmp_vec_lagrange);
Thomas Witkowski's avatar
Thomas Witkowski committed
72

73
    MatMult(data->subSolver->getMatCoarseInt(), data->tmp_vec_b, data->tmp_vec_primal);
74
75

    wtime01 = MPI::Wtime();
76
    KSPSolve(*(data->ksp_schur_primal), data->tmp_vec_primal, data->tmp_vec_primal);
77
78
    FetiTimings::fetiSolve02 += (MPI::Wtime() - wtime01);

79
    MatMult(data->subSolver->getMatIntCoarse(), data->tmp_vec_primal, data->tmp_vec_b);
80
81

    wtime01 = MPI::Wtime();
82
    data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
83
84
    FetiTimings::fetiSolve01 += (MPI::Wtime() - wtime01);

Thomas Witkowski's avatar
Thomas Witkowski committed
85
    MatMult(*(data->mat_lagrange), data->tmp_vec_b, y);
86

Thomas Witkowski's avatar
Thomas Witkowski committed
87
    VecAXPBY(y, 1.0, 1.0, data->tmp_vec_lagrange);
88

89
90
    FetiTimings::fetiSolve += (MPI::Wtime() - wtime);

91
92
93
94
    return 0;
  }


95
  // y = PC * x
96
  PetscErrorCode petscApplyFetiDirichletPrecon(PC pc, Vec x, Vec y)
97
  {
98
99
    double wtime = MPI::Wtime();

100
    // Get data for the preconditioner
101
102
    void *ctx;
    PCShellGetContext(pc, &ctx);
103
    FetiDirichletPreconData* data = static_cast<FetiDirichletPreconData*>(ctx);
104

105
    // Multiply with scaled Lagrange constraint matrix.
106
107
108
    MatMultTranspose(*(data->mat_lagrange_scaled), x, data->tmp_vec_b);


109
    // === Restriction of the B nodes to the boundary nodes. ===
110

111
112
113
114
    int nLocalB;
    int nLocalDuals;
    VecGetLocalSize(data->tmp_vec_b, &nLocalB);
    VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals);
115

116
117
118
119
    PetscScalar *local_b, *local_duals;
    VecGetArray(data->tmp_vec_b, &local_b);
    VecGetArray(data->tmp_vec_duals0, &local_duals);

120
121
122
    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_duals[it->second] = local_b[it->first];
123
124

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


128
    // === K_DD - K_DI inv(K_II) K_ID ===
129

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

132
    MatMult(*(data->mat_interior_duals), data->tmp_vec_duals0, data->tmp_vec_interior);
133
    KSPSolve(*(data->ksp_interior), data->tmp_vec_interior, data->tmp_vec_interior);
134
135
136
137
138
139
140
141
142
143
    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);

144
145
146
    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_b[it->first] = local_duals[it->second];
147
148
149
150
151
152
153
154

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

155
156
    FetiTimings::fetiPreconditioner += (MPI::Wtime() - wtime);

157
158
159
160
161
162
163
164
165
166
167
168
169
170
    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);
171
172


173
    // === Restriction of the B nodes to the boundary nodes. ===
174

175
176
177
178
    int nLocalB;
    int nLocalDuals;
    VecGetLocalSize(data->tmp_vec_b, &nLocalB);
    VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals);
179

180
    PetscScalar *local_b, *local_duals;
181
    VecGetArray(data->tmp_vec_b, &local_b);
182
    VecGetArray(data->tmp_vec_duals0, &local_duals);
183

Thomas Witkowski's avatar
Thomas Witkowski committed
184
185
186
    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_duals[it->second] = local_b[it->first];
187
188

    VecRestoreArray(data->tmp_vec_b, &local_b);
189
190
191
192
193
194
195
    VecRestoreArray(data->tmp_vec_duals0, &local_duals);


    // === K_DD ===

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

196

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

199
200
201
    VecGetArray(data->tmp_vec_b, &local_b);
    VecGetArray(data->tmp_vec_duals1, &local_duals);

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

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

209
210

    // Multiply with scaled Lagrange constraint matrix.
211
212
213
214
215
216
    MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b, y);

    return 0;
  }


217
218
  PetscSolverFeti::PetscSolverFeti()
    : PetscSolver(),
219
      schurPrimalSolver(0),
220
      multiLevelTest(false),
Thomas Witkowski's avatar
Thomas Witkowski committed
221
      subdomain(NULL),
222
223
      meshLevel(0),
      rStartInterior(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
224
225
      nGlobalOverallInterior(0),
      printTimings(false)
226
227
228
229
230
231
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

    string preconditionerName = "";
    Parameters::get("parallel->solver->precon", preconditionerName);
    if (preconditionerName == "" || preconditionerName == "none") {
232
      MSG("Create FETI-DP solver with no preconditioner!\n");
233
234
      fetiPreconditioner = FETI_NONE;
    } else if (preconditionerName == "dirichlet") {
235
      MSG("Create FETI-DP solver with Dirichlet preconditioner!\n");
236
237
      fetiPreconditioner = FETI_DIRICHLET;
    } else if (preconditionerName == "lumped") {
238
      MSG("Create FETI-DP solver with lumped preconditioner!\n");
239
240
      fetiPreconditioner = FETI_LUMPED;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
241
242
      ERROR_EXIT("Preconditioner \"%s\" not available!\n", 
		 preconditionerName.c_str());
243
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
244
245
246
247
248

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

    Parameters::get("parallel->multi level test", multiLevelTest);
251
252
    if (multiLevelTest)
      meshLevel = 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
253
254

    Parameters::get("parallel->print timings", printTimings);
255
256
257
  }


258
  void PetscSolverFeti::initialize(vector<const FiniteElemSpace*> feSpaces)
259
  {
260
261
    FUNCNAME("PetscSolverFeti::initialize()");

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

265
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();
266
    vector<const FiniteElemSpace*>& uniqueFe = meshDistributor->getFeSpaces();
267

Thomas Witkowski's avatar
Thomas Witkowski committed
268
269
    if (subdomain == NULL) {
      subdomain = new PetscSolverGlobalMatrix();
270

271
      if (meshLevel == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
272
	subdomain->setMeshDistributor(meshDistributor, 
273
				      mpiCommGlobal, mpiCommLocal);
274
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
275
	subdomain->setMeshDistributor(meshDistributor, 
276
277
				      levelData.getMpiComm(meshLevel - 1),
				      levelData.getMpiComm(meshLevel));
Thomas Witkowski's avatar
Thomas Witkowski committed
278
	subdomain->setLevel(meshLevel);
279
280
      }
    }
281

282
283
284
285
    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
286

287
288
    if (fullInterface != NULL)
      interfaceDofMap.init(levelData, fullInterface);
Thomas Witkowski's avatar
Thomas Witkowski committed
289

290
    if (fetiPreconditioner != FETI_NONE) {
291
292
293
      TEST_EXIT(meshLevel == 0)
	("Dirichlet preconditioner not yet implemented for multilevel FETI-DP\n");

294
      interiorDofMap.init(levelData, feSpaces, uniqueFe, false);
295
    }
296
297
298
299
300
301
302
  }


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

Thomas Witkowski's avatar
Thomas Witkowski committed
303
304
    double timeCounter = MPI::Wtime();

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

309
310
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

311
312
313
314
    primalDofMap.clear();
    dualDofMap.clear();
    lagrangeMap.clear();
    localDofMap.clear();
315
    if (fetiPreconditioner != FETI_NONE)
316
317
      interiorDofMap.clear();

318
319
    primalDofMap.setDofComm(meshDistributor->getDofComm());
    lagrangeMap.setDofComm(meshDistributor->getDofComm());
320

321
322
    primalDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    dualDofMap.setMpiComm(levelData.getMpiComm(0), 0);
323
    lagrangeMap.setMpiComm(levelData.getMpiComm(0), 0);
324
    localDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
325
326
    if (fetiPreconditioner != FETI_NONE)
      interiorDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
327

328
329
330
    if (meshLevel == 0)
      localDofMap.setDofComm(meshDistributor->getDofComm());
    else
331
332
      localDofMap.setDofComm(meshDistributor->getDofCommSd());

333
334
335
336
337
338
    if (fullInterface != NULL) {
      interfaceDofMap.clear();
      interfaceDofMap.setDofComm(meshDistributor->getDofComm());
      interfaceDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    }

339
340
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
341
    
342
343
      createPrimals(feSpace);  

344
345
346
      createDuals(feSpace);

      createInterfaceNodes(feSpace);
347

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
348
      createIndexB(feSpace);     
349
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
350
351
352
353

    primalDofMap.update();
    dualDofMap.update();
    localDofMap.update();
354
    if (fetiPreconditioner != FETI_NONE)
355
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
356

357
358
359
    if (fullInterface != NULL)
      interfaceDofMap.update();

360
361
362
363
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
      createLagrange(feSpace);
    }
364

365
366
367
    lagrangeMap.update();


368
369
370
371
372
373
374
375
376
377
378
379
    // === ===

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

380
      mpi::getDofNumbering(mpiCommGlobal, groupRowsInterior,
381
382
383
384
385
386
387
388
389
			   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);
    }

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

395
      MSG("  nRankPrimals = %d   nOverallPrimals = %d\n", 
396
397
	  primalDofMap[feSpace].nRankDofs, 
	  primalDofMap[feSpace].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
398
      
399
      MSG("  nRankDuals = %d  nOverallDuals = %d\n",
400
401
	  dualDofMap[feSpace].nRankDofs, 
	  dualDofMap[feSpace].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
402

403
      MSG("  nRankLagrange = %d  nOverallLagrange = %d\n",
404
405
	  lagrangeMap[feSpace].nRankDofs, 
	  lagrangeMap[feSpace].nOverallDofs);
406

407
      TEST_EXIT_DBG(localDofMap[feSpace].size() + primalDofMap[feSpace].size() == 
408
409
		    static_cast<unsigned int>(feSpace->getAdmin()->getUsedDofs()))
	("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
410
    }
411

Thomas Witkowski's avatar
Thomas Witkowski committed
412
413
    subdomain->setDofMapping(&localDofMap);
    subdomain->setCoarseSpaceDofMapping(&primalDofMap);
Thomas Witkowski's avatar
Thomas Witkowski committed
414
415

    if (printTimings) {
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
416
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
417
      timeCounter = MPI::Wtime() - timeCounter;
418
      MSG("FETI-DP timing 01: %.5f seconds (creation of basic data structures)\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
419
420
	  timeCounter);
    }
421
422
423
  }


424
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
425
  {
426
    FUNCNAME("PetscSolverFeti::createPrimals()");  
427

428
429
430
    if (feSpace == fullInterface)
      return;

431
432
433
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

Thomas Witkowski's avatar
Thomas Witkowski committed
434
    // Set of DOF indices that are considered to be primal variables.
435
    DofContainerSet& vertices = 
436
      meshDistributor->getBoundaryDofInfo(feSpace, meshLevel).geoDofs[VERTEX];
437
438

    DofIndexSet primals;
439
    for (DofContainerSet::iterator it = vertices.begin(); 
440
441
	 it != vertices.end(); ++it) {
      double e = 1e-8;
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
442
      if (meshLevel == 0) {
443
	primals.insert(**it);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
444
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
445
446
447
	WorldVector<double> c;
	feSpace->getMesh()->getDofIndexCoords(*it, feSpace, c);

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
448
449
450
451
452
453
454
	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);
	}
      }
455
    }
456

457
458
459
460

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

461
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
462
463
      if (meshDistributor->getDofMap()[feSpace].isRankDof(*it))
	primalDofMap[feSpace].insertRankDof(*it);
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
464
      else
465
  	primalDofMap[feSpace].insertNonRankDof(*it);
466
467
468
  }


469
  void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace)
470
471
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
472

473
474
475
    if (feSpace == fullInterface)
      return;

476
477
478
    // === Create global index of the dual nodes on each rank. ===

    DofContainer allBoundaryDofs;
479
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
480
481
482
483

    for (DofContainer::iterator it = allBoundaryDofs.begin();
	 it != allBoundaryDofs.end(); ++it)
      if (!isPrimal(feSpace, **it))
484
485
	if (meshLevel == 0) {
	  dualDofMap[feSpace].insertRankDof(**it);
486
487
488
489
 	} else {
 	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(**it))
 	    dualDofMap[feSpace].insertRankDof(**it);
 	}	  
490
491
492
  }

  
493
494
495
496
497
498
499
500
501
502
503
504
  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)
Thomas Witkowski's avatar
Thomas Witkowski committed
505
506
      if (meshDistributor->getDofMap()[feSpace].isRankDof(**it))
	interfaceDofMap[feSpace].insertRankDof(**it);
507
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
508
	interfaceDofMap[feSpace].insertNonRankDof(**it);
509
510
511
  }


512
513
514
515
  void PetscSolverFeti::createLagrange(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

516
517
518
    if (feSpace == fullInterface)
      return;

519
520
    boundaryDofRanks[feSpace].clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
521
522
523
    // 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
524
    map<int, std::set<DegreeOfFreedom> > sdRankDofs;
525

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
526
    if (meshLevel > 0) {
527
      StdMpi<vector<int> > stdMpi(mpiCommGlobal);
528
529
530
531
532
533
534
535
536

      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
537
	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))
538
	    subdomainRankDofs.push_back(1);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
539
	  else
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
	    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
555
556
557
558
559
560
	   !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());
    }
561

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
562
563
564
565
566
567
568
569
570
571
572
573
574
575
    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);

576
 	  if (meshLevel == 0 ||
Thomas Witkowski's avatar
Thomas Witkowski committed
577
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
578
	    boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());	  
579
580
581
582
	}
      }
    }

583
584
585
586

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

587
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
588

589
    for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), meshLevel, feSpace);
590
591
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
592
	if (!isPrimal(feSpace, it.getDofIndex()))
593
594
595
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
	    stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]);
596
597
598

    stdMpi.updateSendDataSize();

599
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
600
	 !it.end(); it.nextRank()) {
601
      bool recvFromRank = false;
602
      for (; !it.endDofIter(); it.nextDof()) {
603
	if (!isPrimal(feSpace, it.getDofIndex())) {
604
605
606
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
 	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
607
608
609
	    recvFromRank = true;
	    break;
	  }
610
	}
611
      }
612
613

      if (recvFromRank)
614
	stdMpi.recv(it.getRank());
615
    }
616

617
618
    stdMpi.startCommunication();

619
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
620
	 !it.end(); it.nextRank()) {
621
      int i = 0;
622
      for (; !it.endDofIter(); it.nextDof())
623
	if (!isPrimal(feSpace, it.getDofIndex()))
624
625
626
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
 	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex())))	    
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
627
628
	    boundaryDofRanks[feSpace][it.getDofIndex()] = 
	      stdMpi.getRecvData(it.getRank())[i++];
629
630
	  else
	    lagrangeMap[feSpace].insertNonRankDof(it.getDofIndex());
631
632
    }

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
633

634
635
636
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

637
    int nRankLagrange = 0;
638
    DofMap& dualMap = dualDofMap[feSpace].getMap();
639
    for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
640
641
      if (meshDistributor->getDofMap()[feSpace].isRankDof(it->first)) {
	lagrangeMap[feSpace].insertRankDof(it->first, nRankLagrange);
642
	int degree = boundaryDofRanks[feSpace][it->first].size();
643
	nRankLagrange += (degree * (degree - 1)) / 2;
Thomas Witkowski's avatar
Thomas Witkowski committed
644
      } else {
645
	lagrangeMap[feSpace].insertNonRankDof(it->first);
646
647
      }
    }
648
    lagrangeMap[feSpace].nRankDofs = nRankLagrange;
649
650
651
  }


652
  void PetscSolverFeti::createIndexB(const FiniteElemSpace *feSpace)
653
  {
654
    FUNCNAME("PetscSolverFeti::createIndexB()");
655

656
    DOFAdmin* admin = feSpace->getAdmin();
657
658
659
660

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

662
    int nLocalInterior = 0;    
663
    for (int i = 0; i < admin->getUsedSize(); i++) {
664
      if (admin->isDofFree(i) == false && 
665
	  isPrimal(feSpace, i) == false &&
666
667
	  isDual(feSpace, i) == false &&
	  isInterface(feSpace, i) == false) {
668
669
	if (meshLevel == 0) {
	  localDofMap[feSpace].insertRankDof(i, nLocalInterior);
670

671
672
673
674
675
676
677
678
679
	  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);
680
681
682

	  TEST_EXIT_DBG(fetiPreconditioner == FETI_NONE)
	    ("Not yet implemnted!\n");
683
	}
684
      }
685
    }
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
686
    
687
688
    // === And finally, add the global indicies of all dual nodes. ===

689
690
691
692
693
694
695
696
697
698
    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);
      }
699
700
701
  }


702
  void PetscSolverFeti::createMatLagrange(vector<const FiniteElemSpace*> &feSpaces)
703
704
705
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
706
707
    double wtime = MPI::Wtime();

708
709
    // === Create distributed matrix for Lagrange constraints. ===

710
711
712
713
714
    MatCreateAIJ(mpiCommGlobal,
		 lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		 lagrangeMap.getOverallDofs(), nGlobalOverallInterior,
		 2, PETSC_NULL, 2, PETSC_NULL,
		 &mat_lagrange);
715
    MatSetOption(mat_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
716

717
718
719
720
721
722
723
    // === 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
724
    for (unsigned int k = 0; k < feSpaces.size(); k++) {
725
      DofMap &dualMap = dualDofMap[feSpaces[k]].getMap();
726

727
      for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
728
	TEST_EXIT_DBG(boundaryDofRanks[feSpaces[k]].count(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
729
730
731
	  ("Should not happen!\n");
	
	// Global index of the first Lagrange constriant for this node.
732
	int index = lagrangeMap.getMatIndex(k, it->first);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
733

Thomas Witkowski's avatar
Thomas Witkowski committed
734
	// Copy set of all ranks that contain this dual node.
735
736
	vector<int> W(boundaryDofRanks[feSpaces[k]][it->first].begin(), 
		      boundaryDofRanks[feSpaces[k]][it->first].end());
Thomas Witkowski's avatar
Thomas Witkowski committed
737
738
	// Number of ranks that contain this dual node.
	int degree = W.size();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
739
740

	TEST_EXIT_DBG(degree > 1)("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
741
742
743
744
	
	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
745
746
	      int colIndex = 
		localDofMap.getMatIndex(k, it->first) + rStartInterior;
747

748
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
749
	      MatSetValue(mat_lagrange, index, colIndex, value, INSERT_VALUES);
750
	    }
Thomas Witkowski's avatar
Thomas Witkowski committed
751
	    index++;	      
752
753
754
755
756
757
758
	  }
	}
      }
    }

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

Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
760
761
    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
762
763
      MSG("FETI-DP timing 05: %.5f seconds (creation of lagrange constraint matrix)\n", 
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
764
    }
765
766
767
  }


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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
775
      schurPrimalData.subSolver = subdomain;
Thomas Witkowski's avatar
Thomas Witkowski committed
776
      
777
      VecCreateMPI(mpiCommGlobal, 
778
		   localDofMap.getRankDofs(), 
779
		   nGlobalOverallInterior,
Thomas Witkowski's avatar
Thomas Witkowski committed
780
		   &(schurPrimalData.tmp_vec_b));
781
      VecCreateMPI(mpiCommGlobal, 
782
783
		   primalDofMap.getRankDofs(), 
		   primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
784
785
		   &(schurPrimalData.tmp_vec_primal));

786
      MatCreateShell(mpiCommGlobal,
787
788
789
790
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getOverallDofs(), 
		     primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
791
792
793
794
795
		     &schurPrimalData, 
		     &mat_schur_primal);
      MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			   (void(*)(void))petscMultMatSchurPrimal);
      
796
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
797
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
798
799
      KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
      KSPSetType(ksp_schur_primal, KSPGMRES);
Thomas Witkowski's avatar
Thomas Witkowski committed
800
801
      KSPSetFromOptions(ksp_schur_primal);
    } else {
802
803
      MSG("Create direct schur primal solver!\n");

804
805
      double wtime = MPI::Wtime();

806
807
808
      TEST_EXIT_DBG(meshLevel == 0)
	("Does not support for multilevel, check usage of localDofMap.\n");

809
810
811
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
      int nRowsRankB = localDofMap.getRankDofs();
812

Thomas Witkowski's avatar
Thomas Witkowski committed
813
      Mat matBPi;
814
815
816
      MatCreateAIJ(mpiCommGlobal,
		   nRowsRankB, nRowsRankPrimal, 
		   nGlobalOverallInterior, nRowsOverallPrimal,
817
818
819
		   150, PETSC_NULL, 150, PETSC_NULL, &matBPi);
      MatSetOption(matBPi, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

820

821
822
823
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
824

825
826
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

827
      for (int i = 0; i < nRowsRankB; i++) {
828
	PetscInt row = localDofMap.getStartDofs() + i;
Thomas Witkowski's avatar
Thomas Witkowski committed
829
	MatGetRow(subdomain->getMatIntCoarse(), row, &nCols, &cols, &values);
830
831
832
833
834

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

Thomas Witkowski's avatar
Thomas Witkowski committed
835
	MatRestoreRow(subdomain->getMatIntCoarse(), row, &nCols, &cols, &values);
836
837
      }

838
      TEST_EXIT(static_cast<int>(mat_b_primal_cols.size()) == 
839
		primalDofMap.getLocalDofs())
840
	("Should not happen!\n");
841

842
843
844
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
845
	VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);
846
847
848
849
850
851
852
853

 	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
854
	subdomain->solve(tmpVec, tmpVec);
855

Thomas Witkowski's avatar
Thomas Witkowski committed
856
	PetscScalar *tmpValues;
857
	VecGetArray(tmpVec, &tmpValues);
858
	for (int i  = 0; i < nRowsRankB; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
859
	  MatSetValue(matBPi, 
860
		      localDofMap.getStartDofs() + i,
Thomas Witkowski's avatar
Thomas Witkowski committed
861
862
863
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
864
865
866
867
868
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
869
870
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
871

Thomas Witkowski's avatar
Thomas Witkowski committed
872
      MatDuplicate(subdomain->getMatCoarseCoarse(), MAT_COPY_VALUES, &mat_schur_primal);
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
873
874

      Mat matPrimal;
Thomas Witkowski's avatar
Thomas Witkowski committed
875
      MatMatMult(subdomain->getMatCoarseInt(), matBPi, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &matPrimal);
876
      MatAXPY(mat_schur_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);
877
878

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
879
      MatDestroy(&matBPi);
880

881
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
882
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal,
883
		      SAME_NONZERO_PATTERN);
884
885
886
887
888
889
      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);
890
      KSPSetFromOptions(ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
891

Thomas Witkowski's avatar
Thomas Witkowski committed
892
      if (printTimings) {
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
893
	MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
894
895
896
897
898
899
900
901
902
903
	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
904
	MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
905
906
907
	MSG("FETI-DP timing 07: %.5f seconds (factorization of primal schur matrix).\n",
	    MPI::Wtime() - wtime);
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
908
    }
909
910
911
912
913
914
915
  }


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

916
    if (schurPrimalSolver == 0) {
917
      schurPrimalData.subSolver = NULL;
918

919
920
921
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
    }
922
923
924
    
    MatDestroy(&mat_schur_primal);
    KSPDestroy(&ksp_schur_primal);